44 const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;
54 void Foam::enrichedPatch::calcCutFaces()
const 56 if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
59 <<
"Cut faces addressing already calculated." 75 DynamicList<face> cf(2*lf.size());
76 DynamicList<label> cfMaster(2*lf.size());
77 DynamicList<label> cfSlave(2*lf.size());
105 DynamicList<bool> usedFaceEdges(256);
106 CircularBuffer<edge> edgeSeeds(256);
110 const face& curLocalFace = lf[facei];
111 const face& curGlobalFace = enFaces[facei];
148 usedFaceEdges.resize_nocopy(curLocalFace.size());
149 usedFaceEdges =
false;
153 edgeSeeds.push_back(curLocalFace.edges());
156 const vector normal = curLocalFace.unitNormal(lp);
158 while (!edgeSeeds.empty())
162 const edge curEdge = edgeSeeds.front();
163 edgeSeeds.pop_front();
166 const label curEdgeWhich = curLocalFace.which(curEdge.start());
173 && curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
178 if (usedFaceEdges[curEdgeWhich])
continue;
180 usedFaceEdges[curEdgeWhich] =
true;
184 if (edgesUsedTwice.found(curEdge))
continue;
193 DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
194 DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
197 label prevPointLabel = curEdge.start();
198 cutFaceGlobalPoints.append(
mp[prevPointLabel]);
199 cutFaceLocalPoints.append(prevPointLabel);
203 label curPointLabel = curEdge.
end();
204 point curPoint = lp[curPointLabel];
205 cutFaceGlobalPoints.append(
mp[curPointLabel]);
206 cutFaceLocalPoints.append(curPointLabel);
208 bool completedCutFace =
false;
210 label faceSizeDebug = maxFaceSizeDebug_;
225 vector ahead = curPoint - lp[prevPointLabel];
236 scalar atanTurn = -GREAT;
237 label bestAtanPoint = -1;
243 if (nextPoints[nextI] != prevPointLabel)
249 vector newDir = lp[nextPoints[nextI]] - curPoint;
253 scalar magNewDir =
mag(newDir);
257 if (magNewDir < SMALL)
260 <<
"projection error: slave patch probably " 261 <<
"does not project onto master. " 262 <<
"Please switch on " 263 <<
"enriched patch debug for more info" 269 const scalar curAtanTurn =
270 atan2(newDir & right, newDir & ahead);
274 if (curAtanTurn > atanTurn)
276 bestAtanPoint = nextPoints[nextI];
277 atanTurn = curAtanTurn;
295 const label whichNextPoint = curLocalFace.which(curPointLabel);
300 && curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
301 && usedFaceEdges[whichNextPoint]
310 completedCutFace =
true;
314 if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
318 completedCutFace =
true;
323 if (completedCutFace)
continue;
326 if (
mp[bestAtanPoint] == cutFaceGlobalPoints[0])
330 completedCutFace =
true;
334 Pout<<
" local: " << cutFaceLocalPoints
335 <<
" one side: " << facei;
340 cutFaceGlobal.transfer(cutFaceGlobalPoints);
342 cf.append(cutFaceGlobal);
345 cutFaceLocal.transfer(cutFaceLocalPoints);
355 forAll(cutFaceLocal, cuti)
357 const edge curCutFaceEdge
360 cutFaceLocal.nextLabel(cuti)
364 auto euoIter = edgesUsedOnce.find(curCutFaceEdge);
371 edgesUsedOnce.insert(curCutFaceEdge);
378 edgesUsedOnce.erase(euoIter);
379 edgesUsedTwice.insert(curCutFaceEdge);
382 const label curCutFaceEdgeWhich = curLocalFace.which
384 curCutFaceEdge.start()
389 curCutFaceEdgeWhich > -1
390 && curLocalFace.nextLabel(curCutFaceEdgeWhich)
391 == curCutFaceEdge.end()
401 usedFaceEdges[curCutFaceEdgeWhich] =
true;
411 edgeSeeds.push_back(curCutFaceEdge.reverseEdge());
435 if (facei < slavePatch_.size())
437 const auto mpfAddrIter =
438 masterPointFaceAddr.cfind(cutFaceGlobal[0]);
440 bool otherSideFound =
false;
442 if (mpfAddrIter.good())
447 const labelList& masterFacesOfPZero = mpfAddrIter();
449 labelList hits(masterFacesOfPZero.size(), 1);
454 pointi < cutFaceGlobal.size();
458 const auto mpfAddrPointIter =
459 masterPointFaceAddr.cfind
461 cutFaceGlobal[pointi]
464 if (!mpfAddrPointIter.good())
476 for (
const label facei : curMasterFaces)
478 forAll(masterFacesOfPZero, j)
480 if (facei == masterFacesOfPZero[j])
494 if (hits[pointi] == cutFaceGlobal.size())
497 otherSideFound =
true;
501 masterFacesOfPZero[pointi]
504 cfSlave.append(facei);
512 Pout<<
" other side: " 513 << masterFacesOfPZero[pointi]
521 if (!otherSideFound || miss)
529 cfSlave.append(facei);
541 cfSlave.append(facei);
551 cfMaster.append(facei - slavePatch_.size());
558 prevPointLabel = curPointLabel;
559 curPointLabel = bestAtanPoint;
560 curPoint = lp[curPointLabel];
561 cutFaceGlobalPoints.append(
mp[curPointLabel]);
562 cutFaceLocalPoints.append(curPointLabel);
564 if (
debug || cutFaceGlobalPoints.size() > faceSizeDebug)
569 forAll(cutFaceGlobalPoints, checkI)
573 label checkJ = checkI + 1;
574 checkJ < cutFaceGlobalPoints.size();
580 cutFaceGlobalPoints[checkI]
581 == cutFaceGlobalPoints[checkJ]
585 cutFaceLocalPoints.shrink();
589 if (facei < slavePatch_.size())
591 origFace = slavePatch_[facei];
599 [facei - slavePatch_.size()];
603 [facei - slavePatch_.size()];
607 <<
"Duplicate point found in cut face. " 608 <<
"Error in the face cutting " 609 <<
"algorithm for global face " 610 << origFace <<
" local face " 611 << origFaceLocal <<
nl 612 <<
"Slave size: " << slavePatch_.size()
614 << masterPatch_.size()
615 <<
" index: " << facei <<
".\n" 616 <<
"Face: " << curGlobalFace <<
nl 617 <<
"Cut face: " << cutFaceGlobalPoints
618 <<
" local: " << cutFaceLocalPoints
620 << face(cutFaceLocalPoints).points(lp)
627 }
while (!completedCutFace);
632 Pout<<
" Finished face " << facei <<
endl;
638 cutFacesPtr_.reset(
new faceList(std::move(cf)));
639 cutFaceMasterPtr_.reset(
new labelList(std::move(cfMaster)));
640 cutFaceSlavePtr_.reset(
new labelList(std::move(cfSlave)));
644 void Foam::enrichedPatch::clearCutFaces()
646 cutFacesPtr_.reset(
nullptr);
647 cutFaceMasterPtr_.reset(
nullptr);
648 cutFaceSlavePtr_.reset(
nullptr);
661 return *cutFacesPtr_;
667 if (!cutFaceMasterPtr_)
672 return *cutFaceMasterPtr_;
678 if (!cutFaceSlavePtr_)
683 return *cutFaceSlavePtr_;
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList & cutFaceMaster() const
Return cut face master list.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
HashSet< edge, Hash< edge > > edgeHashSet
A HashSet with edge for its key. Hashing (and ==) on an edge is symmetric.
iterator end() noexcept
Return an iterator to end of VectorSpace.
const faceList & cutFaces() const
Return list of cut faces.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
static const unsigned edgesPerPoint_
Estimated number of edges per point.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
const Map< labelList > & masterPointFaces() const
Master point face addressing.
List< labelList > labelListList
List of labelList.
#define forAll(list, i)
Loop across all elements in list.
const faceList & enrichedFaces() const
Return enriched faces.
List< face > faceList
List of faces.
vectorField pointField
pointField is a vectorField.
const labelList & cutFaceSlave() const
Return cut face slave list.
errorManip< error > abort(error &err)
const labelListList & pointPoints() const
Return point-point addressing.
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
int debug
Static debugging option.
const pointField & localPoints() const
Return local points.
vector point
Point is a vector.
const faceList & localFaces() const
Return local faces.
List< label > labelList
A List of labels.
const labelList & meshPoints() const
Return mesh points.
Vector< Cmpt > & removeCollinear(const Vector< Cmpt > &unitVec)
Inplace removal of components that are collinear to the given unit vector.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const dimensionedScalar mp
Proton mass.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())