From db8751c521e3e55df08edac925ac48d8317a59ae Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 31 Oct 2017 08:59:21 +0000 Subject: [PATCH] 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. --- .../AMIInterpolation/AMIInterpolation.C | 54 +++++++++++--- .../AMIInterpolation/AMIInterpolation.H | 10 ++- .../faceAreaIntersect/faceAreaIntersect.C | 74 ++++++++++++------- .../faceAreaIntersect/faceAreaIntersect.H | 18 +++-- .../faceAreaIntersect/faceAreaIntersectI.H | 25 +------ 5 files changed, 111 insertions(+), 70 deletions(-) diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C index 292eb6aeb..825de3f33 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C @@ -124,6 +124,48 @@ Foam::AMIInterpolation::wordTointerpolationMethod } +template +template +Foam::tmp +Foam::AMIInterpolation::patchMagSf +( + const Patch& patch, + const faceAreaIntersect::triangulationMode triMode +) +{ + tmp 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 @@ -849,16 +891,8 @@ void Foam::AMIInterpolation::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); diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H index a2a144970..4840f82fb 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H @@ -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 | ------------------------------------------------------------------------------- License @@ -104,6 +104,14 @@ public: const word& method ); + //- Calculate the patch face magnitudes for the given tri-mode + template + static tmp patchMagSf + ( + const Patch& patch, + const faceAreaIntersect::triangulationMode triMode + ); + private: diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C index 617a85a0a..69c415233 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C @@ -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 @@ -325,6 +325,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, @@ -334,32 +377,9 @@ Foam::scalar Foam::faceAreaIntersect::calc ) { // split faces into triangles - DynamicList trisA; - DynamicList 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; diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H index a8dd73915..5b4f12a10 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -106,13 +106,6 @@ private: FixedList& tris ) const; - //- Decompose face into triangle fan - inline void triangleFan - ( - const face& f, - DynamicList& faces - ) const; - //- Return point of intersection between plane and triangle edge inline point planeIntersection ( @@ -163,6 +156,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 ( diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H index 87fab781e..b7bfb3836 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H @@ -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 | ------------------------------------------------------------------------------- License @@ -67,29 +67,6 @@ inline Foam::faceAreaIntersect::triPoints Foam::faceAreaIntersect::getTriPoints } -inline void Foam::faceAreaIntersect::triangleFan -( - const face& f, - DynamicList& 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& d,