Merge branch 'feature-ami-face-area-intersect' into 'develop'

AMI improvements

See merge request !468

src/OpenFOAM/primitives/Vector2D/Vector2D.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,13 +116,18 @@ public:
         //- Access to the vector y component
         inline Cmpt& y();
 
-
         //- Normalise the vector by its magnitude
         inline Vector2D<Cmpt>& normalise();
 
-
         //- Perp dot product (dot product with perpendicular vector)
         inline scalar perp(const Vector2D<Cmpt>& b) const;
+
+        //- Return true if vector is within tol
+        inline bool isClose
+        (
+            const Vector2D<Cmpt>& b,
+            const scalar tol = 1e-10
+        ) const;
 };
 
 

src/OpenFOAM/primitives/Vector2D/Vector2DI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -126,7 +126,18 @@ operator&(const Vector2D<Cmpt>& v1, const Vector2D<Cmpt>& v2)
 template<class Cmpt>
 inline scalar Vector2D<Cmpt>::perp(const Vector2D<Cmpt>& b) const
 {
-    return x()*b.y()-y()*b.x();
+    return x()*b.y() - y()*b.x();
+}
+
+
+template<class Cmpt>
+inline bool Vector2D<Cmpt>::isClose
+(
+    const Vector2D<Cmpt>& b,
+    const scalar tol
+) const
+{
+    return (mag(x() - b.x()) < tol && mag(y() - b.y()) < tol);
 }
 
 

src/functionObjects/field/mapFields/mapFields.H

@@ -98,7 +98,6 @@ Usage
     \verbatim
       nearestFaceAMI
       faceAreaWeightAMI
-      partialFaceAreaWeightAMI
     \endverbatim
 
     The inherited entries are elaborated in:

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -1246,7 +1246,7 @@ void Foam::AMIInterpolation::write(Ostream& os) const
 
     if (reverseTarget_)
     {
-        os.writeEntry("flipNormals", reverseTarget_);
+        os.writeEntry("reverseTarget", reverseTarget_);
     }
 
     if (lowWeightCorrection_ > 0)

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -364,6 +364,9 @@ public:
             //- Access to the requireMatch flag
             inline bool setRequireMatch(const bool flag);
 
+            //- Return true if requireMatch and lowWeightCorrectionin active
+            inline bool mustMatchFaces() const;
+
             //- Access to the reverseTarget flag
             inline bool reverseTarget() const;
 

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,7 +73,7 @@ inline bool Foam::AMIInterpolation::distributed() const
 
 inline bool Foam::AMIInterpolation::requireMatch() const
 {
-    return requireMatch_ && lowWeightCorrection_ < 0;
+    return requireMatch_;
 }
 
 
@@ -84,6 +84,12 @@ inline bool Foam::AMIInterpolation::setRequireMatch(const bool flag)
 }
 
 
+inline bool Foam::AMIInterpolation::mustMatchFaces() const
+{
+    return requireMatch_ && !applyLowWeightCorrection();
+}
+
+
 inline bool Foam::AMIInterpolation::reverseTarget() const
 {
     return reverseTarget_;

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C

@@ -55,7 +55,7 @@ void Foam::advancingFrontAMI::checkPatches() const
     }
 
 
-    if (conformal())
+    if (requireMatch_)
     {
         const scalar maxBoundsError = 0.05;
 
@@ -345,6 +345,27 @@ void Foam::advancingFrontAMI::triangulatePatch
 }
 
 
+void Foam::advancingFrontAMI::nonConformalCorrection()
+{
+    if (!requireMatch_ && distributed())
+    {
+        scalarList newTgtMagSf(std::move(tgtMagSf_));
+
+        // Assign default sizes. Override selected values with calculated
+        // values. This is to support ACMI where some of the target faces
+        // are never used (so never get sent over and hence never assigned
+        // to)
+        tgtMagSf_ = tgtPatch0().magFaceAreas();
+
+        for (const labelList& smap : this->extendedTgtMapPtr_->subMap())
+        {
+            UIndirectList<scalar>(tgtMagSf_, smap) =
+                UIndirectList<scalar>(newTgtMagSf, smap);
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::advancingFrontAMI::advancingFrontAMI
@@ -421,7 +442,8 @@ bool Foam::advancingFrontAMI::calculate
 {
     if (AMIInterpolation::calculate(srcPatch, tgtPatch, surfPtr))
     {
-        // Create a representation of the target patch that covers the source patch
+        // Create a representation of the target patch that covers the source
+        // patch
         if (distributed())
         {
             createExtendedTgtPatch();
@@ -454,10 +476,4 @@ bool Foam::advancingFrontAMI::calculate
 }
 
 
-bool Foam::advancingFrontAMI::conformal() const
-{
-    return true;
-}
-
-
 // ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.H

@@ -199,6 +199,9 @@ protected:
                 List<scalar>& magSf
             ) const;
 
+            //- Correction for non-conformal interpolations, e.g. for ACMI
+            virtual void nonConformalCorrection();
+
 
 public:
 
@@ -249,10 +252,6 @@ public:
         //- Labels of faces that are not overlapped by any target faces
         //  Note: this should be empty for correct functioning
         inline const labelList& srcNonOverlap() const;
-
-        //- Flag to indicate that interpolation patches are conformal, i.e.
-        //- should fully cover each other
-        virtual bool conformal() const;
 };
 
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C

@@ -115,6 +115,9 @@ void Foam::faceAreaWeightAMI::calcAddressing
     // Reset starting seed
     label startSeedi = 0;
 
+    // Should all faces be matched?
+    const bool mustMatch = mustMatchFaces();
+
     bool continueWalk = true;
     DynamicList<label> nonOverlapFaces;
     do
@@ -154,7 +157,7 @@ void Foam::faceAreaWeightAMI::calcAddressing
             mapFlag,
             seedFaces,
             visitedFaces,
-            requireMatch_ && (lowWeightCorrection_ < 0)
+            mustMatch
             // pass in nonOverlapFaces for failed tree search?
         );
     } while (continueWalk);
@@ -675,7 +678,7 @@ bool Foam::faceAreaWeightAMI::calculate
         }
 
         // Check for badly covered faces
-        if (restartUncoveredSourceFace_) // && requireMatch_???
+        if (restartUncoveredSourceFace_) //  && mustMatchFaces())
         {
             restartUncoveredSourceFace
             (
@@ -775,7 +778,9 @@ bool Foam::faceAreaWeightAMI::calculate
     }
 
     // Convert the weights from areas to normalised values
-    normaliseWeights(conformal(), true);
+    normaliseWeights(requireMatch_, true);
+
+    nonConformalCorrection();
 
     upToDate_ = true;
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.H

@@ -30,6 +30,8 @@ Class
 Description
     Face area weighted Arbitrary Mesh Interface (AMI) method
 
+    Searching is performed using an advancing front.
+
 SourceFiles
     faceAreaWeightAMI.C
 
@@ -61,7 +63,6 @@ private:
         const bool restartUncoveredSourceFace_;
 
 
-
 protected:
 
     // Protected Member Functions
@@ -69,15 +70,6 @@ protected:
         //- No copy assignment
         void operator=(const faceAreaWeightAMI&) = delete;
 
-        //- Initialise the geometry
-        void initGeom
-        (
-            const primitivePatch& srcPatch,
-            const primitivePatch& tgtPatch,
-            const globalIndex& globalTgtFaces,
-            labelList& extendedTgtFaceIDs
-        );
-
 
         // Marching front
 

src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "faceAreaWeightAMI2D.H"
+#include "profiling.H"
+#include "OBJstream.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triangle2D.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(faceAreaWeightAMI2D, 0);
+    addToRunTimeSelectionTable(AMIInterpolation, faceAreaWeightAMI2D, dict);
+    addToRunTimeSelectionTable
+    (
+        AMIInterpolation,
+        faceAreaWeightAMI2D,
+        component
+    );
+}
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+void Foam::faceAreaWeightAMI2D::writeNoMatch
+(
+    const label srcFacei,
+    const labelList& tgtFaceCandidates,
+    const boundBox& srcFaceBb
+) const
+{
+    Info<< "NO MATCH for source face " << srcFacei << endl;
+    {
+        const auto& src = this->srcPatch();
+        const auto& tgt = this->tgtPatch(); // might be the extended patch!
+
+        OFstream os("no_match_" + Foam::name(srcFacei) + ".obj");
+
+        const pointField& srcPoints = src.points();
+        const pointField& tgtPoints = tgt.points();
+
+        label np = 0;
+
+        // Write source face
+        const face& srcF = src[srcFacei];
+        string faceStr = "f";
+        for (const label pointi : srcF)
+        {
+            const point& p = srcPoints[pointi];
+            os  << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
+            ++np;
+            faceStr += " " + Foam::name(np);
+        }
+        os << faceStr.c_str() << " " << (np - srcF.size() + 1) << nl;
+
+        // Write target faces as lines
+        for (const label tgtFacei : tgtFaceCandidates)
+        {
+            const face& tgtF = tgt[tgtFacei];
+            forAll(tgtF, pointi)
+            {
+                const point& p = tgtPoints[tgtF[pointi]];
+                os  << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
+                ++np;
+                if (pointi)
+                {
+                    os << "l " << np-1 << " " << np << nl;
+                }
+            }
+            os << "l " << (np - tgtF.size() + 1) << " " << np << nl;
+        }
+    }
+
+    {
+        OFstream os("no_match_" + Foam::name(srcFacei) + "_bb.obj");
+
+        const pointField points(srcFaceBb.points());
+        for (const point& p : points)
+        {
+            os  << "v " << p.x() << " " << p.y() << " " << p.z() << endl;
+        }
+        os  << "l 1 2" << nl;
+        os  << "l 2 4" << nl;
+        os  << "l 4 3" << nl;
+        os  << "l 3 1" << nl;
+        os  << "l 5 6" << nl;
+        os  << "l 6 8" << nl;
+        os  << "l 8 7" << nl;
+        os  << "l 7 5" << nl;
+        os  << "l 5 1" << nl;
+        os  << "l 6 2" << nl;
+        os  << "l 8 4" << nl;
+        os  << "l 7 3" << nl;
+    }
+}
+
+
+void Foam::faceAreaWeightAMI2D::storeInterArea
+(
+    const label srcFacei,
+    const label tgtFacei,
+    DynamicList<label>& srcAddr,
+    DynamicList<scalar>& srcWght,
+    DynamicList<vector>& srcCtr,
+    DynamicList<label>& tgtAddr,
+    DynamicList<scalar>& tgtWght
+) const
+{
+    addProfiling(ami, "faceAreaWeightAMI2D::calcInterArea");
+
+    // Quick reject if either face has zero area
+    if
+    (
+        (srcMagSf_[srcFacei] < ROOTVSMALL)
+     || (tgtMagSf_[tgtFacei] < ROOTVSMALL)
+    )
+    {
+        return;
+    }
+
+    const auto& srcPatch = this->srcPatch();
+    const auto& tgtPatch = this->tgtPatch();
+
+    const pointField& srcPoints = srcPatch.points();
+    const pointField& tgtPoints = tgtPatch.points();
+
+    const auto& srcTris = srcTris_[srcFacei];
+    const auto& tgtTris = tgtTris_[tgtFacei];
+
+    const auto& srcFaceNormals = srcPatch.faceNormals();
+
+    //triangle2D::debug = 1;
+
+    scalar area = 0;
+    vector centroid = Zero;
+
+    for (const auto& tris : srcTris)
+    {
+        const vector& origin = srcPoints[tris[0]];
+        const vector p10(srcPoints[tris[1]] - origin);
+        const vector p20(srcPoints[tris[2]] - origin);
+        const vector axis1(p10/(mag(p10) + ROOTVSMALL));
+        const vector axis2(srcFaceNormals[srcFacei]^axis1);
+
+        triangle2D s
+        (
+            vector2D(0, 0),
+            vector2D(axis1&p10, axis2&p10),
+            vector2D(axis1&p20, axis2&p20)
+        );
+
+        for (const auto& trit : tgtTris)
+        {
+            // Triangle t has opposite orientation wrt triangle s
+            triangle2D t
+            (
+                tgtPoints[trit[0]] - origin,
+                tgtPoints[trit[2]] - origin,
+                tgtPoints[trit[1]] - origin,
+                axis1,
+                axis2
+            );
+
+            scalar da = 0;
+            vector2D c(Zero);
+
+            if (t.snapClosePoints(s) == 3)
+            {
+                c = s.centre();
+                da = mag(s.area());
+            }
+            else
+            {
+                s.interArea(t, c, da);
+            }
+
+            area += da;
+            centroid += da*(origin + c.x()*axis1 + c.y()*axis2);
+        }
+    }
+
+    //triangle2D::debug = 0;
+
+    if (area/srcMagSf_[srcFacei] > triangle2D::relTol)
+    {
+        centroid /= area + ROOTVSMALL;
+
+        srcAddr.append(tgtFacei);
+        srcWght.append(area);
+        srcCtr.append(centroid);
+
+        tgtAddr.append(srcFacei);
+        tgtWght.append(area);
+    }
+}
+
+
+Foam::labelList Foam::faceAreaWeightAMI2D::overlappingTgtFaces
+(
+    const AABBTree<face>& tree,
+    const List<boundBox>& tgtFaceBbs,
+    const boundBox& srcFaceBb
+) const
+{
+    labelHashSet faceIds(16);
+
+    const auto& treeBb = tree.boundBoxes();
+    const auto& treeAddr = tree.addressing();
+
+    forAll(treeBb, boxi)
+    {
+        const auto& tbb = treeBb[boxi];
+
+        if (srcFaceBb.overlaps(tbb))
+        {
+            const auto& boxAddr = treeAddr[boxi];
+
+            for (const auto& tgtFacei : boxAddr)
+            {
+                if (srcFaceBb.overlaps(tgtFaceBbs[tgtFacei]))
+                {
+                    faceIds.insert(tgtFacei);
+                }
+            }
+        }
+    }
+
+    return faceIds.toc();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceAreaWeightAMI2D::faceAreaWeightAMI2D
+(
+    const dictionary& dict,
+    const bool reverseTarget
+)
+:
+    advancingFrontAMI(dict, reverseTarget),
+    Cbb_(dict.getCheckOrDefault<scalar>("Cbb", 0.1, scalarMinMax::ge(SMALL)))
+{}
+
+
+Foam::faceAreaWeightAMI2D::faceAreaWeightAMI2D
+(
+    const bool requireMatch,
+    const bool reverseTarget,
+    const scalar lowWeightCorrection,
+    const faceAreaIntersect::triangulationMode triMode,
+    const bool restartUncoveredSourceFace
+)
+:
+    advancingFrontAMI
+    (
+        requireMatch,
+        reverseTarget,
+        lowWeightCorrection,
+        triMode
+    ),
+    Cbb_(0.1)
+{}
+
+
+Foam::faceAreaWeightAMI2D::faceAreaWeightAMI2D
+(
+    const faceAreaWeightAMI2D& ami
+)
+:
+    advancingFrontAMI(ami),
+    Cbb_(ami.Cbb_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::faceAreaWeightAMI2D::calculate
+(
+    const primitivePatch& srcPatch,
+    const primitivePatch& tgtPatch,
+    const autoPtr<searchableSurface>& surfPtr
+)
+{
+    if (upToDate_)