ENH: AMI: Consistency between overlap/normalisation areas

The patch magSf calculation has been changed so that it uses the same
triangulation as the overlap algorithm. This improves consistency and
means that for exactly conforming patches (typically before any mesh
motion) the weights do not require normalisation.
This commit is contained in:
Will Bainbridge
2017-10-31 08:59:21 +00:00
committed by Andrew Heather
parent 293c0c3014
commit d76923d851
5 changed files with 111 additions and 70 deletions

View File

@ -119,6 +119,48 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::wordTointerpolationMethod
}
template<class SourcePatch, class TargetPatch>
template<class Patch>
Foam::tmp<Foam::scalarField>
Foam::AMIInterpolation<SourcePatch, TargetPatch>::patchMagSf
(
const Patch& patch,
const faceAreaIntersect::triangulationMode triMode
)
{
tmp<scalarField> tResult(new scalarField(patch.size(), Zero));
scalarField& result = tResult.ref();
const pointField& patchPoints = patch.localPoints();
faceList patchFaceTris;
forAll(result, patchFacei)
{
faceAreaIntersect::triangulate
(
patch.localFaces()[patchFacei],
patchPoints,
triMode,
patchFaceTris
);
forAll(patchFaceTris, i)
{
result[patchFacei] +=
triPointRef
(
patchPoints[patchFaceTris[i][0]],
patchPoints[patchFaceTris[i][1]],
patchPoints[patchFaceTris[i][2]]
).mag();
}
}
return tResult;
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class SourcePatch, class TargetPatch>
@ -884,16 +926,8 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
<< endl;
// Calculate face areas
srcMagSf_.setSize(srcPatch.size());
forAll(srcMagSf_, facei)
{
srcMagSf_[facei] = srcPatch[facei].mag(srcPatch.points());
}
tgtMagSf_.setSize(tgtPatch.size());
forAll(tgtMagSf_, facei)
{
tgtMagSf_[facei] = tgtPatch[facei].mag(tgtPatch.points());
}
srcMagSf_ = patchMagSf(srcPatch, triMode_);
tgtMagSf_ = patchMagSf(tgtPatch, triMode_);
// Calculate if patches present on multiple processors
singlePatchProc_ = calcDistribution(srcPatch, tgtPatch);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@ -104,6 +104,14 @@ public:
const word& method
);
//- Calculate the patch face magnitudes for the given tri-mode
template<class Patch>
static tmp<scalarField> patchMagSf
(
const Patch& patch,
const faceAreaIntersect::triangulationMode triMode
);
private:

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -323,6 +323,49 @@ Foam::faceAreaIntersect::faceAreaIntersect
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::faceAreaIntersect::triangulate
(
const face& f,
const pointField& points,
const triangulationMode& triMode,
faceList& faceTris
)
{
faceTris.resize(f.nTriangles());
switch (triMode)
{
case tmFan:
{
for (label i = 0; i < f.nTriangles(); ++ i)
{
faceTris[i] = face(3);
faceTris[i][0] = f[0];
faceTris[i][1] = f[i + 1];
faceTris[i][2] = f[i + 2];
}
break;
}
case tmMesh:
{
const label nFaceTris = f.nTriangles();
label nFaceTris1 = 0;
const label nFaceTris2 = f.triangles(points, nFaceTris1, faceTris);
if (nFaceTris != nFaceTris1 || nFaceTris != nFaceTris2)
{
FatalErrorInFunction
<< "The numbers of reported triangles in the face do not "
<< "match that generated by the triangulation"
<< exit(FatalError);
}
}
}
}
Foam::scalar Foam::faceAreaIntersect::calc
(
const face& faceA,
@ -332,32 +375,9 @@ Foam::scalar Foam::faceAreaIntersect::calc
)
{
// split faces into triangles
DynamicList<face> trisA;
DynamicList<face> trisB;
switch (triMode)
{
case tmFan:
{
triangleFan(faceA, trisA);
triangleFan(faceB, trisB);
break;
}
case tmMesh:
{
faceA.triangles(pointsA_, trisA);
faceB.triangles(pointsB_, trisB);
break;
}
default:
{
FatalErrorInFunction
<< "Unknown triangulation mode enumeration"
<< abort(FatalError);
}
}
faceList trisA, trisB;
triangulate(faceA, pointsA_, triMode, trisA);
triangulate(faceB, pointsB_, triMode, trisB);
// intersect triangles
scalar totalArea = 0.0;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -105,13 +105,6 @@ private:
FixedList<triPoints, 10>& tris
) const;
//- Decompose face into triangle fan
inline void triangleFan
(
const face& f,
DynamicList<face>& faces
) const;
//- Return point of intersection between plane and triangle edge
inline point planeIntersection
(
@ -162,6 +155,15 @@ public:
//- Fraction of local length scale to use as intersection tolerance
inline static scalar& tolerance();
//- Triangulate a face using the given triangulation mode
static void triangulate
(
const face& f,
const pointField& points,
const triangulationMode& triMode,
faceList& faceTris
);
//- Return area of intersection of faceA with faceB
scalar calc
(

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
@ -67,29 +67,6 @@ inline Foam::triPoints Foam::faceAreaIntersect::getTriPoints
}
inline void Foam::faceAreaIntersect::triangleFan
(
const face& f,
DynamicList<face>& faces
) const
{
if (f.size() > 2)
{
const label v0 = 0;
labelList indices(3);
for (label i = 1; i < f.size() - 1; i++)
{
indices[0] = f[v0];
indices[1] = f[i];
indices[2] = f[i + 1];
faces.append(face(indices));
}
}
}
inline Foam::point Foam::faceAreaIntersect::planeIntersection
(
const FixedList<scalar, 3>& d,