2312/enrichedPatchCutFaces_8C_source.html

 /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
     Copyright (C) 2017-2022 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/>.

 Description
     Calculating cut faces of the enriched patch, together with the addressing
     into master and slave patch.

 \*---------------------------------------------------------------------------*/

 #include "enrichedPatch.H"
 #include "boolList.H"
 #include "CircularBuffer.H"
 #include "DynamicList.H"
 #include "labelPair.H"
 #include "primitiveMesh.H"
 #include "edgeHashes.H"

 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

 // If the cut face gets more than this number of points, it will be checked
 const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;


 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

 // Index of debug signs:
 // x - skip a point
 // l - left turn
 // r - right turn

 void Foam::enrichedPatch::calcCutFaces() const
 {
     if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
     {
         FatalErrorInFunction
             << "Cut faces addressing already calculated."
             << abort(FatalError);
     }

     const faceList& enFaces = enrichedFaces();
     const labelList& mp = meshPoints();
     const faceList& lf = localFaces();
     const pointField& lp = localPoints();
     const labelListList& pp = pointPoints();
     // Pout<< "enFaces: " << enFaces << endl;
     // Pout<< "lf: " << lf << endl;
     // Pout<< "lp: " << lp << endl;
     // Pout<< "pp: " << pp << endl;
     const Map<labelList>& masterPointFaceAddr = masterPointFaces();

     // Prepare the storage
     DynamicList<face> cf(2*lf.size());
     DynamicList<label> cfMaster(2*lf.size());
     DynamicList<label> cfSlave(2*lf.size());

     // Algorithm
     // Go through all the faces
     // 1) For each face, start at point zero and grab the edge zero.
     //    Both points are added into the cut face.  Mark the first edge
     //    as used and position the current point as the end of the first
     //    edge and previous point as point zero.
     // 2) Grab all possible points from the current point.  Excluding
     //    the previous point find the point giving the biggest right
     //    turn. Add the point to the face and mark the edges as used.
     //    Continue doing this until the face is complete, i.e. the next point
     //    to add is the first point of the face.
     // 3) Once the facet is completed, register its owner from the face
     //    it has been created from (remember that the first lot of faces
     //    inserted into the enriched faces list are the slaves, to allow
     //    matching of the other side).
     // 4) If the facet is created from an original slave face, go
     //    through the master patch and try to identify the master face
     //    this facet belongs to.  This is recognised by the fact that the
     //    faces consists exclusively of the points of the master face and
     //    the points projected onto the face.

     // Create a set of edge usage parameters
     edgeHashSet edgesUsedOnce(pp.size());
     edgeHashSet edgesUsedTwice(pp.size()*primitiveMesh::edgesPerPoint_);


     DynamicList<bool> usedFaceEdges(256);
     CircularBuffer<edge> edgeSeeds(256);

     forAll(lf, facei)
     {
         const face& curLocalFace = lf[facei];
         const face& curGlobalFace = enFaces[facei];

         // Pout<< "Doing face " << facei
         //     << " local: " << curLocalFace
         //     << " or " << curGlobalFace
         //     << endl;

         // if (facei < slavePatch_.size())
         // {
         //     Pout<< "original slave: " << slavePatch_[facei]
         //         << " local: " << slavePatch_.localFaces()[facei] << endl;
         // }
         // else
         // {
         //     Pout<< "original master: "
         //         << masterPatch_[facei - slavePatch_.size()] << " "
         //         << masterPatch_.localFaces()[facei - slavePatch_.size()]
         //         << endl;
         // }
         // {
         //     pointField facePoints = curLocalFace.points(lp);
         //     forAll(curLocalFace, pointi)
         //     {
         //         Pout<< "v " << facePoints[pointi].x() << " "
         //             << facePoints[pointi].y() << " "
         //             << facePoints[pointi].z() << endl;
         //     }
         // }

         // Track the usage of face edges.  When all edges are used, the
         // face decomposition is complete.
         // The seed edges include all the edges of the original face + all edges
         // of other faces that have been used in the construction of the
         // facet.  Edges from other faces can be considered as
         // internal to the current face if used only once.

         // Track the edge usage to avoid duplicate faces and reset it to unused
         usedFaceEdges.resize_nocopy(curLocalFace.size());
         usedFaceEdges = false;

         // Add edges of current face into the seed list.
         edgeSeeds.clear();
         edgeSeeds.push_back(curLocalFace.edges());

         // Grab face normal
         const vector normal = curLocalFace.unitNormal(lp);

         while (!edgeSeeds.empty())
         {
             // Pout<< "edgeSeeds.size(): " << edgeSeeds.size() << endl;

             const edge curEdge = edgeSeeds.front();
             edgeSeeds.pop_front();

             // Locate the edge in current face
             const label curEdgeWhich = curLocalFace.which(curEdge.start());

             // Check if the edge is in current face and if it has been
             // used already.  If so, skip it
             if
             (
                 curEdgeWhich > -1
              && curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
             )
             {
                 // Edge is in the starting face.
                 // If unused, mark it as used; if used, skip it
                 if (usedFaceEdges[curEdgeWhich]) continue;

                 usedFaceEdges[curEdgeWhich] = true;
             }

             // If the edge has already been used twice, skip it
             if (edgesUsedTwice.found(curEdge)) continue;

             // Pout<< "Trying new edge (" << mp[curEdge.start()]
             //     << ", " << mp[curEdge.end()]
             //     << ") seed: " << curEdge
             //     << " used: " << edgesUsedTwice.found(curEdge)
             //     << endl;

             // Estimate the size of cut face as twice the size of original face
             DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
             DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());

             // Found unused edge.
             label prevPointLabel = curEdge.start();
             cutFaceGlobalPoints.append(mp[prevPointLabel]);
             cutFaceLocalPoints.append(prevPointLabel);
             // Pout<< "prevPointLabel: " << mp[prevPointLabel] << endl;

             // Grab current point and append it to the list
             label curPointLabel = curEdge.end();
             point curPoint = lp[curPointLabel];
             cutFaceGlobalPoints.append(mp[curPointLabel]);
             cutFaceLocalPoints.append(curPointLabel);

             bool completedCutFace = false;

             label faceSizeDebug = maxFaceSizeDebug_;

             do
             {
                 // Grab the next point options

                 // Pout<< "curPointLabel: " << mp[curPointLabel] << endl;

                 const labelList& nextPoints = pp[curPointLabel];

                 // Pout<< "nextPoints: "
                 //     << labelUIndList(mp, nextPoints)
                 //     << endl;

                 // Get the vector along the edge and the right vector
                 vector ahead = curPoint - lp[prevPointLabel];
                 ahead.removeCollinear(normal);
                 ahead.normalise();

                 const vector right = normalised(normal ^ ahead);

                 // Pout<< "normal: " << normal
                 //     << " ahead: " << ahead
                 //     << " right: " << right
                 //     << endl;

                 scalar atanTurn = -GREAT;
                 label bestAtanPoint = -1;

                 forAll(nextPoints, nextI)
                 {
                     // Exclude the point we are coming from; there will always
                     // be more than one edge, so this is safe
                     if (nextPoints[nextI] != prevPointLabel)
                     {
                         // Pout<< "cur point: " << curPoint
                         //     << " trying for point: "
                         //     << mp[nextPoints[nextI]]
                         //     << " " << lp[nextPoints[nextI]];
                         vector newDir = lp[nextPoints[nextI]] - curPoint;
                         // Pout<< " newDir: " << newDir
                         //     << " mag: " << mag(newDir) << flush;
                         newDir.removeCollinear(normal);
                         scalar magNewDir = mag(newDir);
                         // Pout<< " corrected: " << newDir
                         //     << " mag: " << mag(newDir) << flush;

                         if (magNewDir < SMALL)
                         {
                             FatalErrorInFunction
                                 << "projection error: slave patch probably "
                                 << "does not project onto master.  "
                                 << "Please switch on "
                                 << "enriched patch debug for more info"
                                 << abort(FatalError);
                         }

                         newDir /= magNewDir;

                         const scalar curAtanTurn =
                             atan2(newDir & right, newDir & ahead);

                         // Pout<< " atan: " << curAtanTurn << endl;

                         if (curAtanTurn > atanTurn)
                         {
                             bestAtanPoint = nextPoints[nextI];
                             atanTurn = curAtanTurn;
                         }
                     } // end of prev point skip
                 } // end of next point selection

                 // Pout<< "   bestAtanPoint: " << bestAtanPoint << " or "
                 //     << mp[bestAtanPoint]
                 //     << endl;

                 // Selected next best point.

                 // Pout<< "cutFaceGlobalPoints: "
                 //     << cutFaceGlobalPoints
                 //     << endl;

                 // Check if the edge about to be added has been used
                 // in the current face or twice in other faces.  If
                 // so, the face is bad.
                 const label whichNextPoint = curLocalFace.which(curPointLabel);

                 if
                 (
                     whichNextPoint > -1
                  && curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
                  && usedFaceEdges[whichNextPoint]
                 )
                 {
                     // This edge is already used in current face
                     // face cannot be good; start on a new one

                     // Pout<< "Double usage in current face, cannot be good"
                     //     << endl;

                     completedCutFace = true;
                 }


                 if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
                 {
                     // This edge is already used -
                     // face cannot be good; start on a new one
                     completedCutFace = true;

                     // Pout<< "Double usage elsewhere, cannot be good" << endl;
                 }

                 if (completedCutFace) continue;

                 // Insert the next best point into the list
                 if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
                 {
                     // Face is completed.  Save it and move on to the next
                     // available edge
                     completedCutFace = true;

                     if (debug)
                     {
                         Pout<< " local: " << cutFaceLocalPoints
                             << " one side: " << facei;
                     }

                     // Append the face
                     face cutFaceGlobal;
                     cutFaceGlobal.transfer(cutFaceGlobalPoints);

                     cf.append(cutFaceGlobal);

                     face cutFaceLocal;
                     cutFaceLocal.transfer(cutFaceLocalPoints);

                     // Pout<< "\ncutFaceLocal: "
                     //     << cutFaceLocal.points(lp)
                     //     << endl;

                     // Go through all edges of the cut faces.
                     // If the edge corresponds to a starting face edge,
                     // mark the starting face edge as true

                     forAll(cutFaceLocal, cuti)
                     {
                         const edge curCutFaceEdge
                         (
                             cutFaceLocal[cuti],
                             cutFaceLocal.nextLabel(cuti)
                         );

                         // Increment the usage count using two hash sets
                         auto euoIter = edgesUsedOnce.find(curCutFaceEdge);

                         if (!euoIter.good())
                         {
                             // Pout<< "Found edge not used before: "
                             //     << curCutFaceEdge
                             //     << endl;
                             edgesUsedOnce.insert(curCutFaceEdge);
                         }
                         else
                         {
                             // Pout<< "Found edge used once: "
                             //     << curCutFaceEdge
                             //     << endl;
                             edgesUsedOnce.erase(euoIter);
                             edgesUsedTwice.insert(curCutFaceEdge);
                         }

                         const label curCutFaceEdgeWhich = curLocalFace.which
                         (
                             curCutFaceEdge.start()
                         );

                         if
                         (
                             curCutFaceEdgeWhich > -1
                          && curLocalFace.nextLabel(curCutFaceEdgeWhich)
                          == curCutFaceEdge.end()
                         )
                         {
                             // Found edge in original face

                             // Pout<< "Found edge in orig face: "
                             //     << curCutFaceEdge << ": "
                             //     << curCutFaceEdgeWhich
                             //     << endl;

                             usedFaceEdges[curCutFaceEdgeWhich] = true;
                         }
                         else
                         {
                             // Edge not in original face.  Add it to seeds

                             // Pout<< "Found new edge seed: "
                             //     << curCutFaceEdge
                             //     << endl;

                             edgeSeeds.push_back(curCutFaceEdge.reverseEdge());
                         }
                     }

                     // Find out what the other side is

                     // Algorithm

                     // If the face is in the slave half of the
                     // enrichedFaces lists, it may be matched against
                     // the master face.  It can be recognised by the
                     // fact that all its points belong to one master
                     // face.  For this purpose:
                     // 1) Grab the master faces around the point zero
                     // of the cut face and collect all master faces
                     // around it.
                     // 2) Loop through the rest of cut face points and
                     // if the point refers to one of the faces marked
                     // by point zero, increment its count.
                     // 3) When completed, the face whose count is
                     // equal to the number of points in the cut face
                     // is the other side.  If this is not the case, there is no
                     // face on the other side.

                     if (facei < slavePatch_.size())
                     {
                         const auto mpfAddrIter =
                             masterPointFaceAddr.cfind(cutFaceGlobal[0]);

                         bool otherSideFound = false;

                         if (mpfAddrIter.good())
                         {
                             bool miss = false;

                             // Create the label count using point zero
                             const labelList& masterFacesOfPZero = mpfAddrIter();

                             labelList hits(masterFacesOfPZero.size(), 1);

                             for
                             (
                                 label pointi = 1;
                                 pointi < cutFaceGlobal.size();
                                 pointi++
                             )
                             {
                                 const auto mpfAddrPointIter =
                                     masterPointFaceAddr.cfind
                                     (
                                         cutFaceGlobal[pointi]
                                     );

                                 if (!mpfAddrPointIter.good())
                                 {
                                     // Point is off the master patch. Skip
                                     miss = true;
                                     break;
                                 }

                                 const labelList& curMasterFaces =
                                     mpfAddrPointIter();

                                 // For every current face, try to find it in the
                                 // zero-list
                                 for (const label facei : curMasterFaces)
                                 {
                                     forAll(masterFacesOfPZero, j)
                                     {
                                         if (facei == masterFacesOfPZero[j])
                                         {
                                             hits[j]++;
                                             break;
                                         }
                                     }
                                 }
                             }

                             // If all point are found attempt matching
                             if (!miss)
                             {
                                 forAll(hits, pointi)
                                 {
                                     if (hits[pointi] == cutFaceGlobal.size())
                                     {
                                         // Found other side.
                                         otherSideFound = true;

                                         cfMaster.append
                                         (
                                             masterFacesOfPZero[pointi]
                                         );

                                         cfSlave.append(facei);

                                         // Reverse the face such that it
                                         // points out of the master patch
                                         cf.last().flip();

                                         if (debug)
                                         {
                                             Pout<< " other side: "
                                                 << masterFacesOfPZero[pointi]
                                                 << endl;
                                         }
                                     } // end of hits
                                 } // end of for all hits

                             } // end of not miss

                             if (!otherSideFound || miss)
                             {
                                 if (debug)
                                 {
                                     Pout<< " solo slave A" << endl;
                                 }

                                 cfMaster.append(-1);
                                 cfSlave.append(facei);
                             }
                         }
                         else
                         {
                             // First point not in master patch
                             if (debug)
                             {
                                 Pout<< " solo slave B" << endl;
                             }

                             cfMaster.append(-1);
                             cfSlave.append(facei);
                         }
                     }
                     else
                     {
                         if (debug)
                         {
                             Pout<< " master side" << endl;
                         }

                         cfMaster.append(facei - slavePatch_.size());
                         cfSlave.append(-1);
                     }
                 }
                 else
                 {
                     // Face not completed.  Prepare for the next point search
                     prevPointLabel = curPointLabel;
                     curPointLabel = bestAtanPoint;
                     curPoint = lp[curPointLabel];
                     cutFaceGlobalPoints.append(mp[curPointLabel]);
                     cutFaceLocalPoints.append(curPointLabel);

                     if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
                     {
                         faceSizeDebug *= 2;

                         // Check for duplicate points in the face
                         forAll(cutFaceGlobalPoints, checkI)
                         {
                             for
                             (
                                 label checkJ = checkI + 1;
                                 checkJ < cutFaceGlobalPoints.size();
                                 checkJ++
                             )
                             {
                                 if
                                 (
                                     cutFaceGlobalPoints[checkI]
                                  == cutFaceGlobalPoints[checkJ]
                                 )
                                 {
                                     // Shrink local points for debugging
                                     cutFaceLocalPoints.shrink();

                                     face origFace;
                                     face origFaceLocal;
                                     if (facei < slavePatch_.size())
                                     {
                                         origFace = slavePatch_[facei];
                                         origFaceLocal =
                                             slavePatch_.localFaces()[facei];
                                     }
                                     else
                                     {
                                         origFace =
                                             masterPatch_
                                             [facei - slavePatch_.size()];

                                         origFaceLocal =
                                             masterPatch_.localFaces()
                                             [facei - slavePatch_.size()];
                                     }

                                     FatalErrorInFunction
                                         << "Duplicate point found in cut face. "
                                         << "Error in the face cutting "
                                         << "algorithm for global face "
                                         << origFace << " local face "
                                         << origFaceLocal << nl
                                         << "Slave size: " << slavePatch_.size()
                                         << " Master size: "
                                         << masterPatch_.size()
                                         << " index: " << facei << ".\n"
                                         << "Face: " << curGlobalFace << nl
                                         << "Cut face: " << cutFaceGlobalPoints
                                         << " local: " << cutFaceLocalPoints
                                         << nl << "Points: "
                                         << face(cutFaceLocalPoints).points(lp)
                                         << abort(FatalError);
                                 }
                             }
                         }
                     } // end of debug
                 }
             } while (!completedCutFace);
         } // end of used edges

         if (debug)
         {
             Pout<< " Finished face " << facei << endl;
         }

     } // end of local faces

     // Re-pack lists into compact storage
     cutFacesPtr_.reset(new faceList(std::move(cf)));
     cutFaceMasterPtr_.reset(new labelList(std::move(cfMaster)));
     cutFaceSlavePtr_.reset(new labelList(std::move(cfSlave)));
 }


 void Foam::enrichedPatch::clearCutFaces()
 {
     cutFacesPtr_.reset(nullptr);
     cutFaceMasterPtr_.reset(nullptr);
     cutFaceSlavePtr_.reset(nullptr);
 }


 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

 const Foam::faceList& Foam::enrichedPatch::cutFaces() const
 {
     if (!cutFacesPtr_)
     {
         calcCutFaces();
     }

     return *cutFacesPtr_;
 }


 const Foam::labelList& Foam::enrichedPatch::cutFaceMaster() const
 {
     if (!cutFaceMasterPtr_)
     {
         calcCutFaces();
     }

     return *cutFaceMasterPtr_;
 }


 const Foam::labelList& Foam::enrichedPatch::cutFaceSlave() const
 {
     if (!cutFaceSlavePtr_)
     {
         calcCutFaces();
     }

     return *cutFaceSlavePtr_;
 }


 // ************************************************************************* //
Foam::PrimitivePatch::localFaces
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
Definition: PrimitivePatch.C:310

Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)

Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...

FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:598

labelPair.H

Foam::enrichedPatch::cutFaceMaster
const labelList & cutFaceMaster() const
Return cut face master list.
Definition: enrichedPatchCutFaces.C:659

Foam::List< face >

edgeHashes.H

Foam::nl
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:50

primitiveMesh.H

Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:531

Foam::edgeHashSet
HashSet< edge, Hash< edge > > edgeHashSet
A HashSet with edge for its key. Hashing (and ==) on an edge is symmetric.
Definition: edgeHashes.H:48

Foam::VectorSpace::end
iterator end() noexcept
Return an iterator to end of VectorSpace.
Definition: VectorSpaceI.H:199

Foam::enrichedPatch::cutFaces
const faceList & cutFaces() const
Return list of cut faces.
Definition: enrichedPatchCutFaces.C:648

Foam::Vector::normalise
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
Definition: VectorI.H:114

Foam::primitiveMesh::edgesPerPoint_
static const unsigned edgesPerPoint_
Estimated number of edges per point.
Definition: primitiveMesh.H:545

Foam::normalised
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:674

Foam::enrichedPatch::masterPointFaces
const Map< labelList > & masterPointFaces() const
Master point face addressing.
Definition: enrichedPatchMasterPoints.C:96

Foam::labelListList
List< labelList > labelListList
List of labelList.
Definition: labelList.H:38

boolList.H

forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:421

Foam::enrichedPatch::enrichedFaces
const faceList & enrichedFaces() const
Return enriched faces.
Definition: enrichedPatchFaces.C:380

Foam::faceList
List< face > faceList
List of faces.
Definition: faceListFwd.H:39

Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38

Foam::enrichedPatch::cutFaceSlave
const labelList & cutFaceSlave() const
Return cut face slave list.
Definition: enrichedPatchCutFaces.C:670

Foam::vector
Vector< scalar > vector
Definition: vector.H:57

Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:139

Foam::enrichedPatch::pointPoints
const labelListList & pointPoints() const
Return point-point addressing.
Definition: enrichedPatch.C:191

Foam::atan2
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
Definition: dimensionedScalar.C:305

Foam::ensightOutput::debug
int debug
Static debugging option.

Foam::enrichedPatch::localPoints
const pointField & localPoints() const
Return local points.
Definition: enrichedPatch.C:180

DynamicList.H

Foam::point
vector point
Point is a vector.
Definition: point.H:37

enrichedPatch.H

Foam::enrichedPatch::localFaces
const faceList & localFaces() const
Return local faces.
Definition: enrichedPatch.C:169

Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:62

Foam::enrichedPatch::meshPoints
const labelList & meshPoints() const
Return mesh points.
Definition: enrichedPatch.C:158

Foam::Vector::removeCollinear
Vector< Cmpt > & removeCollinear(const Vector< Cmpt > &unitVec)
Inplace removal of components that are collinear to the given unit vector.
Definition: VectorI.H:136

Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.

Foam::constant::atomic::mp
const dimensionedScalar mp
Proton mass.

pp
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())

CircularBuffer.H