50 partialFaceAreaWeightAMI,
111 label nFacesRemaining = srcAddr.size();
120 labelList seedFaces(nFacesRemaining, -1);
121 seedFaces[srcFacei] = tgtFacei;
124 bitSet mapFlag(nFacesRemaining,
true);
127 label startSeedi = 0;
132 bool continueWalk =
true;
138 visitedFaces.
clear();
156 mapFlag.
unset(srcFacei);
160 nonOverlapFaces.
append(srcFacei);
175 }
while (continueWalk);
183 const label srcFacei,
184 const label tgtStartFacei,
187 DynamicList<label>& nbrFaces,
189 DynamicList<label>& visitedFaces,
192 List<DynamicList<label>>& srcAddr,
193 List<DynamicList<scalar>>& srcWght,
194 List<DynamicList<point>>& srcCtr,
195 List<DynamicList<label>>& tgtAddr,
196 List<DynamicList<scalar>>& tgtWght
199 addProfiling(ami,
"faceAreaWeightAMI::processSourceFace");
201 if (tgtStartFacei == -1)
206 const auto& tgtPatch = this->tgtPatch();
209 nbrFaces.append(tgtStartFacei);
210 appendNbrFaces(tgtStartFacei, tgtPatch, visitedFaces, nbrFaces);
212 bool faceProcessed =
false;
214 label maxNeighbourFaces = nbrFaces.size();
216 while (!nbrFaces.empty())
219 label tgtFacei = nbrFaces.back();
221 visitedFaces.push_back(tgtFacei);
223 scalar interArea = 0;
225 calcInterArea(srcFacei, tgtFacei, interArea, interCentroid);
230 srcAddr[srcFacei].append(tgtFacei);
231 srcWght[srcFacei].append(interArea);
232 srcCtr[srcFacei].append(interCentroid);
234 tgtAddr[tgtFacei].append(srcFacei);
235 tgtWght[tgtFacei].append(interArea);
237 appendNbrFaces(tgtFacei, tgtPatch, visitedFaces, nbrFaces);
239 faceProcessed =
true;
241 maxNeighbourFaces =
max(maxNeighbourFaces, nbrFaces.size());
250 return faceProcessed;
259 const bitSet& mapFlag,
261 const DynamicList<label>& visitedFaces,
262 const bool errorOnNotFound
267 if (mapFlag.count() == 0)
273 const auto& srcPatch = this->srcPatch();
280 bool valuesSet =
false;
281 for (label faceS: srcNbrFaces)
283 if (mapFlag.test(faceS) && seedFaces[faceS] == -1)
285 for (label faceT : visitedFaces)
287 const scalar threshold =
291 if (overlaps(faceS, faceT, threshold))
293 seedFaces[faceS] = faceT;
313 label facei = startSeedi;
314 if (!mapFlag.test(startSeedi))
316 facei = mapFlag.find_next(facei);
318 const label startSeedi0 = facei;
320 bool foundNextSeed =
false;
326 foundNextSeed =
true;
329 if (seedFaces[facei] != -1)
332 tgtFacei = seedFaces[facei];
337 facei = mapFlag.find_next(facei);
343 Pout<<
"Advancing front stalled: searching for new " 344 <<
"target face" <<
endl;
351 tgtFacei = findTargetFace(srcFacei, visitedFaces);
358 facei = mapFlag.find_next(facei);
364 <<
"Unable to set target face for source face " << srcFacei
365 <<
" with centre: " << srcPatch.
faceCentres()[srcFacei]
375 const label srcFacei,
376 const label tgtFacei,
384 if (!isCandidate(srcFacei, tgtFacei))
389 const auto& srcPatch = this->srcPatch();
390 const auto& tgtPatch = this->tgtPatch();
396 faceAreaIntersect inter
416 scalar magN =
mag(
n);
418 const face& src = srcPatch[srcFacei];
419 const face& tgt = tgtPatch[tgtFacei];
421 if (magN > ROOTVSMALL)
423 inter.calc(src, tgt,
n/magN,
area, centroid);
428 <<
"Invalid normal for source face " << srcFacei
429 <<
" points " << UIndirectList<point>(srcPoints, src)
430 <<
" target face " << tgtFacei
431 <<
" points " << UIndirectList<point>(tgtPoints, tgt)
437 static OBJstream tris(
"intersectionTris.obj");
438 const auto& triPts = inter.triangles();
439 for (
const auto& tp : triPts)
441 tris.write(
triPointRef(tp[0], tp[1], tp[2]),
false);
447 writeIntersectionOBJ(
area, src, tgt, srcPoints, tgtPoints);
454 const label srcFacei,
455 const label tgtFacei,
456 const scalar threshold
460 if (!isCandidate(srcFacei, tgtFacei))
465 const auto& srcPatch = this->srcPatch();
466 const auto& tgtPatch = this->tgtPatch();
471 faceAreaIntersect inter
491 scalar magN =
mag(
n);
493 const face& src = srcPatch[srcFacei];
494 const face& tgt = tgtPatch[tgtFacei];
496 if (magN > ROOTVSMALL)
498 return inter.overlaps(src, tgt,
n/magN, threshold);
503 <<
"Invalid normal for source face " << srcFacei
504 <<
" points " << UIndirectList<point>(srcPoints, src)
505 <<
" target face " << tgtFacei
506 <<
" points " << UIndirectList<point>(tgtPoints, tgt)
516 List<DynamicList<label>>& srcAddr,
517 List<DynamicList<scalar>>& srcWght,
518 List<DynamicList<point>>& srcCtr,
519 List<DynamicList<label>>& tgtAddr,
520 List<DynamicList<scalar>>& tgtWght
523 addProfiling(ami,
"faceAreaWeightAMI::restartUncoveredSourceFace");
527 label nBelowMinWeight = 0;
528 const scalar minWeight = 0.95;
531 DynamicList<label> nbrFaces(10);
534 DynamicList<label> visitedFaces(10);
536 const auto& srcPatch = this->srcPatch();
540 const scalar
s =
sum(srcWght[srcFacei]);
541 const scalar t =
s/srcMagSf_[srcFacei];
547 const face&
f = srcPatch[srcFacei];
551 const label tgtFacei =
552 findTargetFace(srcFacei, srcAddr[srcFacei], fpi);
557 visitedFaces = srcAddr[srcFacei];
559 (void)processSourceFace
578 if (
debug && nBelowMinWeight)
581 <<
"Restarted search on " << nBelowMinWeight
582 <<
" faces since sum of weights < " << minWeight
592 const dictionary&
dict,
593 const bool reverseTarget
596 advancingFrontAMI(
dict, reverseTarget),
597 restartUncoveredSourceFace_
599 dict.getOrDefault(
"restartUncoveredSourceFace", true)
606 const bool requireMatch,
607 const bool reverseTarget,
608 const scalar lowWeightCorrection,
610 const bool restartUncoveredSourceFace
620 restartUncoveredSourceFace_(restartUncoveredSourceFace)
627 restartUncoveredSourceFace_(ami.restartUncoveredSourceFace_)
652 bool ok = initialiseWalk(srcFacei, tgtFacei);
654 srcCentroids_.setSize(srcAddress_.size());
656 const auto& src = this->srcPatch();
657 const auto& tgt = this->tgtPatch();
660 List<DynamicList<label>> srcAddr(src.size());
661 List<DynamicList<scalar>> srcWght(srcAddr.size());
662 List<DynamicList<point>> srcCtr(srcAddr.size());
663 List<DynamicList<label>> tgtAddr(tgt.size());
664 List<DynamicList<scalar>> tgtWght(tgtAddr.size());
679 if (
debug && !srcNonOverlap_.empty())
681 Pout<<
" AMI: " << srcNonOverlap_.size()
682 <<
" non-overlap faces identified" 687 if (restartUncoveredSourceFace_)
689 restartUncoveredSourceFace
703 srcAddress_[i].transfer(srcAddr[i]);
704 srcWeights_[i].transfer(srcWght[i]);
705 srcCentroids_[i].transfer(srcCtr[i]);
710 tgtAddress_[i].transfer(tgtAddr[i]);
711 tgtWeights_[i].transfer(tgtWght[i]);
720 globalIndex globalSrcFaces(srcPatch0.size());
721 globalIndex globalTgtFaces(tgtPatch0.size());
723 for (
labelList& addressing : srcAddress_)
725 for (label& addr : addressing)
727 addr = extendedTgtFaceIDs_[addr];
731 for (
labelList& addressing : tgtAddress_)
733 globalSrcFaces.inplaceToGlobal(addressing);
744 extendedTgtMapPtr_->constructMap(),
746 extendedTgtMapPtr_->subMap(),
750 ListOps::appendEqOp<label>(),
759 extendedTgtMapPtr_->constructMap(),
761 extendedTgtMapPtr_->subMap(),
765 ListOps::appendEqOp<scalar>(),
770 extendedTgtMapPtr_->reverseDistribute(tgtPatch0.size(), tgtMagSf_);
773 List<Map<label>> cMapSrc;
776 new mapDistribute(globalSrcFaces, tgtAddress_, cMapSrc)
779 List<Map<label>> cMapTgt;
782 new mapDistribute(globalTgtFaces, srcAddress_, cMapTgt)
787 normaliseWeights(requireMatch_,
true);
789 nonConformalCorrection();
801 if (restartUncoveredSourceFace_)
805 "restartUncoveredSourceFace",
806 restartUncoveredSourceFace_
faceAreaWeightAMI(const dictionary &dict, const bool reverseTarget=false)
Construct from dictionary.
List< scalar > scalarList
List of scalar.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
#define addProfiling(name, descr)
Define profiling trigger with specified name and description string.
static scalar & tolerance()
Fraction of local length scale to use as intersection tolerance.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
bool mustMatchFaces() const
Return true if requireMatch and lowWeightCorrectionin active.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
virtual bool setNextFaces(label &startSeedi, label &srcFacei, label &tgtFacei, const bitSet &mapFlag, labelList &seedFaces, const DynamicList< label > &visitedFaces, const bool errorOnNotFound=true) const
Set the source and target seed faces.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
static const List< T > & null()
Return a null List.
static bool cacheIntersections_
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
Face area weighted Arbitrary Mesh Interface (AMI) method.
virtual void write(Ostream &os) const
Write.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
addAliasToRunTimeSelectionTable(AMIInterpolation, faceAreaWeightAMI, dict, faceAreaWeightAMI, partialFaceAreaWeightAMI, 2012)
A list of faces which address into the list of points.
const labelListList & faceFaces() const
Return face-face addressing.
vectorField pointField
pointField is a vectorField.
bitSet & unset(const bitSet &other)
Unset (subtract) the bits specified in the other bitset, which is a set difference corresponds to the...
const wordList area
Standard area field types (scalar, vector, tensor, etc)
const Field< point_type > & faceCentres() const
Return face centres for patch.
virtual void calcInterArea(const label srcFacei, const label tgtFacei, scalar &area, vector ¢roid) const
Area of intersection between source and target faces.
virtual bool overlaps(const label srcFacei, const label tgtFacei, const scalar threshold) const
Return true if faces overlap.
const Field< point_type > & points() const noexcept
Return reference to global points.
void append(const T &val)
Copy append an element to the end of this list.
errorManip< error > abort(error &err)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
int debug
Static debugging option.
virtual void write(Ostream &os) const
Write.
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
virtual bool calculate(const primitivePatch &srcPatch, const primitivePatch &tgtPatch, const autoPtr< searchableSurface > &surfPtr=nullptr)
Update addressing, weights and (optional) centroids.
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
virtual bool calculate(const primitivePatch &srcPatch, const primitivePatch &tgtPatch, const autoPtr< searchableSurface > &surfPtr=nullptr)
Update addressing, weights and (optional) centroids.
virtual void restartUncoveredSourceFace(List< DynamicList< label >> &srcAddr, List< DynamicList< scalar >> &srcWght, List< DynamicList< point >> &srcCtr, List< DynamicList< label >> &tgtAddr, List< DynamicList< scalar >> &tgtWght)
Attempt to re-evaluate source faces that have not been included.
virtual bool processSourceFace(const label srcFacei, const label tgtStartFacei, DynamicList< label > &nbrFaces, DynamicList< label > &visitedFaces, List< DynamicList< label >> &srcAddr, List< DynamicList< scalar >> &srcWght, List< DynamicList< point >> &srcCtr, List< DynamicList< label >> &tgtAddr, List< DynamicList< scalar >> &tgtWght)
Determine overlap contributions for source face srcFacei.
virtual void calcAddressing(List< DynamicList< label >> &srcAddress, List< DynamicList< scalar >> &srcWeights, List< DynamicList< point >> &srcCentroids, List< DynamicList< label >> &tgtAddress, List< DynamicList< scalar >> &tgtWeights, label srcFacei, label tgtFacei)
Calculate addressing, weights and centroids using temporary storage.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
static void distribute(const Pstream::commsTypes commsType, const List< labelPair > &schedule, const label constructSize, const labelListList &subMap, const bool subHasFlip, const labelListList &constructMap, const bool constructHasFlip, List< T > &field, const NegateOp &negOp, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Distribute data with specified negate operator (for flips).
Base class for Arbitrary Mesh Interface (AMI) methods.
triangle< point, const point & > triPointRef
A triangle using referred points.
"nonBlocking" : (MPI_Isend, MPI_Irecv)
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
#define DebugVar(var)
Report a variable name and value.
List< label > labelList
A List of labels.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
labelList srcNonOverlap_
Labels of faces that are not overlapped by any target faces (should be empty for correct functioning ...
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
static constexpr const zero Zero
Global zero (0)