From 94679fa88d9b764ed83736ee96630a9e1d00d102 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Tue, 12 Apr 2022 14:56:46 +0100 Subject: [PATCH] meshTools: Added patchToPatch class for patch mapping --- src/meshTools/Make/files | 6 + .../PrimitiveOldTimePatch.C | 235 ++++ .../PrimitiveOldTimePatch.H | 230 ++++ .../primitiveOldTimePatch.H | 53 + .../uindirectPrimitiveOldTimePatch.H | 53 + .../intersection/intersectionPatchToPatch.C | 937 +++++++++++++++ .../intersection/intersectionPatchToPatch.H | 424 +++++++ .../intersection/intersectionPatchToPatchI.H | 60 + .../inverseDistancePatchToPatch.C | 284 +++++ .../inverseDistancePatchToPatch.H | 178 +++ .../nearest/nearestPatchToPatch.C | 303 +++++ .../nearest/nearestPatchToPatch.H | 157 +++ .../patchToPatch/patchToPatch/patchToPatch.C | 1063 +++++++++++++++++ .../patchToPatch/patchToPatch/patchToPatch.H | 485 ++++++++ .../patchToPatch/patchToPatch/patchToPatchI.H | 125 ++ .../patchToPatch/patchToPatchParallelOps.C | 469 ++++++++ .../patchToPatch/rays/raysPatchToPatch.C | 332 +++++ .../patchToPatch/rays/raysPatchToPatch.H | 204 ++++ 18 files changed, 5598 insertions(+) create mode 100644 src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.C create mode 100644 src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.H create mode 100644 src/meshTools/PrimitiveOldTimePatch/primitiveOldTimePatch.H create mode 100644 src/meshTools/PrimitiveOldTimePatch/uindirectPrimitiveOldTimePatch.H create mode 100644 src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C create mode 100644 src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.H create mode 100644 src/meshTools/patchToPatch/intersection/intersectionPatchToPatchI.H create mode 100644 src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.C create mode 100644 src/meshTools/patchToPatch/inverseDistance/inverseDistancePatchToPatch.H create mode 100644 src/meshTools/patchToPatch/nearest/nearestPatchToPatch.C create mode 100644 src/meshTools/patchToPatch/nearest/nearestPatchToPatch.H create mode 100644 src/meshTools/patchToPatch/patchToPatch/patchToPatch.C create mode 100644 src/meshTools/patchToPatch/patchToPatch/patchToPatch.H create mode 100644 src/meshTools/patchToPatch/patchToPatch/patchToPatchI.H create mode 100644 src/meshTools/patchToPatch/patchToPatch/patchToPatchParallelOps.C create mode 100644 src/meshTools/patchToPatch/rays/raysPatchToPatch.C create mode 100644 src/meshTools/patchToPatch/rays/raysPatchToPatch.H diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index bf2b8ec1f9..317a015345 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -256,6 +256,12 @@ $(AMIOverlapPatches)/cyclicRepeatAMIPolyPatch/cyclicRepeatAMIPolyPatch.C $(AMIOverlapPatches)/cyclicRepeatAMIPointPatch/cyclicRepeatAMIPointPatch.C $(AMIOverlapPatches)/cyclicRepeatAMIPointPatchField/cyclicRepeatAMIPointPatchFields.C +patchToPatch/patchToPatch/patchToPatch.C +patchToPatch/patchToPatch/patchToPatchParallelOps.C +patchToPatch/nearest/nearestPatchToPatch.C +patchToPatch/inverseDistance/inverseDistancePatchToPatch.C +patchToPatch/intersection/intersectionPatchToPatch.C + mappedPatches/mappedPolyPatch/mappedPatchBase.C mappedPatches/mappedPolyPatch/mappedPolyPatch.C mappedPatches/mappedPolyPatch/mappedWallPolyPatch.C diff --git a/src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.C b/src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.C new file mode 100644 index 0000000000..3fa039db98 --- /dev/null +++ b/src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.C @@ -0,0 +1,235 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "PrimitiveOldTimePatch.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::PrimitiveOldTimePatch::PrimitiveOldTimePatch +( + const FaceList& faces, + const Field& points, + const Field& points0 +) +: + PrimitivePatch(faces, points), + points0_(points0), + patch0Ptr_(new patch0Type(faces, points0_)), + localPoints0Ptr_(nullptr) +{} + + +template +Foam::PrimitiveOldTimePatch::PrimitiveOldTimePatch +( + const PrimitivePatch& patch, + const Field& points0 +) +: + PrimitivePatch(patch), + points0_(points0), + patch0Ptr_(new patch0Type(patch, points0_)), + localPoints0Ptr_(nullptr) +{} + + +template +Foam::PrimitiveOldTimePatch::PrimitiveOldTimePatch +( + const FaceList& faces, + const Field& points +) +: + PrimitivePatch(faces, points), + points0_(NullObjectRef>()), + patch0Ptr_(nullptr), + localPoints0Ptr_(nullptr) +{} + + +template +Foam::PrimitiveOldTimePatch::PrimitiveOldTimePatch +( + const PrimitivePatch& patch +) +: + PrimitivePatch(patch), + points0_(NullObjectRef>()), + patch0Ptr_(nullptr), + localPoints0Ptr_(nullptr) +{} + + +template +Foam::PrimitiveOldTimePatch::PrimitiveOldTimePatch +( + const PrimitiveOldTimePatch& patch +) +: + PrimitivePatch(patch), + points0_(patch.points0_), + patch0Ptr_(patch.patch0Ptr_, false), + localPoints0Ptr_(nullptr) +{} + + +template +Foam::PrimitiveOldTimePatch::PrimitiveOldTimePatch +( + PrimitiveOldTimePatch&& patch +) +: + PrimitivePatch(move(patch)), + points0_(move(patch.points0_)), + patch0Ptr_(patch.patch0Ptr_), + localPoints0Ptr_(nullptr) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::PrimitiveOldTimePatch::~PrimitiveOldTimePatch() +{ + clearOut(); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +template +const Foam::Field +< + typename Foam::PrimitiveOldTimePatch::PointType +>& Foam::PrimitiveOldTimePatch::localPoints0() const +{ + // !!! Cannot just call patch0Ptr_->localPoints() as this would generate + // topology in patch0Ptr_() that is already available in the base class. + // For now, we just duplicate the implementation in PrimitivePatch. + + if (!localPoints0Ptr_) + { + const labelList& meshPts = this->meshPoints(); + + localPoints0Ptr_ = new Field(meshPts.size()); + + Field& locPts = *localPoints0Ptr_; + + forAll(meshPts, pointi) + { + locPts[pointi] = points0_[meshPts[pointi]]; + } + } + + // But, it would be preferable to add a method to PrimitivePatch which + // calculates the local points given a list of points different to those + // that are stored. Then the implementations could be shared and we could + // do this: + /* + if (!localPoints0Ptr_) + { + localPoints0Ptr_ = this->calcLocalPoints(points0_); + } + */ + + return *localPoints0Ptr_; +} + + +template +const Foam::Field +< + typename Foam::PrimitiveOldTimePatch::PointType +>& Foam::PrimitiveOldTimePatch::faceCentres0() const +{ + return patch0Ptr_->faceCentres(); +} + + +template +const Foam::Field +< + typename Foam::PrimitiveOldTimePatch::PointType +>& Foam::PrimitiveOldTimePatch::faceAreas0() const +{ + return patch0Ptr_->faceAreas(); +} + + +template +const Foam::Field +< + typename Foam::PrimitiveOldTimePatch::PointType +>& Foam::PrimitiveOldTimePatch::faceNormals0() const +{ + return patch0Ptr_->faceNormals(); +} + + +template +const Foam::Field +< + typename Foam::PrimitiveOldTimePatch::PointType +>& Foam::PrimitiveOldTimePatch::pointNormals0() const +{ + // !!! See comments in localPoints0. This isn't needed for now. + + NotImplemented; + return NullObjectRef>(); +} + + +template +void Foam::PrimitiveOldTimePatch::clearOut() +{ + PrimitivePatch::clearOut(); + if (has0()) patch0Ptr_->clearOut(); + + deleteDemandDrivenData(localPoints0Ptr_); +} + + +template +void Foam::PrimitiveOldTimePatch::clearGeom() +{ + PrimitivePatch::clearGeom(); + if (has0()) patch0Ptr_->clearGeom(); + + deleteDemandDrivenData(localPoints0Ptr_); +} + + +template +void Foam::PrimitiveOldTimePatch::movePoints0 +( + const Field& +) +{ + if (has0()) patch0Ptr_->clearGeom(); +} + + +// ************************************************************************* // diff --git a/src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.H b/src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.H new file mode 100644 index 0000000000..e08271688d --- /dev/null +++ b/src/meshTools/PrimitiveOldTimePatch/PrimitiveOldTimePatch.H @@ -0,0 +1,230 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::PrimitiveOldTimePatch + +Description + +SourceFiles + PrimitiveOldTimePatch.C + +\*---------------------------------------------------------------------------*/ + +#ifndef PrimitiveOldTimePatch_H +#define PrimitiveOldTimePatch_H + +#include "PrimitivePatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Trait to determine what face list should be used for the old-time patch +template +struct UFaceList; +template +struct UFaceList> { typedef UList type; }; +template +struct UFaceList> : public UFaceList> {}; +template +struct UFaceList> : public UFaceList> {}; +template +struct UFaceList> { typedef UIndirectList type; }; +template +struct UFaceList> : public UFaceList> {}; + +/*---------------------------------------------------------------------------*\ + Class PrimitiveOldTimePatch Declaration +\*---------------------------------------------------------------------------*/ + +template +class PrimitiveOldTimePatch +: + public PrimitivePatch +{ +public: + + // Public Typedefs + + using typename PrimitivePatch::FaceListType; + + using typename PrimitivePatch::FaceType; + + using typename PrimitivePatch::PointFieldType; + + using typename PrimitivePatch::PointType; + + +private: + + // Private Typedefs + + typedef + PrimitivePatch + < + typename UFaceList::type, + const PointFieldType& + > + patch0Type; + + + // Private Data + + //- Reference to global list of old-time points + PointField points0_; + + //- Engine for calculating old-time geometry. Note: Methods that + // generate topology should not be called here. The base patch should + // have all the necessary topology available. The calculation of some + // geometric quantities requires topological data to be available. In + // these cases, special steps need to be taken here to make sure that + // the old-time patch object does not generate duplicate topology. + autoPtr patch0Ptr_; + + //- Points local to patch + mutable Field* localPoints0Ptr_; + + +public: + + // Constructors + + //- Construct from components + PrimitiveOldTimePatch + ( + const FaceList& faces, + const Field& points, + const Field& points0 + ); + + //- Construct from patch and old-time points + PrimitiveOldTimePatch + ( + const PrimitivePatch& patch, + const Field& points0 + ); + + //- Construct from components without old-time points + PrimitiveOldTimePatch + ( + const FaceList& faces, + const Field& points + ); + + //- Construct from patch without old-time points + PrimitiveOldTimePatch + ( + const PrimitivePatch& patch + ); + + //- Copy constructor + PrimitiveOldTimePatch(const PrimitiveOldTimePatch& patch); + + //- Move constructor + PrimitiveOldTimePatch(PrimitiveOldTimePatch&& patch); + + //- Construct and return a clone + autoPtr> clone() const + { + return autoPtr> + ( + new PrimitiveOldTimePatch(*this) + ); + } + + + //- Destructor + virtual ~PrimitiveOldTimePatch(); + + + // Member Functions + + // Access + + //- Return whether or not old-time geometry is available + bool has0() const + { + return patch0Ptr_.valid(); + } + + //- Return reference to old-time global points + const Field& points0() const + { + return points0_; + } + + + // Geometry + + //- Return pointField of old-time points in patch + const Field& localPoints0() const; + + //- Return old-time face centres for patch + const Field& faceCentres0() const; + + //- Return old-time face areas for patch + const Field& faceAreas0() const; + + //- Return old-time face normals for patch + const Field& faceNormals0() const; + + //- Return old-time point normals for patch + const Field& pointNormals0() const; + + + // Edit + + //- ... + void clearOut(); + + //- ... + void clearGeom(); + + //- Correct patch after moving points + virtual void movePoints0(const Field&); + + + // Member Operators + + //- Disallow default bitwise assignment + void operator=(const PrimitiveOldTimePatch&) = delete; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "PrimitiveOldTimePatch.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/PrimitiveOldTimePatch/primitiveOldTimePatch.H b/src/meshTools/PrimitiveOldTimePatch/primitiveOldTimePatch.H new file mode 100644 index 0000000000..27a36e3366 --- /dev/null +++ b/src/meshTools/PrimitiveOldTimePatch/primitiveOldTimePatch.H @@ -0,0 +1,53 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Typedef + Foam::primitiveOldTimePatch + +Description + Addressing for a faceList slice. + +\*---------------------------------------------------------------------------*/ + +#ifndef primitiveOldTimePatch_H +#define primitiveOldTimePatch_H + +#include "PrimitiveOldTimePatch.H" +#include "face.H" +#include "SubList.H" +#include "pointField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef + PrimitiveOldTimePatch, const pointField&> + primitiveOldTimePatch; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/PrimitiveOldTimePatch/uindirectPrimitiveOldTimePatch.H b/src/meshTools/PrimitiveOldTimePatch/uindirectPrimitiveOldTimePatch.H new file mode 100644 index 0000000000..d249d78cc2 --- /dev/null +++ b/src/meshTools/PrimitiveOldTimePatch/uindirectPrimitiveOldTimePatch.H @@ -0,0 +1,53 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Typedef + Foam::uindirectPrimitiveOldTimePatch + +Description + Addressing for a faceList slice. + +\*---------------------------------------------------------------------------*/ + +#ifndef uindirectPrimitiveOldTimePatch_H +#define uindirectPrimitiveOldTimePatch_H + +#include "PrimitiveOldTimePatch.H" +#include "face.H" +#include "UIndirectList.H" +#include "pointField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef + PrimitiveOldTimePatch, const pointField&> + uindirectPrimitiveOldTimePatch; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C b/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C new file mode 100644 index 0000000000..0c4200209f --- /dev/null +++ b/src/meshTools/patchToPatch/intersection/intersectionPatchToPatch.C @@ -0,0 +1,937 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "intersectionPatchToPatch.H" +#include "triIntersect.H" +#include "vtkWritePolyData.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace patchToPatches +{ + defineTypeNameAndDebug(intersection, 0); + addToRunTimeSelectionTable(patchToPatch, intersection, bool); + + int intersection::debugSrcFacei = + debug::debugSwitch((intersection::typeName + "SrcFace").c_str(), -1); + int intersection::debugTgtFacei = + debug::debugSwitch((intersection::typeName + "TgtFace").c_str(), -1); +} +} + + +// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * // + +template +Foam::FixedList +Foam::patchToPatches::intersection::triPointValues +( + const triFace& t, + const UList& values +) +{ + FixedList result; + forAll(t, i) + { + result[i] = values[t[i]]; + } + return result; +} + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::treeBoundBox Foam::patchToPatches::intersection::srcBoxStatic +( + const face& srcFace, + const pointField& srcPoints, + const vectorField& srcPointNormals +) +{ + static DynamicList ps; + + ps.clear(); + + const scalar l = sqrt(mag(srcFace.area(srcPoints))); + + forAll(srcFace, srcFacePointi) + { + const label srcPointi = srcFace[srcFacePointi]; + + const point& p = srcPoints[srcPointi]; + const vector& n = srcPointNormals[srcPointi]; + + ps.append(p - l/2*n); + ps.append(p + l/2*n); + } + + return treeBoundBox(ps); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::treeBoundBox Foam::patchToPatches::intersection::srcBox +( + const face& srcFace, + const pointField& srcPoints, + const vectorField& srcPointNormals +) const +{ + return srcBoxStatic(srcFace, srcPoints, srcPointNormals); +} + + +bool Foam::patchToPatches::intersection::intersectFaces +( + const primitiveOldTimePatch& srcPatch, + const vectorField& srcPointNormals, + const vectorField& srcPointNormals0, + const primitiveOldTimePatch& tgtPatch, + const label srcFacei, + const label tgtFacei +) +{ + // Quick rejection based on bound box + const treeBoundBox srcFaceBox = + srcBox + ( + srcPatch.localFaces()[srcFacei], + srcPatch.localPoints(), + srcPointNormals + ); + const treeBoundBox tgtFaceBox(tgtPatch.points(), tgtPatch[tgtFacei]); + if (!srcFaceBox.overlaps(tgtFaceBox)) return false; + + // Construct face triangulations on demand + if (srcTriPoints_[srcFacei].empty()) + { + triEngine_.triangulate + ( + UIndirectList + ( + srcPatch.localPoints(), + srcPatch.localFaces()[srcFacei] + ) + ); + + srcTriPoints_[srcFacei] = + triEngine_.triPoints(srcPatch.localFaces()[srcFacei]); + srcTriFaceEdges_[srcFacei] = triEngine_.triEdges(); + } + if (tgtTriPoints_[tgtFacei].empty()) + { + triEngine_.triangulate + ( + UIndirectList + ( + tgtPatch.localPoints(), + tgtPatch.localFaces()[tgtFacei] + ) + ); + + tgtTriPoints_[tgtFacei] = + triEngine_.triPoints(tgtPatch.localFaces()[tgtFacei]); + tgtTriFaceEdges_[tgtFacei] = triEngine_.triEdges(); + } + + // Construct and initialise workspace + bool srcCouples = false; + couple srcCouple; + srcFaceEdgePart_.resize(srcPatch[srcFacei].size()); + forAll(srcFaceEdgePart_, srcFaceEdgei) + { + const edge e = + srcPatch.localFaces()[srcFacei].faceEdge(srcFaceEdgei); + const vector eC = e.centre(srcPatch.localPoints()); + srcFaceEdgePart_[srcFaceEdgei] = part(Zero, eC); + } + + bool tgtCouples = false; + couple tgtCouple; + tgtFaceEdgePart_.resize(tgtPatch[tgtFacei].size()); + forAll(tgtFaceEdgePart_, tgtFaceEdgei) + { + const edge e = + tgtPatch.localFaces()[tgtFacei].faceEdge(tgtFaceEdgei); + const vector eC = e.centre(tgtPatch.localPoints()); + tgtFaceEdgePart_[tgtFaceEdgei] = part(Zero, eC); + } + + part errorPart(Zero, srcPatch.faceCentres()[srcFacei]); + + // Cache the face area magnitudes + const scalar srcMagA = mag(srcPatch.faceAreas()[srcFacei]); + const scalar tgtMagA = mag(tgtPatch.faceAreas()[tgtFacei]); + + // Determine whether or not to debug this tri intersection + const bool debugTriIntersect = + (debugSrcFacei != -1 || debugTgtFacei != -1) + && (debugSrcFacei == -1 || debugSrcFacei == srcFacei) + && (debugTgtFacei == -1 || debugTgtFacei == tgtFacei); + + // Loop the face triangles and compute the intersections + bool anyCouples = false; + forAll(srcTriPoints_[srcFacei], srcFaceTrii) + { + const triFace& srcT = srcTriPoints_[srcFacei][srcFaceTrii]; + + const FixedList srcPs = + triPointValues(srcT, srcPatch.localPoints()); + const FixedList srcNs = + triPointValues(srcT, srcPointNormals); + + forAll(tgtTriPoints_[tgtFacei], tgtFaceTrii) + { + const triFace tgtT = + reverse() + ? tgtTriPoints_[tgtFacei][tgtFaceTrii].reverseFace() + : tgtTriPoints_[tgtFacei][tgtFaceTrii]; + + const FixedList tgtPs = + triPointValues(tgtT, tgtPatch.localPoints()); + + // Do tri-intersection + ictSrcPoints_.clear(); + ictSrcPointNormals_.clear(); + ictTgtPoints_.clear(); + ictPointLocations_.clear(); + triIntersect::intersectTris + ( + srcPs, + srcNs, + {false, false, false}, + {-1, -1, -1}, + tgtPs, + {false, false, false}, + {-1, -1, -1}, + ictSrcPoints_, + ictSrcPointNormals_, + ictTgtPoints_, + ictPointLocations_, + debugTriIntersect, + debugTriIntersect + ? word + ( + typeName + + "_srcFace=" + Foam::name(srcFacei) + + "_tgtFace=" + Foam::name(tgtFacei) + + "_intersection=" + Foam::name + (srcFaceTrii*tgtTriPoints_[tgtFacei].size() + tgtFaceTrii) + ) + : word::null + ); + + // If there is no intersection then continue + if (ictPointLocations_.empty()) + { + continue; + } + + // Mark that there has been an intersection + anyCouples = true; + + // Compute the intersection geometry + const part ictSrcPart(ictSrcPoints_); + const part ictTgtPart(ictTgtPoints_); + + // If the intersection is below tolerance then continue + if + ( + mag(ictSrcPart.area) < small*srcMagA + || mag(ictTgtPart.area) < small*tgtMagA + ) + { + continue; + } + + // Mark that the source and target faces intersect + srcCouples = tgtCouples = true; + + // Store the intersection geometry + srcCouple += ictSrcPart; + srcCouple.nbr += ictTgtPart; + if (reverse()) + { + tgtCouple += ictTgtPart; + tgtCouple.nbr += ictSrcPart; + } + else + { + tgtCouple -= ictTgtPart; + tgtCouple.nbr -= ictSrcPart; + } + + // Store the intersection polygons for debugging + const label debugSrcPoint0 = debugPoints_.size(); + const label debugTgtPoint0 = + debugPoints_.size() + ictSrcPoints_.size(); + if (debug) + { + debugPoints_.append(ictSrcPoints_); + debugPoints_.append(ictTgtPoints_); + debugFaces_.append + ( + debugSrcPoint0 + identity(ictSrcPoints_.size()) + ); + debugFaceSrcFaces_.append(srcFacei); + debugFaceTgtFaces_.append(tgtFacei); + debugFaceSides_.append(1); + debugFaces_.append + ( + debugTgtPoint0 + identity(ictTgtPoints_.size()) + ); + debugFaceSrcFaces_.append(srcFacei); + debugFaceTgtFaces_.append(tgtFacei); + debugFaceSides_.append(-1); + } + + // Store edge and error areas + forAll(ictPointLocations_, i0) + { + const label i1 = ictPointLocations_.fcIndex(i0); + + // Get the locations on each end of this edge of the + // intersection polygon + const triIntersect::location l0 = ictPointLocations_[i0]; + const triIntersect::location l1 = ictPointLocations_[i1]; + + // Get the geometry for the projection of this edge + const part ictEdgePart + ( + FixedList + ({ + ictSrcPoints_[i0], + ictSrcPoints_[i1], + ictTgtPoints_[i1], + ictTgtPoints_[i0] + }) + ); + + // Store the "side" of the intersection that this edge + // corresponds to + label ictEdgeSide = -labelMax; + + // If this edge corresponds to an edge of the source + // triangle + if + ( + l0.isSrcNotTgtPoint() + || l1.isSrcNotTgtPoint() + || ( + l0.isIntersection() + && l1.isIntersection() + && l0.srcEdgei() == l1.srcEdgei() + ) + ) + { + const label srcEi = + l0.isSrcPoint() ? l0.srcPointi() + : l1.isSrcPoint() ? (l1.srcPointi() + 2) % 3 + : l0.srcEdgei(); + + const label srcFaceEdgei = + srcTriFaceEdges_[srcFacei][srcFaceTrii][srcEi]; + + if (srcFaceEdgei < srcPatch[srcFacei].size()) + { + srcFaceEdgePart_[srcFaceEdgei] += ictEdgePart; + ictEdgeSide = 1; + } + else + { + errorPart += ictEdgePart; + ictEdgeSide = 0; + } + } + + // If this edge corresponds to an edge of the target + // triangle + else if + ( + l0.isTgtNotSrcPoint() + || l1.isTgtNotSrcPoint() + || ( + l0.isIntersection() + && l1.isIntersection() + && l0.tgtEdgei() == l1.tgtEdgei() + ) + ) + { + const label tgtEi = + l0.isTgtPoint() ? (l0.tgtPointi() + 2) % 3 + : l1.isTgtPoint() ? l1.tgtPointi() + : l0.tgtEdgei(); + + const label tgtFaceEdgei = + tgtTriFaceEdges_[tgtFacei][tgtFaceTrii] + [reverse() ? 2 - tgtEi : tgtEi]; + + if (tgtFaceEdgei < tgtPatch[tgtFacei].size()) + { + if (reverse()) + { + tgtFaceEdgePart_[tgtFaceEdgei] += ictEdgePart; + } + else + { + tgtFaceEdgePart_[tgtFaceEdgei] -= ictEdgePart; + } + ictEdgeSide = -1; + } + else + { + errorPart += ictEdgePart; + ictEdgeSide = 0; + } + } + + // No other location combinations should be possible for an + // intersection without any shared points + else + { + FatalErrorInFunction + << "Tri-intersection topology not recognised. " + << "This is a bug." << exit(FatalError); + } + + // Store the projected edge quadrilateral for debugging + if (debug) + { + debugFaces_.append + ( + labelList + ({ + debugSrcPoint0 + i0, + debugSrcPoint0 + i1, + debugTgtPoint0 + i1, + debugTgtPoint0 + i0 + }) + ); + debugFaceSrcFaces_.append(srcFacei); + debugFaceTgtFaces_.append(tgtFacei); + debugFaceSides_.append(ictEdgeSide); + } + } + } + } + + // If the source face couples the target, then store the intersection + if (srcCouples) + { + srcLocalTgtFaces_[srcFacei].append(tgtFacei); + srcCouples_[srcFacei].append(srcCouple); + } + + // If any intersection has occurred then store the edge and error parts + if (anyCouples) + { + forAll(srcFaceEdgeParts_[srcFacei], srcFaceEdgei) + { + srcFaceEdgeParts_[srcFacei][srcFaceEdgei] += + srcFaceEdgePart_[srcFaceEdgei]; + } + srcErrorParts_[srcFacei] -= sum(tgtFaceEdgePart_); + srcErrorParts_[srcFacei] += errorPart; + } + + // If the target face couples the source, then store in the intersection + if (tgtCouples) + { + tgtLocalSrcFaces_[tgtFacei].append(srcFacei); + tgtCouples_[tgtFacei].append(tgtCouple); + } + + return anyCouples; +} + + +void Foam::patchToPatches::intersection::initialise +( + const primitiveOldTimePatch& srcPatch, + const vectorField& srcPointNormals, + const vectorField& srcPointNormals0, + const primitiveOldTimePatch& tgtPatch +) +{ + patchToPatch::initialise + ( + srcPatch, + srcPointNormals, + srcPointNormals0, + tgtPatch + ); + + srcCouples_.resize(srcPatch.size()); + forAll(srcLocalTgtFaces_, i) + { + srcCouples_[i].clear(); + } + + srcEdgeParts_.resize(srcPatch.nEdges()); + forAll(srcEdgeParts_, srcEdgei) + { + const edge& e = srcPatch.edges()[srcEdgei]; + const point c = e.centre(srcPatch.localPoints()); + srcEdgeParts_[srcEdgei] = part(Zero, c); + } + + srcErrorParts_.resize(srcPatch.size()); + forAll(srcErrorParts_, srcFacei) + { + srcErrorParts_[srcFacei] = + part(Zero, srcPatch.faceCentres()[srcFacei]); + } + + tgtCouples_.resize(tgtPatch.size()); + forAll(tgtLocalSrcFaces_, i) + { + tgtCouples_[i].clear(); + } + + srcTriPoints_ = List(srcPatch.size()); + srcTriFaceEdges_ = List>>(srcPatch.size()); + tgtTriPoints_ = List(tgtPatch.size()); + tgtTriFaceEdges_ = List>>(tgtPatch.size()); + + srcFaceEdgeParts_.resize(srcPatch.size()); + forAll(srcFaceEdgeParts_, srcFacei) + { + srcFaceEdgeParts_[srcFacei].resize(srcPatch[srcFacei].size()); + forAll(srcFaceEdgeParts_[srcFacei], srcFaceEdgei) + { + const label srcEdgei = + srcPatch.faceEdges()[srcFacei][srcFaceEdgei]; + srcFaceEdgeParts_[srcFacei][srcFaceEdgei] = srcEdgeParts_[srcEdgei]; + } + } + + if (debug) + { + debugPoints_.clear(); + debugFaces_.clear(); + debugFaceSrcFaces_.clear(); + debugFaceTgtFaces_.clear(); + debugFaceSides_.clear(); + } +} + + +void Foam::patchToPatches::intersection::rDistributeTgt +( + const primitiveOldTimePatch& tgtPatch, + const distributionMap& tgtMap +) +{ + patchToPatch::rDistributeTgt(tgtPatch, tgtMap); + + rDistributeListList(tgtPatch.size(), tgtMap, tgtCouples_); +} + + +Foam::label Foam::patchToPatches::intersection::finalise +( + const primitiveOldTimePatch& srcPatch, + const vectorField& srcPointNormals, + const vectorField& srcPointNormals0, + const primitiveOldTimePatch& tgtPatch, + const transformer& tgtToSrc +) +{ + const label nCouples = + patchToPatch::finalise + ( + srcPatch, + srcPointNormals, + srcPointNormals0, + tgtPatch, + tgtToSrc + ); + + // Convert face-edge-parts to edge-parts + labelList srcEdgeNParts(srcEdgeParts_.size(), 0); + forAll(srcEdgeParts_, srcEdgei) + { + const edge& e = srcPatch.edges()[srcEdgei]; + + srcEdgeParts_[srcEdgei] = part(); + + forAll(srcPatch.edgeFaces()[srcEdgei], i) + { + const label srcFacei = srcPatch.edgeFaces()[srcEdgei][i]; + const label srcFaceEdgei = + findIndex(srcPatch.faceEdges()[srcFacei], srcEdgei); + + const edge fe = + srcPatch.localFaces()[srcFacei].faceEdge(srcFaceEdgei); + + if (edge::compare(e, fe) > 0) + { + srcEdgeParts_[srcEdgei] += + srcFaceEdgeParts_[srcFacei][srcFaceEdgei]; + } + else + { + srcEdgeParts_[srcEdgei] -= + srcFaceEdgeParts_[srcFacei][srcFaceEdgei]; + } + + srcEdgeNParts[srcEdgei] ++; + } + } + forAll(srcEdgeParts_, srcEdgei) + { + srcEdgeParts_[srcEdgei].area /= srcEdgeNParts[srcEdgei]; + } + + // Add the difference between the face-edge-part and the edge-part into the + // face-error-parts + forAll(srcEdgeParts_, srcEdgei) + { + const edge& e = srcPatch.edges()[srcEdgei]; + + forAll(srcPatch.edgeFaces()[srcEdgei], i) + { + const label srcFacei = srcPatch.edgeFaces()[srcEdgei][i]; + const label srcFaceEdgei = + findIndex(srcPatch.faceEdges()[srcFacei], srcEdgei); + + const edge fe = + srcPatch.localFaces()[srcFacei].faceEdge(srcFaceEdgei); + + if (edge::compare(e, fe) > 0) + { + srcErrorParts_[srcFacei] -= srcEdgeParts_[srcEdgei]; + } + else + { + srcErrorParts_[srcFacei] += srcEdgeParts_[srcEdgei]; + } + + srcErrorParts_[srcFacei] += + srcFaceEdgeParts_[srcFacei][srcFaceEdgei]; + } + } + + // Transform the target couples back to the target side + if (!isNull(tgtToSrc)) + { + forAll(tgtCouples_, tgtFacei) + { + forAll(tgtCouples_[tgtFacei], i) + { + couple& c = tgtCouples_[tgtFacei][i]; + + c.area = tgtToSrc.invTransform(c.area); + c.centre = tgtToSrc.invTransformPosition(c.centre); + c.nbr.area = tgtToSrc.invTransform(c.nbr.area); + c.nbr.centre = tgtToSrc.invTransformPosition(c.nbr.centre); + } + } + } + + // Clear the triangulation workspace + srcTriPoints_.clear(); + srcTriFaceEdges_.clear(); + tgtTriPoints_.clear(); + tgtTriFaceEdges_.clear(); + + // Clear face-edge-parts + srcFaceEdgePart_.clear(); + tgtFaceEdgePart_.clear(); + srcFaceEdgeParts_.clear(); + + // Checking and reporting + if (nCouples != 0) + { + scalar srcArea = 0, srcCoupleArea = 0; + scalarField srcCoverage(srcPatch.size()); + scalarField srcOpenness(srcPatch.size()); + scalarField srcError(srcPatch.size()); + scalarField srcDepth(srcPatch.size()); + scalarField srcAngle(srcPatch.size()); + forAll(srcPatch, srcFacei) + { + const vector& a = srcPatch.faceAreas()[srcFacei]; + const scalar magA = mag(a); + const point& c = srcPatch.faceCentres()[srcFacei]; + + couple Cpl(part(Zero, c), part(Zero, c)); + forAll(srcCouples_[srcFacei], srcTgtFacei) + { + const couple& cpl = srcCouples_[srcFacei][srcTgtFacei]; + + Cpl += cpl; + Cpl.nbr += cpl.nbr; + } + const scalar magCA = mag(Cpl.area); + + srcArea += magA; + srcCoupleArea += magCA; + srcCoverage[srcFacei] = magCA/magA; + + vector projectionA = Zero; + scalar projectionV = 0; + forAll(srcCouples_[srcFacei], srcTgtFacei) + { + const couple& cpl = srcCouples_[srcFacei][srcTgtFacei]; + + projectionA += cpl.nbr.area; + projectionV += + - (cpl.area/3 & (cpl.centre - Cpl.centre)) + + (cpl.nbr.area/3 & (cpl.nbr.centre - Cpl.centre)); + } + forAll(srcPatch.faceEdges()[srcFacei], srcFaceEdgei) + { + const label srcEdgei = + srcPatch.faceEdges()[srcFacei][srcFaceEdgei]; + + const edge& e = srcPatch.edges()[srcEdgei]; + const edge fe = + srcPatch.localFaces()[srcFacei].faceEdge(srcFaceEdgei); + + const scalar sign = edge::compare(e, fe); + + projectionA += sign*srcEdgeParts_[srcEdgei].area; + projectionV += + sign*srcEdgeParts_[srcEdgei].area/3 + & (srcEdgeParts_[srcEdgei].centre - Cpl.centre); + } + projectionA += srcErrorParts_[srcFacei].area; + projectionV += + srcErrorParts_[srcFacei].area/3 + & (srcErrorParts_[srcFacei].centre - Cpl.centre); + + const vector aHat = normalised(a); + const vector aOppHat = normalised(a - Cpl.area + Cpl.nbr.area); + srcAngle[srcFacei] = + radToDeg(acos(min(max(aHat & aOppHat, -1), +1))); + srcOpenness[srcFacei] = mag(projectionA - Cpl.area)/magA; + srcError[srcFacei] = mag(srcErrorParts_[srcFacei].area)/magA; + srcDepth[srcFacei] = mag(projectionV)/pow3(sqrt(magA)); + } + + reduce(srcArea, sumOp()); + reduce(srcCoupleArea, sumOp()); + + scalar tgtArea = 0, tgtCoupleArea = 0; + scalarField tgtCoverage(tgtPatch.size()); + forAll(tgtPatch, tgtFacei) + { + const scalar magA = mag(tgtPatch.faceAreas()[tgtFacei]); + + vector aCouple = Zero; + forAll(tgtCouples_[tgtFacei], tgtSrcFacei) + { + aCouple += tgtCouples_[tgtFacei][tgtSrcFacei].area; + } + const scalar magACouple = mag(aCouple); + + tgtArea += magA; + tgtCoupleArea += magACouple; + tgtCoverage[tgtFacei] = magACouple/magA; + } + + reduce(tgtArea, sumOp()); + reduce(tgtCoupleArea, sumOp()); + + Info<< indent << "Source min/average/max coverage = " + << gMin(srcCoverage) << '/' << srcCoupleArea/srcArea << '/' + << gMax(srcCoverage) << endl + << indent << "Target min/average/max coverage = " + << gMin(tgtCoverage) << '/' << tgtCoupleArea/tgtArea << '/' + << gMax(tgtCoverage) << endl + << indent << "Source average openness/error/depth/angle = " + << gAverage(srcOpenness) << '/' << gAverage(srcError) << '/' + << gAverage(srcDepth) << '/' << gAverage(srcAngle) << endl + << indent << "Source max openness/error/depth/angle = " + << gMax(srcOpenness) << '/' << gMax(srcError) << '/' + << gMax(srcDepth) << '/' << gMax(srcAngle) << endl; + + if (debug) + { + word name = patchToPatch::typeName + '_' + typeName; + + if (Pstream::parRun()) + { + name += "_proc" + Foam::name(Pstream::myProcNo()); + } + + Info<< indent << "Writing intersected faces to " + << name + ".vtk" << endl; + vtkWritePolyData::write + ( + name + ".vtk", + name, + false, + debugPoints_, + labelList(), + labelListList(), + debugFaces_, + "srcFace", false, Field