zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

patchToPatch: Improve robustness of non-intersection methods

Browse Source 
The nearest, matching and inverseDistance methods are now based on a
shared "nearby" method. This method creates, for each face, a local
stencil of opposing faces for which the bounding spheres overlap. This
has proven far more robust on cases with both conformal and
non-conformal interfaces.
This commit is contained in:
Will Bainbridge
2022-09-23 16:50:47 +01:00
parent 97a00ff520
commit 0203618a91
21 changed files with 1416 additions and 348 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/test/boundSphere/Make/files Normal file
View File

@ -0,0 +1,3 @@
Test-boundSphere.C
EXE = $(FOAM_USER_APPBIN)/Test-boundSphere

5
applications/test/boundSphere/Make/options Normal file
View File

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/fileFormats/lnInclude
EXE_LIBS = \
-lfileFormats

224
applications/test/boundSphere/Test-boundSphere.C Normal file
View File

@ -0,0 +1,224 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "boundSphere.H"
#include "clock.H"
#include "cpuTime.H"
#include "OBJstream.H"
#include "triFace.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::validArgs.append("number of points");
argList::addOption
(
"nTests",
"value",
"number of tests - default is 1"
);
argList::addOption
(
"seed",
"value",
"random number generator seed - default is clock time"
);
Foam::argList args(argc, argv);
Info<< nl;
const label nPs = args.argRead<label>(1);
const label nTests = args.optionLookupOrDefault<label>("nTests", 1);
// Initialise the random number generator
label seed;
if (args.optionFound("seed"))
{
seed = args.optionRead<label>("seed");
}
else
{
seed = clock::getTime();
Info<< "Seeding random number generator with value " << seed
<< nl << endl;
}
// Do tests
List<point> ps;
Tuple2<point, scalar> sphere;
cpuTime time;
for (label testi = 0; testi < nTests; ++ testi)
{
Random rndGen(seed + testi);
// Get random points
ps.resize(nPs);
forAll(ps, i)
{
ps[i] = point::uniform(1);
while (magSqr(ps[i]) > 1)
{
ps[i] = 2*(rndGen.sample01<point>() - point::uniform(0.5));
}
}
// Determine the bounds
sphere = boundSphere(ps, rndGen);
// Check
forAll(ps, i)
{
if
(
magSqr(ps[i] - sphere.first())
> (1 + rootSmall)*sqr(sphere.second())
)
{
FatalErrorInFunction
<< "Sphere does not bound point " << i << '/' << nPs << " ("
<< magSqr(ps[i] - sphere.first())/sqr(sphere.second())
<< ")" << exit(FatalError);
}
}
if (nTests != 1)
{
Info<< "\rSuccessfully computed " << testi + 1 << " bound spheres";
}
}
// Report
if (nTests == 1)
{
Info<< "Bounding sphere with centre = " << sphere.first()
<< ", radius = " << sphere.second() << ", calculated in "
<< time.cpuTimeIncrement() << " s" << nl << endl;
// Write the points
OBJstream osPs(args[0] + "_points.obj");
Info<< "Writing points to " << osPs.name() << endl;
forAll(ps, i)
{
osPs.write(ps[i]);
}
// Write the sphere
OBJstream osSphere(args[0] + "_sphere.obj");
Info<< "Writing sphere to " << osSphere.name() << endl;
using namespace constant::mathematical;
static const label nRings = 20;
static const label nRing1Points = 4;
DynamicList<point> spherePs;
DynamicList<face> sphereTris;
auto appendP = [&](const scalar x, const scalar y, const scalar z)
{
spherePs.append(sphere.first() + sphere.second()*point(x, y, z));
spherePs.append(sphere.first() + sphere.second()*point(x, y, -z));
};
auto appendTri = [&](const label a, const label b, const label c)
{
sphereTris.append(triFace(2*a, 2*b, 2*c));
sphereTris.append(triFace(2*a + 1, 2*b + 1, 2*c + 1));
};
appendP(0, 0, 1);
label nPoints0 = 0, nPoints = 1;
for (label ringi = 1; ringi < nRings; ++ ringi)
{
const scalar rXY = Foam::sin(pi/2*ringi/(nRings - 1));
const label ringPN0 = max((ringi - 1)*nRing1Points, 1);
const label ringPN = ringi*nRing1Points;
label ringPi0 = 0, ringPi = 0;
for (label i = 0; i < nRing1Points; ++ i)
{
const scalar theta0 = 2*pi*i/nRing1Points;
const scalar theta1 = 2*pi*(i + 1)/nRing1Points;
for (label j = 0; j < ringi; ++ j)
{
const scalar f = scalar(j)/ringi;
const scalar theta = (1 - f)*theta0 + f*theta1;
const scalar x = rXY*Foam::cos(theta);
const scalar y = rXY*Foam::sin(theta);
const scalar z = Foam::sqrt(max(1 - sqr(x) - sqr(y), 0));
appendP(x, y, z);
appendTri
(
nPoints0 + (ringPi0 % ringPN0),
nPoints + ringPi,
nPoints + ((ringPi + 1) % ringPN)
);
if (ringi > 1 && j < ringi - 1)
{
appendTri
(
nPoints0 + (ringPi0 % ringPN0),
nPoints + ((ringPi + 1) % ringPN),
nPoints0 + ((ringPi0 + 1) % ringPN0)
);
}
if (j < ringi - 1) ringPi0 ++;
ringPi ++;
}
}
nPoints0 = nPoints;
nPoints += ringPi;
}
osSphere.write(faceList(sphereTris), pointField(spherePs), false);
}
else
{
Info<< nl << nl << "Bounding spheres calculated in "
<< time.cpuTimeIncrement() << " s" << endl;
}
Info<< nl << "End" << nl << endl;
return 0;
}
// ************************************************************************* //

54
applications/test/patchToPatch/Test-patchToPatch.C
View File

@ -39,18 +39,66 @@ int main(int argc, char *argv[])
argList::validArgs.append("source");
argList::validArgs.append("target");
argList::validArgs.append("method");
argList::validArgs.append("reverse");
argList::addOption("sourceCase", "dir", "The case with the source patch");
#include "setRootCase.H"
#include "createTime.H"
#include "createPolyMesh.H"
const polyPatch& srcPatch = mesh.boundaryMesh()[args[1]];
const polyPatch& tgtPatch = mesh.boundaryMesh()[args[2]];
// Optionally read a different mesh for the source
autoPtr<Time> srcRunTimePtr;
autoPtr<polyMesh> srcMeshPtr;
if (args.optionFound("sourceCase"))
{
const string tgtCase = getEnv("FOAM_CASE");
const string tgtCaseName = getEnv("FOAM_CASENAME");
fileName sourceCase = args["sourceCase"];
sourceCase.clean();
const fileName sourceCaseName =
Pstream::parRun()
? fileName(sourceCase.name())/args.caseName().name()
: fileName(sourceCase.name());
setEnv("FOAM_CASE", sourceCase, true);
setEnv("FOAM_CASENAME", sourceCase.name(), true);
srcRunTimePtr.set
(
new Time(sourceCase.path(), sourceCaseName)
);
setEnv("FOAM_CASE", tgtCase, true);
setEnv("FOAM_CASENAME", tgtCaseName, true);
srcMeshPtr.set
(
new polyMesh
(
Foam::IOobject
(
Foam::polyMesh::defaultRegion,
srcRunTimePtr->timeName(),
srcRunTimePtr(),
Foam::IOobject::MUST_READ
)
)
);
}
const polyMesh& srcMesh = srcMeshPtr.valid() ? srcMeshPtr() : mesh;
const polyMesh& tgtMesh = mesh;
const polyPatch& srcPatch = srcMesh.boundaryMesh()[args[1]];
const polyPatch& tgtPatch = tgtMesh.boundaryMesh()[args[2]];
const word& method = args[3];
const bool reverse = args.argRead<bool>(4);
cpuTime time;
patchToPatch::New(method, false)->update
patchToPatch::New(method, reverse)->update
(
srcPatch,
srcPatch.pointNormals(),

6
src/OpenFOAM/Make/files
View File

@ -715,16 +715,12 @@ $(interpolationWeights)/interpolationWeights/interpolationWeights.C
$(interpolationWeights)/linearInterpolationWeights/linearInterpolationWeights.C
$(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
algorithms/indexedOctree/indexedOctreeName.C
algorithms/indexedOctree/treeDataCell.C
algorithms/indexedOctree/volumeType.C
algorithms/polygonTriangulate/polygonTriangulate.C
algorithms/dynamicIndexedOctree/dynamicIndexedOctreeName.C
algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C
algorithms/polygonTriangulate/polygonTriangulate.C
meshes/data/data.C
meshes/Residuals/residuals.C

99
src/OpenFOAM/algorithms/boundSphere/boundSphere.H Normal file
View File

@ -0,0 +1,99 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Functions for constructing bounding spheres of lists of points
SourceFiles
boundSphere.C
boundSphereTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef boundSphere_H
#define boundSphere_H
#include "point.H"
#include "Random.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//- Return whether or not the given sphere is valid
inline bool isValidBoundSphere(const Tuple2<point, scalar>& sphere);
//- Compute a sphere of four points or less where every point intersects the
// sphere's surface
template<class PointField>
Tuple2<point, scalar> intersectBoundSphere
(
const PointField& ps,
const FixedList<label, 4>& pis,
const label nPs
);
//- Compute a bounding sphere of four points or less
template<class PointField>
Tuple2<point, scalar> trivialBoundSphere
(
const PointField& ps,
const FixedList<label, 4>& pis,
const label nPs
);
//- Compute a bounding sphere for an arbitrary number of points recursively
// using Weizl's algorithm
template<class PointField>
Tuple2<point, scalar> weizlBoundSphere
(
const PointField& ps,
List<label>& pis,
const label nPs,
FixedList<label, 4>& boundaryPis,
const label nBoundaryPs
);
//- Compute a bounding sphere for an arbitrary number of points, and given an
// engine with which to randomise Weizl's algorithm
template<class PointField>
Tuple2<point, scalar> boundSphere(const PointField& ps, Random& rndGen);
//- Compute a bounding sphere for an arbitrary number of points
template<class PointField>
Tuple2<point, scalar> boundSphere(const PointField& ps);
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "boundSphereTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

301
src/OpenFOAM/algorithms/boundSphere/boundSphereTemplates.C Normal file
View File

@ -0,0 +1,301 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "boundSphere.H"
#include "labelPair.H"
#include "triPointRef.H"
#include "tetPointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
inline bool Foam::isValidBoundSphere(const Tuple2<point, scalar>& sphere)
{
return sphere.second() >= 0;
}
template<class PointField>
Foam::Tuple2<Foam::point, Foam::scalar>
Foam::intersectBoundSphere
(
const PointField& ps,
const FixedList<label, 4>& pis,
const label nPs
)
{
switch (nPs)
{
case 0:
return Tuple2<point, scalar>(point::uniform(NaN), -vGreat);
case 1:
return Tuple2<point, scalar>(ps[pis[0]], 0);
case 2:
// Return a mid point-centred sphere
return Tuple2<point, scalar>
(
(ps[pis[0]] + ps[pis[1]])/2,
mag(ps[pis[0]] - ps[pis[1]])/2
);
case 3:
{
// Return the sphere that intersects the triangle
return
triPointRef
(
ps[pis[0]],
ps[pis[1]],
ps[pis[2]]
).circumCircle();
}
case 4:
{
// Return the sphere that intersects the tetrahedron
return
tetPointRef
(
ps[pis[0]],
ps[pis[1]],
ps[pis[2]],
ps[pis[3]]
).circumSphere();
}
default:
FatalErrorInFunction
<< "Cannot compute the intersect bounding sphere of more than "
<< "four points" << exit(FatalError);
return Tuple2<point, scalar>(point::uniform(NaN), NaN);
}
}
template<class PointField>
Foam::Tuple2<Foam::point, Foam::scalar>
Foam::trivialBoundSphere
(
const PointField& ps,
const FixedList<label, 4>& pis,
const label nPs
)
{
static const scalar tol = 1 + 2*sqrt(small*rootSmall);
// Search for spheres that intersect sub-sets of the points that bound all
// the other points
switch (nPs)
{
case 0:
case 1:
case 2:
break;
case 3:
{
// Test the spheres that intersect the edges
for (label i = 0; i < 3; ++ i)
{
const point& p0 = ps[pis[i]];
const point& p1 = ps[pis[(i + 1) % 3]];
const point& pOpp = ps[pis[(i + 2) % 3]];
const point c = (p0 + p1)/2;
const scalar rSqr = magSqr(p0 - p1)/4;
if (magSqr(pOpp - c) <= tol*rSqr)
{
return Tuple2<point, scalar>(c, Foam::sqrt(rSqr));
}
}
break;
}
case 4:
{
// Test the spheres that intersect the edges
static const FixedList<labelPair, 6> p01s =
{{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}};
static const FixedList<labelPair, 6> pOpp01s =
{{2, 3}, {1, 3}, {1, 2}, {0, 3}, {0, 2}, {0, 1}};
for (label i = 0; i < 6; ++ i)
{
const point& p0 = ps[pis[p01s[i].first()]];
const point& p1 = ps[pis[p01s[i].second()]];
const point& pOpp0 = ps[pis[pOpp01s[i].first()]];
const point& pOpp1 = ps[pis[pOpp01s[i].second()]];
const point c = (p0 + p1)/2;
const scalar rSqr = magSqr(p0 - p1)/4;
if
(
magSqr(pOpp0 - c) <= tol*rSqr
&& magSqr(pOpp1 - c) <= tol*rSqr
)
{
return Tuple2<point, scalar>(c, Foam::sqrt(rSqr));
}
}
// Test the spheres that intersect the triangles
point minC = point::uniform(NaN);
scalar minR = vGreat;
for (label i = 0; i < 4; ++ i)
{
const triPointRef tri
(
ps[pis[i]],
ps[pis[(i + 1) % 4]],
ps[pis[(i + 2) % 4]]
);
const point& pOpp = ps[pis[(i + 3) % 4]];
const Tuple2<point, scalar> circ = tri.circumCircle();
const point& c = circ.first();
const scalar r = circ.second();
if (magSqr(pOpp - c) <= tol*sqr(r) && r < minR)
{
minC = c;
minR = r;
}
}
if (minR != vGreat)
{
return Tuple2<point, scalar>(minC, minR);
}
break;
}
default:
FatalErrorInFunction
<< "Cannot compute the trivial bounding sphere of more than "
<< "four points" << exit(FatalError);
return Tuple2<point, scalar>(point::uniform(NaN), NaN);
}
// Return the sphere that intersects the points
return intersectBoundSphere(ps, pis, nPs);
}
template<class PointField>
Foam::Tuple2<Foam::point, Foam::scalar> Foam::weizlBoundSphere
(
const PointField& ps,
List<label>& pis,
const label nPs,
FixedList<label, 4>& boundaryPis,
const label nBoundaryPs
)
{
static const scalar tol = 1 + 2*sqrt(small*rootSmall);
Tuple2<point, scalar> sphere;
if (nBoundaryPs != 0)
{
sphere = intersectBoundSphere(ps, boundaryPis, nBoundaryPs);
}
if (nBoundaryPs == 4)
{
return sphere;
}
for (label i = 0; i < nPs; ++ i)
{
if
(
(nBoundaryPs == 0 && i == 0)
|| (magSqr(ps[pis[i]] - sphere.first()) > tol*sqr(sphere.second()))
)
{
boundaryPis[nBoundaryPs] = pis[i];
sphere = weizlBoundSphere(ps, pis, i, boundaryPis, nBoundaryPs + 1);
// Move the limiting point to the start of the list so that the
// sphere grows as quickly as possible in the recursive calls
Swap(pis[0], pis[i]);
}
}
return sphere;
}
template<class PointField>
Foam::Tuple2<Foam::point, Foam::scalar>
Foam::boundSphere(const PointField& ps, Random& rndGen)
{
if (ps.size() <= 4)
{
static const FixedList<label, 4> pis({0, 1, 2, 3});
return trivialBoundSphere(ps, pis, ps.size());
}
else
{
labelList pis(identity(ps.size()));
forAll(pis, i)
{
Swap(pis[i], pis[rndGen.sampleAB<label>(i, ps.size())]);
}
FixedList<label, 4> boundaryPis({-1, -1, -1, -1});
return weizlBoundSphere(ps, pis, ps.size(), boundaryPis, 0);
}
}
template<class PointField>
Foam::Tuple2<Foam::point, Foam::scalar>
Foam::boundSphere(const PointField& ps)
{
if (ps.size() <= 4)
{
static const FixedList<label, 4> pis({0, 1, 2, 3});
return trivialBoundSphere(ps, pis, ps.size());
}
else
{
labelList pis(identity(ps.size()));
Random rndGen(0);
forAll(pis, i)
{
Swap(pis[i], pis[rndGen.sampleAB<label>(i, ps.size())]);
}
FixedList<label, 4> boundaryPis({-1, -1, -1, -1});
return weizlBoundSphere(ps, pis, ps.size(), boundaryPis, 0);
}
}
// ************************************************************************* //

1
src/meshTools/Make/files
View File

@ -220,6 +220,7 @@ patchToPatchTools/patchToPatchTools.C
patchToPatch/patchToPatch/patchToPatch.C
patchToPatch/patchToPatch/patchToPatchParallelOps.C
patchToPatch/matching/matchingPatchToPatch.C
patchToPatch/nearby/nearbyPatchToPatch.C
patchToPatch/nearest/nearestPatchToPatch.C
patchToPatch/inverseDistance/inverseDistancePatchToPatch.C
patchToPatch/intersection/intersectionPatchToPatch.C

12
src/meshTools/meshToMesh/meshToMesh.C
View File

@ -265,11 +265,11 @@ void Foam::meshToMesh::constructNoCuttingPatches
srcToTgtPatchIDs_.transfer(srcToTgtPatchIDs);
}
// calculate cell addressing and weights
calculateCellToCells(methodName);
// calculate patch addressing and weights
calculatePatchToPatches(methodName);
// calculate cell addressing and weights
calculateCellToCells(methodName);
}
@ -295,12 +295,12 @@ void Foam::meshToMesh::constructFromCuttingPatches
);
}
// calculate cell addressing and weights
calculateCellToCells(methodName);
// calculate patch addressing and weights
calculatePatchToPatches(methodName);
// calculate cell addressing and weights
calculateCellToCells(methodName);
// set IDs of cutting patches on target mesh
tgtCuttingPatchIDs_.setSize(tgtCuttingPatches.size());
forAll(tgtCuttingPatchIDs_, i)

15
src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C
View File

@ -542,13 +542,22 @@ void Foam::patchToPatches::intersection::initialise
}
Foam::labelList Foam::patchToPatches::intersection::trimLocalTgt
Foam::labelList Foam::patchToPatches::intersection::finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
)
{
const labelList newToOldLocalTgtFace =
patchToPatch::trimLocalTgt(localTgtPatch);
patchToPatch::finaliseLocal
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch
);
tgtCouples_ = List<DynamicList<couple>>(tgtCouples_, newToOldLocalTgtFace);

12
src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.H
View File

@ -328,11 +328,15 @@ private:
const primitiveOldTimePatch& tgtPatch
);
//- Trim the local target patch so that only parts that actually
// intersect the source remain
virtual labelList trimLocalTgt
//- Finalise the intersection locally. Trims the local target patch
// so that only parts that actually intersect the source remain.
// Returns new-to-old map for trimming target-associated data.
virtual labelList finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
);
//- Send data that resulted from an intersection between the source

261
src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.C
View File

@ -39,31 +39,15 @@ namespace patchToPatches
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
Foam::treeBoundBox Foam::patchToPatches::inverseDistance::srcBox
bool Foam::patchToPatches::inverseDistance::rayHitsFace
(
const face& srcFace,
const pointField& srcPoints,
const vectorField& srcPointNormals
) const
{
const treeBoundBox bb(srcPoints, srcFace);
const point c = bb.midpoint();
const scalar l = bb.maxDim();
return treeBoundBox(c - l*vector::one, c + l*vector::one);
}
bool Foam::patchToPatches::inverseDistance::inside
(
const face& f,
const pointField& ps,
const point& p,
const vector& r
) const
const vector& r,
const face& f,
const pointField& ps
)
{
using namespace constant::mathematical;
@ -75,100 +59,23 @@ bool Foam::patchToPatches::inverseDistance::inside
{
const vector& a = T & (ps[f[i]] - p);
const vector& b = T & (ps[f[f.fcIndex(i)]] - p);
const scalar magAB = sqrt(magSqr(a)*magSqr(b));
angle +=
(reverse() ? +1 : -1)
*sign(r & (a ^ b))
*acos(min(max(-1, (a & b)/magAB), +1));
const scalar meanMagSqrAB = (magSqr(a) + magSqr(b))/2;
const scalar geometricMeanMagSqrAB = sqrt(magSqr(a)*magSqr(b));
// This indicates that we have hit a point to within round off error
if (geometricMeanMagSqrAB < small*meanMagSqrAB) return true;
angle -=
sign(r & (a ^ b))
*acos(min(max(-1, (a & b)/geometricMeanMagSqrAB), +1));
}
return pi < angle && angle < 3*pi;
}
bool Foam::patchToPatches::inverseDistance::intersectFaces
(
const primitivePatch& patch,
const primitivePatch& otherPatch,
const label facei,
const label otherFacei,
DynamicList<label>& faceOtherFaces,
DynamicList<scalar>& faceWeights
) const
{
const face& f = otherPatch[otherFacei];
const pointField& ps = otherPatch.points();
const point& p = patch.faceCentres()[facei];
const vector& r = patch.faceNormals()[facei];
bool intersectsOther = inside(f, ps, p, r);
if (!intersectsOther)
{
forAll(otherPatch.faceFaces()[otherFacei], otherFaceFacei)
{
const label otherFacej =
otherPatch.faceFaces()[otherFacei][otherFaceFacei];
const face& g = otherPatch[otherFacej];
if (inside(g, ps, p, r))
{
intersectsOther = true;
break;
}
}
}
if (intersectsOther)
{
faceOtherFaces.append(otherFacei);
faceWeights.append
(
1/max(mag(p - otherPatch.faceCentres()[otherFacei]), vSmall)
);
}
return intersectsOther;
}
bool Foam::patchToPatches::inverseDistance::intersectFaces
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch,
const label srcFacei,
const label tgtFacei
)
{
const bool srcCouples =
intersectFaces
(
srcPatch,
tgtPatch,
srcFacei,
tgtFacei,
srcLocalTgtFaces_[srcFacei],
srcWeights_[srcFacei]
);
const bool tgtCouples =
intersectFaces
(
tgtPatch,
srcPatch,
tgtFacei,
srcFacei,
tgtLocalSrcFaces_[tgtFacei],
tgtWeights_[tgtFacei]
);
return srcCouples || tgtCouples;
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::patchToPatches::inverseDistance::initialise
(
@ -178,7 +85,7 @@ void Foam::patchToPatches::inverseDistance::initialise
const primitiveOldTimePatch& tgtPatch
)
{
patchToPatch::initialise
nearby::initialise
(
srcPatch,
srcPointNormals,
@ -187,17 +94,122 @@ void Foam::patchToPatches::inverseDistance::initialise
);
srcWeights_.resize(srcPatch.size());
forAll(srcWeights_, i)
{
srcWeights_[i].clear();
}
tgtWeights_.resize(tgtPatch.size());
forAll(tgtWeights_, i)
{
tgtWeights_[i].clear();
}
}
Foam::labelList Foam::patchToPatches::inverseDistance::trimLocalTgt
void Foam::patchToPatches::inverseDistance::generateWeights
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const primitiveOldTimePatch& tgtPatch
)
{
auto generate = []
(
const primitiveOldTimePatch& patch,
const primitiveOldTimePatch& otherPatch,
const bool reverse,
List<DynamicList<label>>& otherFaces,
List<DynamicList<scalar>>& weights
)
{
forAll(otherFaces, facei)
{
if (otherFaces[facei].empty()) continue;
label otherFacei = -1;
// Find the other face that "contains" this face's centre
forAll(otherFaces[facei], i)
{
if
(
rayHitsFace
(
patch.faceCentres()[facei],
(reverse ? -1 : +1)*patch.faceNormals()[facei],
otherPatch[otherFaces[facei][i]],
otherPatch.points()
)
)
{
otherFacei = otherFaces[facei][i];
break;
}
}
const point& c = patch.faceCentres()[facei];
// If the above failed, find the closest
if (otherFacei == -1)
{
scalar minDistSqr = vGreat;
forAll(otherFaces[facei], i)
{
const point& otherC =
otherPatch.faceCentres()[otherFaces[facei][i]];
const scalar distSqr = magSqr(c - otherC);
if (distSqr < minDistSqr)
{
minDistSqr = distSqr;
otherFacei = otherFaces[facei][i];
}
}
}
// Remove all faces
otherFaces[facei].clear();
// Add the found face and all its neighbours
const point& otherC = otherPatch.faceCentres()[otherFacei];
otherFaces[facei].append(otherFacei);
weights[facei].append(1/(mag(c - otherC) + rootVSmall));
forAll(otherPatch.faceFaces()[otherFacei], i)
{
const label otherFacej = otherPatch.faceFaces()[otherFacei][i];
const point& otherC = otherPatch.faceCentres()[otherFacej];
otherFaces[facei].append(otherFacej);
weights[facei].append(1/(mag(c - otherC) + rootVSmall));
}
}
};
generate(srcPatch, tgtPatch, reverse_, srcLocalTgtFaces_, srcWeights_);
generate(tgtPatch, srcPatch, reverse_, tgtLocalSrcFaces_, tgtWeights_);
}
Foam::labelList Foam::patchToPatches::inverseDistance::finaliseLocal
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
)
{
// Transfer weight addressing into actual addressing
generateWeights(srcPatch, tgtPatch);
const labelList newToOldLocalTgtFace =
patchToPatch::trimLocalTgt(localTgtPatch);
nearby::finaliseLocal
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch
);
tgtWeights_ = List<DynamicList<scalar>>(tgtWeights_, newToOldLocalTgtFace);
@ -210,8 +222,10 @@ void Foam::patchToPatches::inverseDistance::rDistributeTgt
const primitiveOldTimePatch& tgtPatch
)
{
patchToPatch::rDistributeTgt(tgtPatch);
// Let the base class reverse distribute the addressing
nearby::rDistributeTgt(tgtPatch);
// Reverse distribute the weights
patchToPatchTools::rDistributeListList
(
tgtPatch.size(),
@ -240,20 +254,25 @@ Foam::label Foam::patchToPatches::inverseDistance::finalise
tgtToSrc
);
// Transfer weight addressing into actual addressing (if not done in the
// finaliseLocal method above)
if (isSingleProcess())
{
generateWeights(srcPatch, tgtPatch);
}
// Normalise the weights
forAll(srcWeights_, srcFacei)
{
const scalar w = sum(srcWeights_[srcFacei]);
forAll(srcWeights_[srcFacei], i)
{
srcWeights_[srcFacei][i] /= max(w, vSmall);
}
}
forAll(tgtWeights_, tgtFacei)
{
const scalar w = sum(tgtWeights_[tgtFacei]);
forAll(tgtWeights_[tgtFacei], i)
{
tgtWeights_[tgtFacei][i] /= max(w, vSmall);
@ -279,8 +298,10 @@ Foam::label Foam::patchToPatches::inverseDistance::finalise
return result;
};
Info<< "Number of source faces by number of target connections = "
Info<< indent
<< "Number of source faces by number of target connections = "
<< histogram(srcLocalTgtFaces_) << nl
<< indent
<< "Number of target faces by number of source connections = "
<< histogram(tgtLocalSrcFaces_) << endl;
}
@ -307,7 +328,7 @@ Foam::patchToPatches::inverseDistance::tgtWeights() const
Foam::patchToPatches::inverseDistance::inverseDistance(const bool reverse)
:
patchToPatch(reverse),
nearby(reverse),
srcWeights_(),
tgtWeights_()
{}

70
src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.H
View File

@ -38,7 +38,7 @@ SourceFiles
#ifndef inverseDistancePatchToPatch_H
#define inverseDistancePatchToPatch_H
#include "patchToPatch.H"
#include "nearbyPatchToPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -53,7 +53,7 @@ namespace patchToPatches
class inverseDistance
:
public patchToPatch
public nearby
{
// Private Member Data
@ -64,47 +64,21 @@ class inverseDistance
List<DynamicList<scalar>> tgtWeights_;
// Private Member Functions
// Private Static Member Functions
//- Get the bound box for a source face
virtual treeBoundBox srcBox
//- Return whether a ray starting at point `p` with direction `r` hits
// a face `f` with points `ps`
static bool rayHitsFace
(
const face& srcFace,
const pointField& srcPoints,
const vectorField& srcPointNormals
) const;
//- Return whether or not the face contains a point
bool inside
(
const face& f,
const pointField& ps,
const point& p,
const vector& r
) const;
//- Intersect two faces
virtual bool intersectFaces
(
const primitivePatch& patch,
const primitivePatch& otherPatch,
const label facei,
const label otherFacei,
DynamicList<label>& faceOtherFaces,
DynamicList<scalar>& faceWeights
) const;
//- Intersect two faces
virtual bool intersectFaces
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch,
const label srcFacei,
const label tgtFacei
const vector& r,
const face& f,
const pointField& ps
);
// Private Member Functions
//- Initialisation
virtual void initialise
(
@ -114,11 +88,23 @@ class inverseDistance
const primitiveOldTimePatch& tgtPatch
);
//- Trim the local target patch so that only parts that actually
// intersect the source remain
virtual labelList trimLocalTgt
//- Convert the addressing generated by the base class into weights
// (and addressing) appropriate for inverse distance interpolation
void generateWeights
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const primitiveOldTimePatch& tgtPatch
);
//- Finalise the intersection locally. Trims the local target patch
// so that only parts that actually intersect the source remain.
// Returns new-to-old map for trimming target-associated data.
virtual labelList finaliseLocal
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
);
//- Send data that resulted from an intersection between the source

145
src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.C Normal file
View File

@ -0,0 +1,145 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "nearbyPatchToPatch.H"
#include "boundSphere.H"
#include "OFstream.H"
#include "OBJstream.H"
#include "vtkWritePolyData.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace patchToPatches
{
defineTypeNameAndDebug(nearby, 0);
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::treeBoundBox Foam::patchToPatches::nearby::srcBox
(
const face& srcFace,
const pointField& srcPoints,
const vectorField& srcPointNormals
) const
{
const treeBoundBox bb(srcPoints, srcFace);
const point c = bb.midpoint();
const scalar l = bb.maxDim();
return treeBoundBox(c - l*vector::one, c + l*vector::one);
}
bool Foam::patchToPatches::nearby::intersectFaces
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch,
const label srcFacei,
const label tgtFacei
)
{
const point& srcC = srcSpheres_[srcFacei].first();
const scalar srcR = srcSpheres_[srcFacei].second();
const point& tgtC = tgtSpheres_[tgtFacei].first();
const scalar tgtR = tgtSpheres_[tgtFacei].second();
if (magSqr(srcC - tgtC) < sqr(srcR + tgtR))
{
srcLocalTgtFaces_[srcFacei].append(tgtFacei);
tgtLocalSrcFaces_[tgtFacei].append(srcFacei);
return true;
}
else
{
return false;
}
}
void Foam::patchToPatches::nearby::initialise
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
)
{
patchToPatch::initialise
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch
);
srcSpheres_.resize(srcPatch.size());
forAll(srcPatch, srcFacei)
{
srcSpheres_[srcFacei] =
boundSphere
(
UIndirectList<point>(srcPatch.points(), srcPatch[srcFacei])
);
}
tgtSpheres_.resize(tgtPatch.size());
forAll(tgtPatch, tgtFacei)
{
tgtSpheres_[tgtFacei] =
boundSphere
(
UIndirectList<point>(tgtPatch.points(), tgtPatch[tgtFacei])
);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::patchToPatches::nearby::nearby(const bool reverse)
:
patchToPatch(reverse),
srcSpheres_(),
tgtSpheres_()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::patchToPatches::nearby::~nearby()
{}
// ************************************************************************* //

124
src/meshTools/patchToPatch/nearby/nearbyPatchToPatch.H Normal file
View File

@ -0,0 +1,124 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2021-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::patchToPatches::nearby
Description
Class to generate patchToPatch coupling geometry. Couples a face to the
single nearby opposite face only.
SourceFiles
nearby.C
\*---------------------------------------------------------------------------*/
#ifndef nearbyPatchToPatch_H
#define nearbyPatchToPatch_H
#include "patchToPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace patchToPatches
{
/*---------------------------------------------------------------------------*\
Class nearby Declaration
\*---------------------------------------------------------------------------*/
class nearby
:
public patchToPatch
{
protected:
// Private Member Data
//- For each source face, the bounding sphere
List<Tuple2<point, scalar>> srcSpheres_;
//- For each target face, the bounding sphere
List<Tuple2<point, scalar>> tgtSpheres_;
// Private Member Functions
//- Get the bound box for a source face
virtual treeBoundBox srcBox
(
const face& srcFace,
const pointField& srcPoints,
const vectorField& srcPointNormals
) const;
//- Intersect two faces
virtual bool intersectFaces
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch,
const label srcFacei,
const label tgtFacei
);
//- Initialisation
virtual void initialise
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
);
public:
//- Runtime type information
TypeName("nearby");
// Constructors
//- Construct from components
nearby(const bool reverse);
//- Destructor
~nearby();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace patchToPatches
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

325
src/meshTools/patchToPatch/nearest/nearestPatchToPatch.C
View File

@ -25,6 +25,11 @@ License
#include "nearestPatchToPatch.H"
#include "addToRunTimeSelectionTable.H"
#include "boundSphere.H"
#include "OFstream.H"
#include "OBJstream.H"
#include "vtkWritePolyData.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -40,105 +45,6 @@ namespace patchToPatches
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::treeBoundBox Foam::patchToPatches::nearest::srcBox
(
const face& srcFace,
const pointField& srcPoints,
const vectorField& srcPointNormals
) const
{
const treeBoundBox bb(srcPoints, srcFace);
const point c = bb.midpoint();
const scalar l = bb.maxDim();
return treeBoundBox(c - l*vector::one, c + l*vector::one);
}
bool Foam::patchToPatches::nearest::intersectFaces
(
const primitivePatch& patch,
const primitivePatch& otherPatch,
const label facei,
const label otherFacei,
DynamicList<label>& faceOtherFaces,
scalar& faceDistance
) const
{
auto closest = [&patch,&otherPatch]
(
const label facei,
const label otherFacei
)
{
const point& c = patch.faceCentres()[facei];
const point& otherC = otherPatch.faceCentres()[otherFacei];
const scalar distSqr = magSqr(c - otherC);
forAll(otherPatch.faceEdges()[otherFacei], otherFaceEdgei)
{
const label otherEdgei =
otherPatch.faceEdges()[otherFacei][otherFaceEdgei];
point otherNbrC;
if (otherPatch.edgeFaces()[otherEdgei].size() == 2)
{
const label facej =
otherPatch.edgeFaces()[otherEdgei]
[otherPatch.edgeFaces()[otherEdgei][0] == otherFacei];
otherNbrC = otherPatch.faceCentres()[facej];
}
else
{
const edge& e = otherPatch.edges()[otherEdgei];
const point& p = otherPatch.localPoints()[e[0]];
const vector dp = e.vec(otherPatch.localPoints());
const vector n = otherPatch.faceNormals()[otherFacei] ^ dp;
otherNbrC = p + ((tensor::I - 2*sqr(n)) & (otherC - p));
}
if (magSqr(c - otherNbrC) < distSqr)
{
return false;
}
}
return true;
};
if (closest(facei, otherFacei))
{
const point& c = patch.faceCentres()[facei];
const point& otherC = otherPatch.faceCentres()[otherFacei];
const scalar distSqr = magSqr(c - otherC);
if (faceOtherFaces.empty() || faceDistance > distSqr)
{
faceOtherFaces.clear();
faceOtherFaces.append(otherFacei);
faceDistance = distSqr;
}
return true;
}
const labelList& otherFaceFaces = otherPatch.faceFaces()[otherFacei];
forAll(otherFaceFaces, otherFaceFacei)
{
if (closest(facei, otherFaceFaces[otherFaceFacei]))
{
return true;
}
}
return false;
}
bool Foam::patchToPatches::nearest::intersectFaces
(
const primitiveOldTimePatch& srcPatch,
@ -149,30 +55,52 @@ bool Foam::patchToPatches::nearest::intersectFaces
const label tgtFacei
)
{
const bool srcCouples =
intersectFaces
if
(
nearby::intersectFaces
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch,
srcFacei,
tgtFacei,
srcLocalTgtFaces_[srcFacei],
srcDistances_[srcFacei]
);
tgtFacei
)
)
{
const scalar dSqr =
magSqr
(
srcPatch.faceCentres()[srcFacei]
- tgtPatch.faceCentres()[tgtFacei]
);
const bool tgtCouples =
intersectFaces
(
tgtPatch,
srcPatch,
tgtFacei,
srcFacei,
tgtLocalSrcFaces_[tgtFacei],
tgtDistances_[tgtFacei]
);
if (dSqr < srcDistances_[srcFacei])
{
srcDistances_[srcFacei] = dSqr;
Swap
(
srcLocalTgtFaces_[srcFacei].first(),
srcLocalTgtFaces_[srcFacei].last()
);
}
return srcCouples || tgtCouples;
if (dSqr < tgtDistances_[tgtFacei])
{
tgtDistances_[tgtFacei] = dSqr;
Swap
(
tgtLocalSrcFaces_[tgtFacei].first(),
tgtLocalSrcFaces_[tgtFacei].last()
);
}
return true;
}
else
{
return false;
}
}
@ -184,7 +112,7 @@ void Foam::patchToPatches::nearest::initialise
const primitiveOldTimePatch& tgtPatch
)
{
patchToPatch::initialise
nearby::initialise
(
srcPatch,
srcPointNormals,
@ -200,13 +128,22 @@ void Foam::patchToPatches::nearest::initialise
}
Foam::labelList Foam::patchToPatches::nearest::trimLocalTgt
Foam::labelList Foam::patchToPatches::nearest::finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
)
{
const labelList newToOldLocalTgtFace =
patchToPatch::trimLocalTgt(localTgtPatch);
nearby::finaliseLocal
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch
);
tgtDistances_ = List<scalar>(tgtDistances_, newToOldLocalTgtFace);
@ -219,6 +156,22 @@ void Foam::patchToPatches::nearest::rDistributeTgt
const primitiveOldTimePatch& tgtPatch
)
{
// Keep only the closest opposing face
forAll(srcLocalTgtFaces_, srcFacei)
{
srcLocalTgtFaces_[srcFacei].resize
(
min(srcLocalTgtFaces_[srcFacei].size(), 1)
);
}
forAll(tgtLocalSrcFaces_, tgtFacei)
{
tgtLocalSrcFaces_[tgtFacei].resize
(
min(tgtLocalSrcFaces_[tgtFacei].size(), 1)
);
}
// Create a list-list of distances to match the addressing
List<List<scalar>> tgtDistances(tgtLocalSrcFaces_.size());
forAll(tgtLocalSrcFaces_, tgtFacei)
@ -230,7 +183,7 @@ void Foam::patchToPatches::nearest::rDistributeTgt
}
// Let the base class reverse distribute the addressing
patchToPatch::rDistributeTgt(tgtPatch);
nearby::rDistributeTgt(tgtPatch);
// Reverse distribute the distances
patchToPatchTools::rDistributeListList
@ -259,6 +212,132 @@ void Foam::patchToPatches::nearest::rDistributeTgt
}
Foam::label Foam::patchToPatches::nearest::finalise
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch,
const transformer& tgtToSrc
)
{
// Keep only the closest opposing face
forAll(srcLocalTgtFaces_, srcFacei)
{
srcLocalTgtFaces_[srcFacei].resize
(
min(srcLocalTgtFaces_[srcFacei].size(), 1)
);
}
forAll(tgtLocalSrcFaces_, tgtFacei)
{
tgtLocalSrcFaces_[tgtFacei].resize
(
min(tgtLocalSrcFaces_[tgtFacei].size(), 1)
);
}
const label nCouples =
nearby::finalise
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch,
tgtToSrc
);
if (debug)
{
auto countNonEmpty = [](const List<DynamicList<label>>& ll)
{
label result = 0;
forAll(ll, i)
{
result += !ll[i].empty();
}
return returnReduce(result, sumOp<label>());
};
Info<< indent
<< "Coupled " << countNonEmpty(srcLocalTgtFaces_)
<< "/" << returnReduce(srcLocalTgtFaces_.size(), sumOp<label>())
<< " source faces and " << countNonEmpty(tgtLocalSrcFaces_)
<< "/" << returnReduce(tgtLocalSrcFaces_.size(), sumOp<label>())
<< " target faces" << endl;
}
if (debug && !Pstream::parRun())
{
auto writeConnections = []
(
const primitivePatch& patch,
const primitivePatch& otherPatch,
const bool isSrc,
const List<DynamicList<label>>& faceLocalOtherFaces
)
{
const word name =
typeName + "_" + (isSrc ? "src" : "tgt") + "Connections";
OBJstream obj(name + ".obj");
forAll(faceLocalOtherFaces, facei)
{
const point& p = patch.faceCentres()[facei];
forAll(faceLocalOtherFaces[facei], i)
{
const label otherFacei = faceLocalOtherFaces[facei][i];
const point& q = otherPatch.faceCentres()[otherFacei];
obj.write(linePointRef(p, q));
}
}
};
writeConnections(srcPatch, tgtPatch, true, srcLocalTgtFaces_);
writeConnections(tgtPatch, srcPatch, false, tgtLocalSrcFaces_);
auto writeNotConnected = []
(
const primitivePatch& patch,
const List<DynamicList<label>>& faceLocalOtherFaces,
const bool isSrc
)
{
DynamicList<label> unconnected;
forAll(faceLocalOtherFaces, facei)
{
if (faceLocalOtherFaces[facei].empty())
{
unconnected.append(facei);
}
}
const word name =
typeName + "_" + (isSrc ? "src" : "tgt") + "NotConnected";
vtkWritePolyData::write
(
name + ".vtk",
name,
false,
patch.localPoints(),
labelList(),
labelListList(),
UIndirectList<face>(patch.localFaces(), unconnected)
);
};
writeNotConnected(srcPatch, srcLocalTgtFaces_, true);
writeNotConnected(tgtPatch, tgtLocalSrcFaces_, false);
}
return nCouples;
}
Foam::tmpNrc<Foam::List<Foam::DynamicList<Foam::scalar>>>
Foam::patchToPatches::nearest::srcWeights() const
{
@ -307,7 +386,9 @@ Foam::patchToPatches::nearest::tgtWeights() const
Foam::patchToPatches::nearest::nearest(const bool reverse)
:
patchToPatch(reverse)
nearby(reverse),
srcDistances_(),
tgtDistances_()
{}

45
src/meshTools/patchToPatch/nearest/nearestPatchToPatch.H
View File

@ -36,7 +36,7 @@ SourceFiles
#ifndef nearestPatchToPatch_H
#define nearestPatchToPatch_H
#include "patchToPatch.H"
#include "nearbyPatchToPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -51,7 +51,7 @@ namespace patchToPatches
class nearest
:
public patchToPatch
public nearby
{
protected:
@ -66,25 +66,6 @@ protected:
// Private Member Functions
//- Get the bound box for a source face
virtual treeBoundBox srcBox
(
const face& srcFace,
const pointField& srcPoints,
const vectorField& srcPointNormals
) const;
//- Intersect two faces
bool intersectFaces
(
const primitivePatch& patch,
const primitivePatch& otherPatch,
const label facei,
const label otherFacei,
DynamicList<label>& faceOtherFaces,
scalar& faceDistance
) const;
//- Intersect two faces
virtual bool intersectFaces
(
@ -105,11 +86,15 @@ protected:
const primitiveOldTimePatch& tgtPatch
);
//- Trim the local target patch so that only parts that actually
// intersect the source remain
virtual labelList trimLocalTgt
//- Finalise the intersection locally. Trims the local target patch
// so that only parts that actually intersect the source remain.
// Returns new-to-old map for trimming target-associated data.
virtual labelList finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
);
//- Send data that resulted from an intersection between the source
@ -117,6 +102,16 @@ protected:
// original target processes
virtual void rDistributeTgt(const primitiveOldTimePatch& tgtPatch);
//- Finalise the intersection
virtual label finalise
(
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch,
const transformer& tgtToSrc
);
//- For each source face, the coupled target weights
virtual tmpNrc<List<DynamicList<scalar>>> srcWeights() const;

17
src/meshTools/patchToPatch/patchToPatch/patchToPatch.C
View File

@ -657,13 +657,16 @@ void Foam::patchToPatch::initialise
}
Foam::labelList Foam::patchToPatch::trimLocalTgt
Foam::labelList Foam::patchToPatch::finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
)
{
// Determine which local target faces are actually used
boolList oldLocalTgtFaceIsUsed(localTgtPatch.size(), false);
boolList oldLocalTgtFaceIsUsed(tgtPatch.size(), false);
forAll(srcLocalTgtFaces_, srcFacei)
{
forAll(srcLocalTgtFaces_[srcFacei], i)
@ -961,7 +964,13 @@ void Foam::patchToPatch::update
}
// Trim the local target patch
trimLocalTgt(localTTgtPatch);
finaliseLocal
(
srcPatch,
srcPointNormals,
srcPointNormals,
localTTgtPatch
);
// Distribute the source patch
srcMapPtr_ =

12
src/meshTools/patchToPatch/patchToPatch/patchToPatch.H
View File

@ -215,11 +215,15 @@ protected:
const primitiveOldTimePatch& tgtPatch
);
//- Trim the local target patch so that only parts that actually
// intersect the source remain
virtual labelList trimLocalTgt
//- Finalise the intersection locally. Trims the local target patch
// so that only parts that actually intersect the source remain.
// Returns new-to-old map for trimming target-associated data.
virtual labelList finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
);
//- Distribute the source patch so that any information needed by

21
src/meshTools/patchToPatch/rays/raysPatchToPatch.C
View File

@ -86,21 +86,30 @@ bool Foam::patchToPatches::rays::intersectFaces
}
Foam::labelList Foam::patchToPatches::rays::trimLocalTgt
Foam::labelList Foam::patchToPatches::rays::finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
)
{
const labelList newToOldLocalTgtFace =
patchToPatch::trimLocalTgt(localTgtPatch);
patchToPatch::finaliseLocal
(
srcPatch,
srcPointNormals,
srcPointNormals0,
tgtPatch
);
localTgtPatchPtr_.reset
(
new PrimitiveOldTimePatch<faceList, pointField>
(
faceList(localTgtPatch, newToOldLocalTgtFace),
localTgtPatch.points(),
localTgtPatch.points0()
faceList(tgtPatch, newToOldLocalTgtFace),
tgtPatch.points(),
tgtPatch.points0()
)
);

12
src/meshTools/patchToPatch/rays/raysPatchToPatch.H
View File

@ -92,11 +92,15 @@ class rays
const label tgtFacei
);
//- Subset the local target patch so that only parts that actually
// intersect the source remain
virtual labelList trimLocalTgt
//- Finalise the intersection locally. Trims the local target patch
// so that only parts that actually intersect the source remain.
// Returns new-to-old map for trimming target-associated data.
virtual labelList finaliseLocal
(
const primitiveOldTimePatch& localTgtPatch
const primitiveOldTimePatch& srcPatch,
const vectorField& srcPointNormals,
const vectorField& srcPointNormals0,
const primitiveOldTimePatch& tgtPatch
);
//- Distribute the source patch so that any information needed by