46 namespace cellCellStencils
70 static bool hasWarned =
false;
79 const fvPatch& fvp = pbm[patchi];
94 faceBb.
grow(smallVec);
114 const fvPatch& fvp = pbm[patchi];
117 if (isA<oversetFvPatch>(fvp))
120 const polyPatch& pp = fvp.patch();
125 patchCellTypes[cellMap[fc[i]]] = patchCellType::OVERSET;
128 boundBox faceBb(pp.points(), pp[i],
false);
129 faceBb.grow(smallVec);
137 <<
" : face at " << faceCentres[i]
138 <<
" is not inside search bounding box of" 139 <<
" voxel mesh " << bb <<
endl 140 <<
" Is your 'searchBox' specification correct?" 171 Pout<<
"Voxel mesh too coarse. Bounding box " 173 <<
" contains both non-overset and overset patches" 174 <<
". Refining voxel mesh to " << nDivs <<
endl;
185 patchCellType::OVERSET
210 const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
217 if (srcPatchBb.
overlaps(tgtPatchBb))
222 forAll(tgtCellMap, tgtCelli)
224 label celli = tgtCellMap[tgtCelli];
240 allCellTypes[celli] = HOLE;
254 const treeBoundBox& tgtPatchBb = tgtPatchBbs[procI];
256 if (srcPatchBb.overlaps(tgtPatchBb))
260 UOPstream
os(procI, pBufs);
261 os << srcPatchBb << patchDivisions[srcI] << patchParts[srcI];
270 const treeBoundBox& srcPatchBb = srcPatchBbs[procI];
273 if (srcPatchBb.overlaps(tgtPatchBb))
275 UIPstream is(procI, pBufs);
277 treeBoundBox receivedBb(is);
278 if (srcPatchBb != receivedBb)
282 <<
" srcPatchBb:" << srcPatchBb
283 <<
" receivedBb:" << receivedBb
288 const PackedList<2> srcPatchTypes(is);
290 forAll(tgtCellMap, tgtCelli)
292 label celli = tgtCellMap[tgtCelli];
293 boundBox cBb(mesh_.cellBb(celli));
308 allCellTypes[celli] = HOLE;
319 PstreamBuffers& pBufs,
320 const List<treeBoundBoxList>&
meshBb,
321 const PtrList<voxelMeshSearch>& meshSearches,
333 const fvMesh& srcMesh = meshParts_[srcI].subMesh();
334 const labelList& srcCellMap = meshParts_[srcI].cellMap();
335 const fvMesh& tgtMesh = meshParts_[tgtI].subMesh();
336 const pointField& tgtCc = tgtMesh.cellCentres();
337 const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
344 forAll(tgtCellMap, tgtCelli)
346 const label srcCelli = tgtToSrcAddr[tgtCelli];
351 && allCellTypes[tgtCellMap[tgtCelli]] != HOLE
354 label celli = tgtCellMap[tgtCelli];
358 if (betterDonor(tgtI, allDonor[celli], srcI))
360 allStencil[celli].setSize(1);
361 allStencil[celli][0] =
362 globalCells_.toGlobal(srcCellMap[srcCelli]);
363 allDonor[celli] = srcI;
384 tgtOverlapProcs.append(procI);
388 srcOverlapProcs.append(procI);
397 forAll(srcOverlapProcs, i)
399 label procI = srcOverlapProcs[i];
400 tgtSendCells[procI].reserve(tgtMesh.nCells()/srcOverlapProcs.size());
404 forAll(tgtCellMap, tgtCelli)
406 label celli = tgtCellMap[tgtCelli];
407 if (srcOverlapProcs.size())
409 treeBoundBox subBb(mesh_.cellBb(celli));
410 subBb.grow(smallVec_);
412 forAll(srcOverlapProcs, i)
414 label procI = srcOverlapProcs[i];
415 if (subBb.overlaps(srcBbs[procI]))
417 tgtSendCells[procI].append(tgtCelli);
426 forAll(srcOverlapProcs, i)
428 label procI = srcOverlapProcs[i];
431 UOPstream
os(procI, pBufs);
432 os << UIndirectList<point>(tgtCc,
cellIDs);
434 pBufs.finishedSends();
437 (void)srcMesh.tetBasePtIs();
438 forAll(tgtOverlapProcs, i)
440 label procI = tgtOverlapProcs[i];
442 UIPstream is(procI, pBufs);
452 donors[sampleI] = globalCells_.toGlobal(srcCellMap[srcCelli]);
457 UOPstream
os(procI, pBufs);
460 pBufs.finishedSends();
462 forAll(srcOverlapProcs, i)
464 label procI = srcOverlapProcs[i];
467 UIPstream is(procI, pBufs);
470 if (donors.size() !=
cellIDs.size())
478 label globalDonor = donors[donorI];
480 if (globalDonor != -1)
482 label celli = tgtCellMap[
cellIDs[donorI]];
485 if (betterDonor(tgtI, allDonor[celli], srcI))
487 allStencil[celli].setSize(1);
488 allStencil[celli][0] = globalDonor;
489 allDonor[celli] = srcI;
499 Foam::cellCellStencils::trackingInverseDistance::trackingInverseDistance
507 globalCells_(mesh_.nCells())
520 nCellsPerZone[
zoneID[celli]]++;
535 (void)
meshParts_[zonei].subMesh().nGeometricD();
541 Info<< typeName <<
" : detected " << nZones
542 <<
" mesh regions" <<
endl;
544 forAll(nCellsPerZone, zonei)
547 <<
" nCells:" << nCellsPerZone[zonei]
574 scalar layerRelax(dict_.getOrDefault(
"layerRelax", 1.0));
576 label nZones = meshParts_.size();
581 pointField subPoints(mesh_.points(), meshParts_[zonei].pointMap());
583 fvMesh& subMesh = meshParts_[zonei].subMesh();
591 if (searchBoxDivisions.size() != nZones)
594 << searchBoxDivisions.size()
595 <<
" should equal the number of zones " << nZones
599 const boundBox& allBb = mesh_.bounds();
601 List<treeBoundBoxList>
meshBb(nZones);
611 const fvMesh& subMesh = meshParts_[zonei].subMesh();
613 if (subMesh.nPoints())
616 treeBoundBox(subMesh.points());
624 allBb.min() - 2*allBb.span(),
625 allBb.min() - allBb.span()
640 bbs[proci] = procBb[proci][zonei];
651 List<treeBoundBoxList> patchBb(nZones);
652 List<labelVector> patchDivisions(searchBoxDivisions);
653 PtrList<PackedList<2>> patchParts(nZones);
654 labelList allPatchTypes(mesh_.nCells(), OTHER);
657 treeBoundBox globalPatchBb;
658 if (dict_.readIfPresent(
"searchBox", globalPatchBb))
677 new PackedList<2>(
cmptProduct(patchDivisions[zonei]))
680 bool ok = markBoundaries
682 meshParts_[zonei].subMesh(),
686 patchDivisions[zonei],
689 meshParts_[zonei].cellMap(),
703 PtrList<voxelMeshSearch> meshSearches(meshParts_.size());
711 meshParts_[zonei].subMesh(),
713 patchDivisions[zonei],
721 PtrList<fvMesh> voxelMeshes;
724 voxelMeshes.setSize(meshSearches.size());
725 forAll(meshSearches, zonei)
730 mesh_.time().timeName(),
735 Pout<<
"Constructing voxel mesh " << meshIO.objectPath() <<
endl;
736 voxelMeshes.set(zonei, meshSearches[zonei].makeMesh(meshIO));
746 const fvMesh& vm = voxelMeshes[zonei];
747 tmp<volScalarField> tfld
761 labelList allCellTypes(mesh_.nCells(), CALCULATED);
764 labelList allDonorID(mesh_.nCells(), -1);
766 const globalIndex globalCells(mesh_.nCells());
772 for (label srci = 0; srci < meshParts_.size()-1; srci++)
774 for (label tgti = srci+1; tgti < meshParts_.size(); tgti++)
803 for (label srci = 0; srci < meshParts_.size()-1; srci++)
805 for (label tgti = srci+1; tgti < meshParts_.size(); tgti++)
837 if ((
debug & 2) && (mesh_.time().outputTime()))
839 tmp<volScalarField> tfld
841 createField(mesh_,
"allCellTypes", allCellTypes)
845 tmp<volScalarField> tallDonorID
847 createField(mesh_,
"allDonorID", allDonorID)
849 tallDonorID().write();
854 forAll(allPatchTypes, celli)
856 if (allCellTypes[celli] != HOLE)
858 switch (allPatchTypes[celli])
864 if (allStencil[celli].size())
866 allCellTypes[celli] = INTERPOLATED;
870 allCellTypes[celli] = HOLE;
878 if ((
debug & 2) && (mesh_.time().outputTime()))
880 tmp<volScalarField> tfld
882 createField(mesh_,
"allCellTypes_patch", allCellTypes)
886 tmp<volScalarField> tfldOld
888 createField(mesh_,
"allCellTypes_old", cellTypes_)
894 findHoles(globalCells_, mesh_, zoneID, allStencil, allCellTypes);
897 if ((
debug & 2) && (mesh_.time().outputTime()))
899 tmp<volScalarField> tfld
901 createField(mesh_,
"allCellTypes_hole", allCellTypes)
908 stencilSize[celli] = allStencil[celli].size();
910 tmp<volScalarField> tfldsten
912 createField(mesh_,
"allStencil_hole", stencilSize)
918 forAll(allCellTypes, celli)
920 if (allCellTypes[celli] == HOLE && allStencil[celli].size())
922 allStencil[celli].clear();
930 List<Map<label>> compactStencilMap;
931 mapDistribute map(globalCells, compactStencil, compactStencilMap);
934 map.distribute(compactCellVol);
946 dict_.getOrDefault(
"holeLayers", 1),
947 dict_.getOrDefault(
"useLayer", -1)
950 if ((
debug & 2) && (mesh_.time().outputTime()))
952 tmp<volScalarField> tfld
954 createField(mesh_,
"allCellTypes_front", allCellTypes)
964 label nCalculated = 0;
968 if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
970 if (allStencil[celli].size() == 0)
973 allCellTypes[celli] = HOLE;
974 allStencil[celli].clear();
978 allCellTypes[celli] = INTERPOLATED;
986 Pout<<
"Detected " << nCalculated <<
" cells changing from hole" 987 <<
" to calculated. Changed to interpolated" 995 cellTypes_.transfer(allCellTypes);
996 cellStencil_.setSize(mesh_.nCells());
997 cellInterpolationWeights_.setSize(mesh_.nCells());
998 DynamicList<label> interpolationCells;
1001 map.distribute(compactCellTypes);
1003 label nHoleDonors = 0;
1004 forAll(cellTypes_, celli)
1006 if (cellTypes_[celli] == INTERPOLATED)
1008 const labelList& slots = compactStencil[celli];
1013 compactCellTypes[slots[0]] == HOLE
1016 !allowInterpolatedDonors_
1017 && compactCellTypes[slots[0]] == INTERPOLATED
1021 cellTypes_[celli] = POROUS;
1022 cellStencil_[celli].
clear();
1023 cellInterpolationWeights_[celli].clear();
1028 cellStencil_[celli].transfer(allStencil[celli]);
1029 cellInterpolationWeights_[celli].setSize(1);
1030 cellInterpolationWeights_[celli][0] = 1.0;
1031 interpolationCells.append(celli);
1037 cellStencil_[celli].clear();
1038 cellInterpolationWeights_[celli].clear();
1042 reduce(nHoleDonors, sumOp<label>());
1046 interpolationCells_.transfer(interpolationCells);
1048 List<Map<label>> compactMap;
1049 cellInterpolationMap_.reset
1058 cellInterpolationWeight_.transfer(allWeight);
1062 oversetFvPatchField<scalar>
1063 >(cellInterpolationWeight_.boundaryFieldRef(),
false);
1066 if ((
debug & 2) && (mesh_.time().outputTime()))
1072 mkDir(mesh_.time().timePath());
1073 OBJstream str(mesh_.time().timePath()/
"injectionStencil.obj");
1074 Pout<< typeName <<
" : dumping injectionStencil to " 1077 cellInterpolationMap().distribute(cc);
1079 forAll(cellStencil_, celli)
1081 const labelList& slots = cellStencil_[celli];
1084 const point& accCc = mesh_.cellCentres()[celli];
1087 const point& donorCc = cc[slots[i]];
1088 str.writeLine(accCc, 0.1*accCc+0.9*donorCc);
1098 createStencil(globalCells, allowHoleDonors_);
1103 if (allowHoleDonors_)
1105 holeExtrapolationStencil(globalCells_);
1109 if ((
debug & 2) && (mesh_.time().outputTime()))
1112 cellInterpolationWeight_.instance() = mesh_.time().timeName();
1113 cellInterpolationWeight_.write();
1116 mkDir(mesh_.time().timePath());
1117 OBJstream str(mesh_.time().timePath()/
"stencil.obj");
1118 Pout<< typeName <<
" : dumping to " << str.
name() <<
endl;
1120 cellInterpolationMap().distribute(cc);
1122 forAll(cellStencil_, celli)
1124 const labelList& slots = cellStencil_[celli];
1127 const point& accCc = mesh_.cellCentres()[celli];
1130 const point& donorCc = cc[slots[i]];
1131 str.writeLine(accCc, 0.1*accCc+0.9*donorCc);
1145 forAll(interpolationCells_, i)
1147 label celli = interpolationCells_[i];
1148 const labelList& slots = cellStencil_[celli];
1150 bool hasLocal =
false;
1151 bool hasRemote =
false;
1155 if (slots[sloti] >= mesh_.nCells())
1182 Info<<
"Overset analysis : nCells : " 1185 <<
indent <<
"calculated : " << nCells[CALCULATED] <<
nl 1186 <<
indent <<
"interpolated : " << nCells[INTERPOLATED]
1187 <<
" (from local:" <<
returnReduce(nLocal, sumOp<label>())
1188 <<
" mixed local/remote:" <<
returnReduce(nMixed, sumOp<label>())
1189 <<
" remote:" <<
returnReduce(nRemote, sumOp<label>()) <<
")" <<
nl 1190 <<
indent <<
"hole : " << nCells[HOLE] <<
nl virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
A List of labelList.
Inverse-distance-weighted interpolation stencil.
Ostream & indent(Ostream &os)
Indent stream.
Cmpt cmptProduct(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
errorManipArg< error, int > exit(error &err, const int errNo=1)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
bool contains(const point &pt) const
Contains point? (inside or on edge)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
bool overlaps(const boundBox &bb) const
Overlaps with other bounding box, sphere etc?
labelList globalDonor_
Current (global) donor cell.
static rangeType allProcs(const label communicator=worldComm)
Range of process indices for all processes.
constexpr char nl
The newline '\n' character (0x0a)
wordList patchTypes(nPatches)
static void calculate(const polyMesh &src, const polyMesh &tgt, labelList &srcToTgtAddr)
Calculate addressing.
scalarField samples(nIntervals, Zero)
Ostream & endl(Ostream &os)
Add newline and flush stream.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
virtual void movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
defineTypeNameAndDebug(cellVolumeWeight, 0)
PtrList< fvMeshSubset > meshParts_
Subset according to zone.
A bounding box defined in terms of min/max extrema points.
List< treeBoundBox > treeBoundBoxList
A List of treeBoundBox.
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
List< point > pointList
A List of points.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
static void allGatherList(List< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses linear/tree communication. ...
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator) is 1 for serial run.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Calculation of interpolation stencils.
void clear()
Clear the list, i.e. set size to zero.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual const labelUList & faceCells() const
Return faceCells.
void finishedSends(const bool wait=true)
Mark sends as done.
List< scalar > scalarList
A List of scalars.
static bool overlaps(const boundBox &bb, const labelVector &nDivs, const boundBox &subBb, const Container &elems, const Type val, const bool isNot=false)
Check if any voxel inside bounding box is set to val or.
void markPatchesAsHoles(PstreamBuffers &pBufs, const List< treeBoundBoxList > &patchBb, const List< labelVector > &patchDivisions, const PtrList< PackedList< 2 >> &patchParts, const label srcI, const label tgtI, labelList &allCellTypes) const
Mark all cells overlapping (a voxel covered by) a src patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
errorManip< error > abort(error &err)
const fvMesh & mesh_
Reference to the mesh.
#define DebugInfo
Report an information message using Foam::Info.
static void correctBoundaryConditions(typename GeoField::Boundary &bfld, const bool typeOnly)
Correct boundary conditions of certain type (typeOnly = true)
virtual bool update()
Update stencils. Return false if nothing changed.
addToRunTimeSelectionTable(cellCellStencil, cellVolumeWeight, mesh)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
Type gMax(const FieldField< Field, Type > &f)
OBJstream os(runTime.globalPath()/outputName)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
static bool constraintType(const word &patchType)
Return true if the given type is a constraint type.
void clear()
Clear individual buffers and reset states.
const labelIOList & zoneID() const
Helper: get reference to registered zoneID. Loads volScalarField.
void markDonors(PstreamBuffers &pBufs, const List< treeBoundBoxList > &meshBb, const PtrList< voxelMeshSearch > &meshSearches, const labelList &allCellTypes, const label srcI, const label tgtI, labelListList &allStencil, labelList &allDonor) const
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
label nCells() const noexcept
Number of mesh cells.
const Time & time() const
Return Time associated with the objectRegistry.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
Mesh data needed to do the Finite Volume discretisation.
static void fill(Container &elems, const boundBox &bb, const labelVector &nDivs, const boundBox &subBb, const Type val)
Fill voxels indicated by bounding box.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
virtual ~trackingInverseDistance()
Destructor.
"nonBlocking" : (MPI_Isend, MPI_Irecv)
Standard boundBox with extra functionality for use in octree.
messageStream Info
Information stream (stdout output on master, null elsewhere)
SubField< vector > subField
Declare type of subField.
Vector< label > labelVector
Vector of labels.
List< label > labelList
A List of labels.
A patch is a list of labels that address the faces in the global face list.
const polyPatch & patch() const noexcept
Return the polyPatch.
virtual bool write(const bool valid=true) const
Write using setting from DB.
bool overlaps(const boundBox &bb) const
Overlaps/touches boundingBox?
Ostream & incrIndent(Ostream &os)
Increment the indent level.
void grow(const scalar delta)
Expand box by adjusting min/max by specified amount in each dimension.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< treeBoundBox > meshBb(1, treeBoundBox(coarseMesh.points()).extend(rndGen, 1e-3))
static bool markBoundaries(const fvMesh &mesh, const vector &smallVec, const boundBox &bb, labelVector &nDivs, PackedList< 2 > &patchTypes, const labelList &cellMap, labelList &patchCellTypes)
Mark voxels of patchTypes with type of patch face.
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
static constexpr const zero Zero
Global zero (0)