diff --git a/src/OpenFOAM/primitives/Vector2D/Vector2D.H b/src/OpenFOAM/primitives/Vector2D/Vector2D.H
index 9b67b6187c372da948280af0d211bb784e3b4689..23664464d3ace407d868ea7fb81d90ce7b5d620e 100644
--- a/src/OpenFOAM/primitives/Vector2D/Vector2D.H
+++ b/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;
};
diff --git a/src/OpenFOAM/primitives/Vector2D/Vector2DI.H b/src/OpenFOAM/primitives/Vector2D/Vector2DI.H
index b7387191cec2c29c3f57696d76535a2708bfbaab..46599dfc4875678a651476c550bc796fb48e20a4 100644
--- a/src/OpenFOAM/primitives/Vector2D/Vector2DI.H
+++ b/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);
}
diff --git a/src/functionObjects/field/mapFields/mapFields.H b/src/functionObjects/field/mapFields/mapFields.H
index a699b313929deb1b952209faa67d6b963ea7d471..873ab0899b15ec9d595f495fadd0d50e0b9e77fa 100644
--- a/src/functionObjects/field/mapFields/mapFields.H
+++ b/src/functionObjects/field/mapFields/mapFields.H
@@ -98,7 +98,6 @@ Usage
\verbatim
nearestFaceAMI
faceAreaWeightAMI
- partialFaceAreaWeightAMI
\endverbatim
The inherited entries are elaborated in:
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index 96ba959fa05c99d168da3d855d66c24e2865c0b8..7f6d4ea5406d924cfad6fa1397a10d1fe24a6ea7 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/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)
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index 7a2259066e8c4d71fecf43dd240c3bdf7999fe65..f691d6ace3c9bb430352ee0ac0f5237702a36e96 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/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;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
index 22823a23ec6d8d57f334b1731d8dba53d92506f9..1effda038534ba3eda4b9915dfbfcea305874db6 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolationI.H
+++ b/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_;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
index dc41c132ea60a990f89c16e796022146944c596d..7c03cbed1cf488a860169f1ee7166837c548cf08 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
+++ b/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;
-}
-
-
// ************************************************************************* //
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.H
index 84031daa733f3e53df6d86d6bdaca30a458ee3db..07e08b68d8c683397c4d8712bf5a2ab2b7be0f43 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.H
+++ b/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;
};
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C
index e6fbbb1816a1a9c7e970bafa02589d0082d0e4bf..395b1830f49d26b5268b3573950ff3b23c2d78b0 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C
+++ b/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;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.H
index 818d26ebea8aec4310209e3d46261a90384a3f4d..5c3d61dbd84305daabcba4ab49a9f74c03e1f859 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.H
+++ b/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
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
new file mode 100644
index 0000000000000000000000000000000000000000..c8270c6d313992f45c7f9b35dadd22ce17b7aeb9
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
@@ -0,0 +1,526 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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_)
+ {
+ return false;
+ }
+
+ addProfiling(ami, "faceAreaWeightAMI2D::calculate");
+
+ advancingFrontAMI::calculate(srcPatch, tgtPatch, surfPtr);
+
+ const auto& src = this->srcPatch();
+ const auto& tgt = this->tgtPatch(); // might be the extended patch!
+
+ bool validSize = true;
+
+ // Check that patch sizes are valid
+ if (!src.size())
+ {
+ validSize = false;
+ }
+ else if (!tgt.size())
+ {
+ WarningInFunction
+ << src.size() << " source faces but no target faces" << endl;
+
+ validSize = false;
+ }
+
+ srcCentroids_.setSize(srcAddress_.size());
+
+ // Temporary storage for addressing and weights
+ List<DynamicList<label>> tgtAddr(tgt.size());
+ List<DynamicList<scalar>> tgtWght(tgt.size());
+
+ // Find an approximate face length scale
+ const scalar lf(Cbb_*Foam::sqrt(gAverage(srcMagSf_)));
+
+ // Expansion to apply to source face bounding box
+ const vector d(lf*vector::one);
+
+ if (validSize)
+ {
+ // Create the tgt tree
+ const bool equalBinSize = true;
+ const label maxLevel = 10;
+ const label minBinSize = 4;
+ AABBTree<face> tree
+ (
+ tgt,
+ tgt.points(),
+ equalBinSize,
+ maxLevel,
+ minBinSize
+ );
+
+ const auto& tgtPoints = tgt.points();
+ List<boundBox> tgtFaceBbs(tgt.size());
+ forAll(tgt, facei)
+ {
+ tgtFaceBbs[facei] = boundBox(tgtPoints, tgt[facei], false);
+ }
+
+ DynamicList<label> nonOverlapFaces;
+
+ const auto& srcPoints = src.points();
+
+ forAll(src, srcFacei)
+ {
+ const face& srcFace = src[srcFacei];
+
+ treeBoundBox srcFaceBb(srcPoints, srcFace);
+ srcFaceBb.min() -= d;
+ srcFaceBb.max() += d;
+
+ const labelList tgtFaces
+ (
+ overlappingTgtFaces(tree, tgtFaceBbs, srcFaceBb)
+ );
+
+ DynamicList<label> srcAddr(tgtFaces.size());
+ DynamicList<scalar> srcWght(tgtFaces.size());
+ DynamicList<point> srcCtr(tgtFaces.size());
+
+ for (const label tgtFacei : tgtFaces)
+ {
+ storeInterArea
+ (
+ srcFacei,
+ tgtFacei,
+ srcAddr,
+ srcWght,
+ srcCtr,
+ tgtAddr[tgtFacei],
+ tgtWght[tgtFacei]
+ );
+ }
+
+ if (mustMatchFaces() && srcAddr.empty())
+ {
+ if (debug) writeNoMatch(srcFacei, tgtFaces, srcFaceBb);
+
+ // FatalErrorInFunction
+ // << "Unable to find match for source face " << srcFace
+ // << exit(FatalError);
+ }
+
+ srcAddress_[srcFacei].transfer(srcAddr);
+ srcWeights_[srcFacei].transfer(srcWght);
+ srcCentroids_[srcFacei].transfer(srcCtr);
+ }
+
+ srcNonOverlap_.transfer(nonOverlapFaces);
+
+ if (debug && !srcNonOverlap_.empty())
+ {
+ Pout<< " AMI: " << srcNonOverlap_.size()
+ << " non-overlap faces identified"
+ << endl;
+ }
+ }
+
+ // Transfer data to persistent storage
+
+ forAll(tgtAddr, i)
+ {
+ tgtAddress_[i].transfer(tgtAddr[i]);
+ tgtWeights_[i].transfer(tgtWght[i]);
+ }
+
+ if (distributed())
+ {
+ const primitivePatch& srcPatch0 = this->srcPatch0();
+ const primitivePatch& tgtPatch0 = this->tgtPatch0();
+
+ // Create global indexing for each original patch
+ globalIndex globalSrcFaces(srcPatch0.size());
+ globalIndex globalTgtFaces(tgtPatch0.size());
+
+ for (labelList& addressing : srcAddress_)
+ {
+ for (label& addr : addressing)
+ {
+ addr = extendedTgtFaceIDs_[addr];
+ }
+ }
+
+ for (labelList& addressing : tgtAddress_)
+ {
+ globalSrcFaces.inplaceToGlobal(addressing);
+ }
+
+ // Send data back to originating procs. Note that contributions
+ // from different processors get added (ListOps::appendEqOp)
+
+ mapDistributeBase::distribute
+ (
+ Pstream::commsTypes::nonBlocking,
+ List<labelPair>(),
+ tgtPatch0.size(),
+ extendedTgtMapPtr_->constructMap(),
+ false, // has flip
+ extendedTgtMapPtr_->subMap(),
+ false, // has flip
+ tgtAddress_,
+ labelList(),
+ ListOps::appendEqOp<label>(),
+ flipOp()
+ );
+
+ mapDistributeBase::distribute
+ (
+ Pstream::commsTypes::nonBlocking,
+ List<labelPair>(),
+ tgtPatch0.size(),
+ extendedTgtMapPtr_->constructMap(),
+ false,
+ extendedTgtMapPtr_->subMap(),
+ false,
+ tgtWeights_,
+ scalarList(),
+ ListOps::appendEqOp<scalar>(),
+ flipOp()
+ );
+
+ // Note: using patch face areas calculated by the AMI method
+ extendedTgtMapPtr_->reverseDistribute(tgtPatch0.size(), tgtMagSf_);
+
+ // Cache maps and reset addresses
+ List<Map<label>> cMapSrc;
+ srcMapPtr_.reset
+ (
+ new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
+ );
+
+ List<Map<label>> cMapTgt;
+ tgtMapPtr_.reset
+ (
+ new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
+ );
+ }
+
+ // Convert the weights from areas to normalised values
+ normaliseWeights(requireMatch_, true);
+
+ nonConformalCorrection();
+
+ upToDate_ = true;
+
+ return upToDate_;
+}
+
+
+void Foam::faceAreaWeightAMI2D::write(Ostream& os) const
+{
+ advancingFrontAMI::write(os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.H
similarity index 59%
rename from src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H
rename to src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.H
index f03f8d155dd4281eb33555170e4b0b1328f48206..3fc20bc7a9935477b466a95cbc3e2879ed0d298d 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.H
@@ -5,8 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2013-2016 OpenFOAM Foundation
- Copyright (C) 2016-2020 OpenCFD Ltd.
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -25,20 +24,22 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
- Foam::partialFaceAreaWeightAMI
+ Foam::faceAreaWeightAMI2D
Description
- Partial face area weighted Arbitrary Mesh Interface (AMI) method
+ Face area weighted Arbitrary Mesh Interface (AMI) method that performs
+ the intersection of src and tgt face area in 2D.
SourceFiles
- partialFaceAreaWeightAMI.C
+ faceAreaWeightAMI2D.C
\*---------------------------------------------------------------------------*/
-#ifndef partialFaceAreaWeightAMI_H
-#define partialFaceAreaWeightAMI_H
+#ifndef faceAreaWeightAMI2D_H
+#define faceAreaWeightAMI2D_H
-#include "faceAreaWeightAMI.H"
+#include "advancingFrontAMI.H"
+#include "AABBTree.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -46,58 +47,80 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class partialFaceAreaWeightAMI Declaration
+ Class faceAreaWeightAMI2D Declaration
\*---------------------------------------------------------------------------*/
-class partialFaceAreaWeightAMI
+class faceAreaWeightAMI2D
:
- public faceAreaWeightAMI
+ public advancingFrontAMI
{
+protected:
-private:
-
- // Private Member Functions
+ // Protected Data
- //- No copy assignment
- void operator=(const partialFaceAreaWeightAMI&) = delete;
+ //- Face bounding box factor
+ scalar Cbb_;
-protected:
-
// Protected Member Functions
- //- Set the source and target seed faces
- virtual bool setNextFaces
+ //- No copy assignment
+ void operator=(const faceAreaWeightAMI2D&) = delete;
+
+ //- Helper function to write non-matched source faces to the set
+ //- of candidate faces
+ void writeNoMatch
(
- label& startSeedi,
- label& srcFacei,
- label& tgtFacei,
- const bitSet& mapFlag,
- labelList& seedFaces,
- const DynamicList<label>& visitedFaces,
- const bool errorOnNotFound = true
+ const label srcFacei,
+ const labelList& tgtFaceCandidates,
+ const boundBox& srcFaceBb
) const;
+ // Evaluation
+
+ //- Calculate and store the area of intersection between source and
+ //- target faces
+ void storeInterArea
+ (
+ const label srcFacei,
+ const label tgtFacei,
+ DynamicList<label>& srcAddr,
+ DynamicList<scalar>& srcWght,
+ DynamicList<vector>& srcCtr,
+ DynamicList<label>& tgtAddr,
+ DynamicList<scalar>& tgtWght
+ ) const;
+
+
+ //- Return the set of tgt face IDs that overlap the src face bb
+ labelList overlappingTgtFaces
+ (
+ const AABBTree<face>& tree,
+ const List<boundBox>& tgtFaceBbs,
+ const boundBox& srcFaceBb
+ ) const;
+
+
public:
//- Runtime type information
- TypeName("partialFaceAreaWeightAMI");
+ TypeName("faceAreaWeightAMI2D");
// Constructors
//- Construct from dictionary
- partialFaceAreaWeightAMI
+ faceAreaWeightAMI2D
(
const dictionary& dict,
const bool reverseTarget = false
);
//- Construct from components
- partialFaceAreaWeightAMI
+ faceAreaWeightAMI2D
(
- const bool requireMatch = false,
+ const bool requireMatch,
const bool reverseTarget = false,
const scalar lowWeightCorrection = -1,
const faceAreaIntersect::triangulationMode triMode =
@@ -106,25 +129,21 @@ public:
);
//- Construct as copy
- partialFaceAreaWeightAMI(const partialFaceAreaWeightAMI& ami);
+ faceAreaWeightAMI2D(const faceAreaWeightAMI2D& ami);
//- Construct and return a clone
virtual autoPtr<AMIInterpolation> clone() const
{
- return
- autoPtr<AMIInterpolation>(new partialFaceAreaWeightAMI(*this));
+ return autoPtr<AMIInterpolation>(new faceAreaWeightAMI2D(*this));
}
//- Destructor
- virtual ~partialFaceAreaWeightAMI() = default;
+ virtual ~faceAreaWeightAMI2D() = default;
// Member Functions
- //- Flag to indicate that interpolation patches are conformal
- virtual bool conformal() const;
-
//- Update addressing, weights and (optional) centroids
virtual bool calculate
(
@@ -132,6 +151,9 @@ public:
const primitivePatch& tgtPatch,
const autoPtr<searchableSurface>& surfPtr = nullptr
);
+
+ //- Write
+ virtual void write(Ostream& os) const;
};
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
deleted file mode 100644
index 61b4197f2b0b5354f02bfa7c3dada538d4a46555..0000000000000000000000000000000000000000
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
+++ /dev/null
@@ -1,159 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | www.openfoam.com
- \\/ M anipulation |
--------------------------------------------------------------------------------
- Copyright (C) 2013-2016 OpenFOAM Foundation
- 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 "partialFaceAreaWeightAMI.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
- defineTypeNameAndDebug(partialFaceAreaWeightAMI, 0);
- addToRunTimeSelectionTable
- (
- AMIInterpolation,
- partialFaceAreaWeightAMI,
- dict
- );
- addToRunTimeSelectionTable
- (
- AMIInterpolation,
- partialFaceAreaWeightAMI,
- component
- );
-}
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-bool Foam::partialFaceAreaWeightAMI::setNextFaces
-(
- label& startSeedi,
- label& srcFacei,
- label& tgtFacei,
- const bitSet& mapFlag,
- labelList& seedFaces,
- const DynamicList<label>& visitedFaces,
- const bool errorOnNotFound
-) const
-{
- return faceAreaWeightAMI::setNextFaces
- (
- startSeedi,
- srcFacei,
- tgtFacei,
- mapFlag,
- seedFaces,
- visitedFaces,
- false // no error on not found
- );
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-Foam::partialFaceAreaWeightAMI::partialFaceAreaWeightAMI
-(
- const dictionary& dict,
- const bool reverseTarget
-)
-:
- faceAreaWeightAMI(dict, reverseTarget)
-{}
-
-
-Foam::partialFaceAreaWeightAMI::partialFaceAreaWeightAMI
-(
- const bool requireMatch,
- const bool reverseTarget,
- const scalar lowWeightCorrection,
- const faceAreaIntersect::triangulationMode triMode,
- const bool restartUncoveredSourceFace
-)
-:
- faceAreaWeightAMI
- (
- requireMatch,
- reverseTarget,
- lowWeightCorrection,
- triMode,
- restartUncoveredSourceFace
- )
-{}
-
-
-Foam::partialFaceAreaWeightAMI::partialFaceAreaWeightAMI
-(
- const partialFaceAreaWeightAMI& ami
-)
-:
- faceAreaWeightAMI(ami)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
-bool Foam::partialFaceAreaWeightAMI::conformal() const
-{
- return false;
-}
-
-
-bool Foam::partialFaceAreaWeightAMI::calculate
-(
- const primitivePatch& srcPatch,
- const primitivePatch& tgtPatch,
- const autoPtr<searchableSurface>& surfPtr
-)
-{
- if (faceAreaWeightAMI::calculate(srcPatch, tgtPatch, surfPtr))
- {
- if (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);
- }
- }
-
- return true;
- }
-
- return false;
-}
-
-
-// ************************************************************************* //
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
index 7813526c858bccc1e212fcf79a1fb0896ab0a700..a12b5ddeaa3e91bb77afc3ec822bdc33bba50739 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
+++ b/src/meshTools/AMIInterpolation/patches/cyclicACMI/cyclicACMIPolyPatch/cyclicACMIPolyPatch.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2013-2016 OpenFOAM Foundation
- Copyright (C) 2018-2020 OpenCFD Ltd.
+ Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -63,7 +63,6 @@ SourceFiles
#include "cyclicAMIPolyPatch.H"
#include "AMIPatchToPatchInterpolation.H"
#include "polyBoundaryMesh.H"
-#include "partialFaceAreaWeightAMI.H"
#include "PatchFunction1.H"
#include "uniformDimensionedFields.H"
#include "vectorList.H"
@@ -197,7 +196,7 @@ public:
const polyBoundaryMesh& bm,
const word& patchType,
const transformType transform = UNKNOWN,
- const word& defaultAMIMethod = partialFaceAreaWeightAMI::typeName
+ const word& defaultAMIMethod = faceAreaWeightAMI::typeName
);
//- Construct from dictionary
@@ -208,7 +207,7 @@ public:
const label index,
const polyBoundaryMesh& bm,
const word& patchType,
- const word& defaultAMIMethod = partialFaceAreaWeightAMI::typeName
+ const word& defaultAMIMethod = faceAreaWeightAMI::typeName
);
//- Construct as copy, resetting the boundary mesh
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
index bc63d3ed1e51d023f45a226164ce7f6e8f66a1fa..a7e652c6b8c93dfa4ebbab2831ac4755c9a08add 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2015-2020 OpenCFD Ltd.
+ Copyright (C) 2015-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
\*---------------------------------------------------------------------------*/
#include "cyclicPeriodicAMIPolyPatch.H"
-#include "partialFaceAreaWeightAMI.H"
#include "addToRunTimeSelectionTable.H"
// For debugging
@@ -560,14 +559,16 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
bm,
patchType,
transform,
- partialFaceAreaWeightAMI::typeName
+ faceAreaWeightAMI::typeName
),
periodicPatchName_(word::null),
periodicPatchID_(-1),
nTransforms_(0),
nSectors_(0),
maxIter_(36)
-{}
+{
+ AMIPtr_->setRequireMatch(false);
+}
Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -586,14 +587,16 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
index,
bm,
patchType,
- partialFaceAreaWeightAMI::typeName
+ faceAreaWeightAMI::typeName
),
periodicPatchName_(dict.lookup("periodicPatch")),
periodicPatchID_(-1),
nTransforms_(dict.getOrDefault<label>("nTransforms", 0)),
nSectors_(dict.getOrDefault<label>("nSectors", 0)),
maxIter_(dict.getOrDefault<label>("maxIter", 36))
-{}
+{
+ AMIPtr_->setRequireMatch(false);
+}
Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -608,7 +611,9 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
nTransforms_(pp.nTransforms_),
nSectors_(pp.nSectors_),
maxIter_(pp.maxIter_)
-{}
+{
+ AMIPtr_->setRequireMatch(false);
+}
Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -627,7 +632,9 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
nTransforms_(pp.nTransforms_),
nSectors_(pp.nSectors_),
maxIter_(pp.maxIter_)
-{}
+{
+ AMIPtr_->setRequireMatch(false);
+}
Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
@@ -645,7 +652,9 @@ Foam::cyclicPeriodicAMIPolyPatch::cyclicPeriodicAMIPolyPatch
nTransforms_(pp.nTransforms_),
nSectors_(pp.nSectors_),
maxIter_(pp.maxIter_)
-{}
+{
+ AMIPtr_->setRequireMatch(false);
+}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
diff --git a/src/meshTools/AMIInterpolation/triangle2D/triangle2D.C b/src/meshTools/AMIInterpolation/triangle2D/triangle2D.C
new file mode 100644
index 0000000000000000000000000000000000000000..d0c92c24eda50e4dd2b3fa6f09825a7554f23fb8
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/triangle2D/triangle2D.C
@@ -0,0 +1,303 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020-2021 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 "triangle2D.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::triangle2D::debug = 0;
+
+Foam::scalar Foam::triangle2D::relTol = 1e-8;
+
+Foam::scalar Foam::triangle2D::absTol = 1e-10;
+
+Foam::FixedList<Foam::vector2D, 7> Foam::triangle2D::work_
+(
+ vector2D::zero
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::triangle2D::triangle2D
+(
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c,
+ const bool orient
+)
+:
+ FixedList<vector2D, 3>({a, b, c}),
+ area_(area(a, b, c))
+{
+ if (orient && area_ < 0)
+ {
+ // Inverted triangle
+ triangle2D& tri = *this;
+ vector2D tmp = tri[0];
+ tri[0] = tri[2];
+ tri[2] = tmp;
+
+ area_ = mag(area_);
+ }
+}
+
+
+Foam::triangle2D::triangle2D
+(
+ const vector& a3d,
+ const vector& b3d,
+ const vector& c3d,
+ const vector& axis1,
+ const vector& axis2,
+ const bool orient
+)
+:
+ triangle2D
+ (
+ vector2D(axis1 & a3d, axis2 & a3d),
+ vector2D(axis1 & b3d, axis2 & b3d),
+ vector2D(axis1 & c3d, axis2 & c3d),
+ orient
+ )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::label Foam::triangle2D::snapClosePoints(const triangle2D& triB)
+{
+ label nSnapped = 0;
+ triangle2D& triA = *this;
+
+ FixedList<bool, 3> match(true);
+
+ for (auto& a : triA)
+ {
+ forAll(triB, tb)
+ {
+ if (match[tb] && a.isClose(triB[tb]))
+ {
+ a = triB[tb];
+ match[tb] = false;
+ ++nSnapped;
+ break;
+ }
+ }
+ }
+
+ return nSnapped;
+}
+
+
+void Foam::triangle2D::interArea
+(
+ const triangle2D& triB,
+ vector2D& centre,
+ scalar& area
+) const
+{
+ const triangle2D& triA = *this;
+
+ // Potential short-cut if the triangles are the same-ish. Happens rarely
+ // for moving mesh cases.
+ // if (nClosePoints(triA, triB) == 3)
+ // {
+ // centre = triA.centre();
+ // area = triA.area();
+ // return;
+ // }
+
+ if (debug)
+ {
+ static int nInter = 0;
+ OFstream os("intersection-tris-"+Foam::name(nInter++)+".obj");
+ writeOBJ(os, triA, 0);
+ writeOBJ(os, triB, 3);
+ Info<< "written " << os.name() << endl;
+ }
+
+
+ // Use work_ to store intersections
+
+ // Current number of intersections
+ int nPoint = 0;
+
+ static FixedList<vector2D, 7> work2;
+
+ // Clipped triangle starts as triA
+ work2[0] = triA[0];
+ work2[1] = triA[1];
+ work2[2] = triA[2];
+
+ int nPoint2 = 3;
+
+
+ vector2D intersection(Zero);
+
+ // Cut triA using triB's edges as clipping planes
+ for (int i0 = 0; i0 <= 2; ++i0)
+ {
+ if (debug)
+ {
+ Info<< "\nstarting points:";
+ for (int i = 0; i < nPoint2; ++i)
+ {
+ Info<< work2[i];
+ }
+ Info<< endl;
+ }
+
+ // Clipping plane points
+ const label i1 = (i0 + 1) % 3;
+ const vector2D& c0 = triB[i0];
+ const vector2D& c1 = triB[i1];
+
+ nPoint = 0;
+
+ // Test against all intersection poly points
+ for (int j = 0; j < nPoint2; ++j)
+ {
+ const vector2D& p0 = work2[j];
+ const vector2D& p1 = work2[(j+1) % nPoint2];
+
+ if (triangle2D(c0, c1, p0).order() == 1)
+ {
+ if (triangle2D(c0, c1, p1).order() == 1)
+ {
+ work_[nPoint++] = p1;
+ }
+ else
+ {
+ if (lineIntersectionPoint(p0, p1, c0, c1, intersection))
+ {
+ work_[nPoint++] = intersection;
+ }
+ }
+ }
+ else
+ {
+ if (triangle2D(c0, c1, p1).order() == 1)
+ {
+ if (lineIntersectionPoint(p0, p1, c0, c1, intersection))
+ {
+ work_[nPoint++] = intersection;
+ }
+ work_[nPoint++] = p1;
+ }
+ }
+ }
+
+ work2 = work_;
+ nPoint2 = nPoint;
+ }
+
+ if (debug)
+ {
+ static int n = 0;
+ OFstream os("intersection-poly-"+Foam::name(n++)+".obj");
+ for (int i = 0; i < nPoint; ++i)
+ {
+ os << "v " << work_[i].x() << " " << work_[i].y() << " 0" << nl;
+ }
+ os << "l";
+ for (int i = 0; i < nPoint; ++i)
+ {
+ os << " " << (i + 1);
+ }
+ os << " 1" << endl;
+ Info<< "written " << os.name() << endl;
+
+ Info<< "Intersection points:" << endl;
+ for (int i = 0; i < nPoint; ++i)
+ {
+ Info<< " " << work_[i] << endl;
+ }
+ }
+
+ // Calculate the area of the clipped triangle
+ scalar twoArea = 0;
+ centre = vector2D::zero;
+ if (nPoint)
+ {
+ for (int i = 0; i < nPoint - 1; ++i)
+ {
+ twoArea += work_[i].x()*work_[i+1].y();
+ twoArea -= work_[i].y()*work_[i+1].x();
+ centre += work_[i];
+ }
+ twoArea += work_[nPoint-1].x()*work_[0].y();
+ twoArea -= work_[nPoint-1].y()*work_[0].x();
+
+ centre += work_[nPoint - 1];
+ centre /= scalar(nPoint);
+ }
+
+ area = 0.5*twoArea;
+}
+
+
+Foam::scalar Foam::triangle2D::interArea(const triangle2D& triB) const
+{
+ vector2D dummyCentre(Zero);
+ scalar area;
+
+ interArea(triB, dummyCentre, area);
+
+ return area;
+}
+
+
+bool Foam::triangle2D::overlaps(const triangle2D& triB) const
+{
+ const triangle2D& triA = *this;
+
+ // True if any of triB's edges intersect a triA edge
+ for (int i = 0; i < 3; ++i)
+ {
+ int i1 = (i + 1) % 3;
+
+ for (int j = 0; j < 3; ++j)
+ {
+ int j1 = (j + 1) % 3;
+
+ if (lineIntersects(triA[i], triA[i1], triB[j], triB[j1]))
+ {
+ return true;
+ }
+ }
+ }
+
+ return
+ (nClosePoints(triA, triB) == 3) // same tri
+ || triA.contains(triB) // triA contains triB
+ || triB.contains(triA); // triB contains triA
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/AMIInterpolation/triangle2D/triangle2D.H b/src/meshTools/AMIInterpolation/triangle2D/triangle2D.H
new file mode 100644
index 0000000000000000000000000000000000000000..c85385d6165aae0f1766ea634aab9299c6f51d39
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/triangle2D/triangle2D.H
@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020-2021 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/>.
+
+Class
+ Foam::triangle2D
+
+Description
+ 2-D triangle and queries
+
+SourceFiles
+ triangle2DI.H
+ triangle2D.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triangle2D_H
+#define triangle2D_H
+
+#include "vector.H"
+#include "vector2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class triangle2D Declaration
+\*---------------------------------------------------------------------------*/
+
+class triangle2D
+:
+ public FixedList<vector2D, 3>
+{
+ // Private Data
+
+ //- Area
+ scalar area_;
+
+ //- Helper for intersections
+ static FixedList<vector2D, 7> work_;
+
+
+public:
+
+ static int debug;
+
+ //- Relative tolerance
+ static scalar relTol;
+
+ //- Absolute tolerance
+ static scalar absTol;
+
+ //- Construct from 3 2-D points
+ triangle2D
+ (
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c,
+ const bool orient = false
+ );
+
+ //- Construct from 3 3-D points and axes
+ triangle2D
+ (
+ const vector& a3d,
+ const vector& b3d,
+ const vector& c3d,
+ const vector& axis1,
+ const vector& axis2,
+ const bool orient = false
+ );
+
+
+ // Member Functions
+
+ //- Returns:
+ //- 1 if points are ordered in anti-clockwise direction
+ //- -1 if points are ordered in clockwise direction
+ //- 0 if the triangle has collapsed to a line
+ inline label order() const;
+
+ //- Write the triangle in OBJ format
+ inline static void writeOBJ
+ (
+ Ostream& os,
+ const triangle2D& tri,
+ label offset
+ );
+
+ //- Return the number of similar points between two triangles
+ inline static label nClosePoints
+ (
+ const triangle2D& triA,
+ const triangle2D& triB
+ );
+
+ //- Return the signed area
+ inline static scalar area
+ (
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c
+ );
+
+ //- Set the intersection between a line and segment
+ //- Return true if lines intersect
+ inline static bool lineSegmentIntersectionPoint
+ (
+ const vector2D& lp1,
+ const vector2D& lp2,
+ const vector2D& sp1,
+ const vector2D& sp2,
+ vector2D& intersection
+ );
+
+ //- Set the intersection between two lines
+ //- Return true if lines intersect
+ inline static bool lineIntersectionPoint
+ (
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c,
+ const vector2D& d,
+ vector2D& intersection
+ );
+
+ //- Return true if lines ab and cd intersect
+ inline static bool lineIntersects
+ (
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c,
+ const vector2D& d
+ );
+
+ //- Snap [this] triangle's points to those of triB if they are within
+ //- absTol
+ // Returns the number of snapped points
+ label snapClosePoints(const triangle2D& triB);
+
+ //- Return the intersection centre and area
+ void interArea
+ (
+ const triangle2D& triB,
+ vector2D& centre,
+ scalar& area
+ ) const;
+
+ //- Return the intersection area
+ scalar interArea(const triangle2D& triB) const;
+
+ //- Return true if triB overlaps
+ bool overlaps(const triangle2D& triB) const;
+
+ //- Return the triangle area
+ inline scalar area() const noexcept;
+
+ //- Return the triangle centre
+ inline vector2D centre() const;
+
+ //- Return true if tri is within this triangle
+ inline bool contains(const triangle2D& tri) const;
+
+ //- Return true if triB is the same as this triangle
+ inline bool isSame(const triangle2D& triB) const;
+
+ //- Return true if t point p is inside this triangle
+ inline bool pointInside(const vector2D& p) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "triangle2DI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H b/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H
new file mode 100644
index 0000000000000000000000000000000000000000..6f73594f60f3996160df5501c80cbb04a91825eb
--- /dev/null
+++ b/src/meshTools/AMIInterpolation/triangle2D/triangle2DI.H
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020-2021 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 "OFstream.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+inline void Foam::triangle2D::writeOBJ
+(
+ Ostream& os,
+ const triangle2D& tri,
+ label offset
+)
+{
+ os << "v " << tri[0].x() << " " << tri[0].y() << " 0" << nl
+ << "v " << tri[1].x() << " " << tri[1].y() << " 0" << nl
+ << "v " << tri[2].x() << " " << tri[2].y() << " 0" << nl
+ << "l"
+ << " " << 1 + offset
+ << " " << 2 + offset
+ << " " << 3 + offset
+ << " " << 1 + offset
+ << endl;
+}
+
+
+inline Foam::label Foam::triangle2D::nClosePoints
+(
+ const triangle2D& triA,
+ const triangle2D& triB
+)
+{
+ label nClose = 0;
+
+ FixedList<bool, 3> match(true);
+
+ for (auto& a : triA)
+ {
+ forAll(triB, tb)
+ {
+ if (match[tb] && a.isClose(triB[tb]))
+ {
+ match[tb] = false;
+ ++nClose;
+ break;
+ }
+ }
+ }
+
+ return nClose;
+}
+
+
+inline Foam::scalar Foam::triangle2D::area
+(
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c
+)
+{
+ const vector2D e1(a - c);
+ const vector2D e2(b - c);
+
+ return 0.5*e1.perp(e2);
+}
+
+
+inline Foam::vector2D Foam::triangle2D::centre() const
+{
+ const triangle2D& tri = *this;
+
+ return (1/3.0)*(tri[0] + tri[1] + tri[2]);
+}
+
+
+inline bool Foam::triangle2D::lineSegmentIntersectionPoint
+(
+ const vector2D& lp1,
+ const vector2D& lp2,
+ const vector2D& sp1,
+ const vector2D& sp2,
+ vector2D& intersection
+)
+{
+ const vector2D r(lp2 - lp1);
+ const vector2D s(sp2 - sp1);
+
+ const scalar rcs = r.perp(s);
+
+ if (mag(rcs) > ROOTVSMALL)
+ {
+ const scalar u = (sp1 - lp1).perp(r)/rcs;
+
+ if (u >= -relTol && u <= 1+relTol)
+ {
+ intersection = sp1 + u*s;
+ return true;
+ }
+ }
+
+/*
+ // Collapsed to line
+ if (mag(triangle2D(lp1, lp2, sp1).area()) < absTol)
+ {
+ intersection = sp1;
+ return true;
+ }
+
+ if (mag(triangle2D(lp1, lp2, sp2).area()) < absTol)
+ {
+ intersection = sp2;
+ return true;
+ }
+*/
+
+ if (debug)
+ {
+ OFstream os("bad-intersection.obj");
+ os << "v " << lp1.x() << " " << lp1.y() << " 0" << nl
+ << "v " << lp2.x() << " " << lp2.y() << " 0" << nl
+ << "v " << sp1.x() << " " << sp1.y() << " 0" << nl
+ << "v " << sp2.x() << " " << sp2.y() << " 0" << nl
+ << "l 1 2"
+ << "l 3 4"
+ << endl;
+ }
+
+ return false;
+}
+
+
+inline bool Foam::triangle2D::lineIntersectionPoint
+(
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c,
+ const vector2D& d,
+ vector2D& intersection
+)
+{
+ return lineSegmentIntersectionPoint(c, d, a, b, intersection);
+}
+
+
+inline bool Foam::triangle2D::lineIntersects
+(
+ const vector2D& a,
+ const vector2D& b,
+ const vector2D& c,
+ const vector2D& d
+)
+{
+ if
+ (
+ (triangle2D(a, c, d).order() == triangle2D(b, c, d).order())
+ || (triangle2D(a, b, c).order() == triangle2D(a, b, d).order())
+ )
+ {
+ return false;
+ }
+
+
+ DebugInfo<< "line " << a << b << " intersects " << c << d << endl;
+
+ return true;
+}
+
+
+inline Foam::scalar Foam::triangle2D::area() const noexcept
+{
+ return area_;
+}
+
+
+inline Foam::label Foam::triangle2D::order() const
+{
+ return mag(area_) < SMALL ? 0 : sign(area_);
+}
+
+
+inline bool Foam::triangle2D::contains(const triangle2D& tri) const
+{
+ return
+ pointInside(tri[0])
+ && pointInside(tri[1])
+ && pointInside(tri[2]);
+}
+
+
+inline bool Foam::triangle2D::isSame(const triangle2D& triB) const
+{
+ const triangle2D& triA = *this;
+
+ return
+ triA[0].isClose(triB[0])
+ && triA[1].isClose(triB[1])
+ && triA[2].isClose(triB[2]);
+}
+
+
+inline bool Foam::triangle2D::pointInside(const vector2D& p) const
+{
+ const triangle2D& triA = *this;
+
+ for (int i = 0; i < 3; ++i)
+ {
+ if ((triA[(i + 1) % 3] - triA[i]).perp(p - triA[i]) < 0)
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 3d32aa57367f12a8cf17512276d4115acf9a3f5c..d92cd130ad11bd503cb2289ca69461f7b9353ea4 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -269,7 +269,6 @@ $(AMI)/AMIInterpolation/AMIInterpolationNew.C
$(AMI)/AMIInterpolation/advancingFrontAMI/advancingFrontAMI.C
$(AMI)/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
$(AMI)/AMIInterpolation/faceAreaWeightAMI/faceAreaWeightAMI.C
-$(AMI)/AMIInterpolation/partialFaceAreaWeightAMI/partialFaceAreaWeightAMI.C
$(AMI)/AMIInterpolation/nearestFaceAMI/nearestFaceAMI.C
$(AMI)/faceAreaIntersect/faceAreaIntersect.C
$(AMI)/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
@@ -277,6 +276,10 @@ $(AMI)/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceFi
$(AMI)/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
$(AMI)/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.C
+$(AMI)/triangle2D/triangle2D.C
+$(AMI)/AMIInterpolation/faceAreaWeightAMI2D/faceAreaWeightAMI2D.C
+
+
AMICycPatches=$(AMI)/patches/cyclicAMI
$(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterface.C
$(AMICycPatches)/cyclicAMILduInterfaceField/cyclicAMILduInterfaceField.C