45 namespace cellCellStencils
69 static bool hasWarned =
false;
78 const fvPatch& fvp = pbm[patchi];
81 if (!fvPatch::constraintType(fvp.type()))
93 faceBb.
min() -= smallVec;
94 faceBb.
max() += smallVec;
114 const fvPatch& fvp = pbm[patchi];
117 if (isA<oversetFvPatch>(fvp))
125 patchCellTypes[cellMap[fc[i]]] = patchCellType::OVERSET;
129 faceBb.
min() -= smallVec;
130 faceBb.
max() += 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?"
150 voxelMeshSearch::overlaps
171 Pout<<
"Voxel mesh too coarse. Bounding box "
173 <<
" contains both non-overset and overset patches"
174 <<
". Refining voxel mesh to " << nDivs <<
endl;
179 voxelMeshSearch::fill
185 patchCellType::OVERSET
213 const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
217 const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
218 const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
220 if (srcPatchBb.
overlaps(tgtPatchBb))
225 forAll(tgtCellMap, tgtCelli)
227 label celli = tgtCellMap[tgtCelli];
229 cBb.
min() -= smallVec_;
230 cBb.
max() += smallVec_;
234 voxelMeshSearch::overlaps
244 allCellTypes[celli] = HOLE;
253 for (label procI = 0; procI < Pstream::nProcs(); procI++)
255 if (procI != Pstream::myProcNo())
257 const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
260 if (srcPatchBb.
overlaps(tgtPatchBb))
265 os << srcPatchBb << patchDivisions[srcI] << patchParts[srcI];
270 for (label procI = 0; procI < Pstream::nProcs(); procI++)
272 if (procI != Pstream::myProcNo())
275 const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
277 if (srcPatchBb.
overlaps(tgtPatchBb))
282 if (srcPatchBb != receivedBb)
286 <<
" srcPatchBb:" << srcPatchBb
287 <<
" receivedBb:" << receivedBb
294 forAll(tgtCellMap, tgtCelli)
296 label celli = tgtCellMap[tgtCelli];
298 cBb.
min() -= smallVec_;
299 cBb.
max() += smallVec_;
303 voxelMeshSearch::overlaps
313 allCellTypes[celli] = HOLE;
338 const fvMesh& srcMesh = meshParts_[srcI].subMesh();
339 const labelList& srcCellMap = meshParts_[srcI].cellMap();
341 const fvMesh& tgtMesh = meshParts_[tgtI].subMesh();
343 const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
347 forAll(tgtCellMap, tgtCelli)
350 if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
352 label celli = tgtCellMap[tgtCelli];
356 if (betterDonor(tgtI, allDonor[celli], srcI))
359 allStencil[celli][0] =
360 globalCells_.toGlobal(srcCellMap[srcCelli]);
361 allDonor[celli] = srcI;
376 for (label procI = 0; procI < Pstream::nProcs(); procI++)
378 if (procI != Pstream::myProcNo())
380 if (tgtBbs[procI].overlaps(srcBbs[Pstream::myProcNo()]))
382 tgtOverlapProcs.
append(procI);
384 if (srcBbs[procI].overlaps(tgtBbs[Pstream::myProcNo()]))
386 srcOverlapProcs.
append(procI);
395 forAll(srcOverlapProcs, i)
397 label procI = srcOverlapProcs[i];
398 tgtSendCells[procI].reserve(tgtMesh.
nCells()/srcOverlapProcs.size());
402 forAll(tgtCellMap, tgtCelli)
404 label celli = tgtCellMap[tgtCelli];
405 if (srcOverlapProcs.size())
408 subBb.
min() -= smallVec_;
409 subBb.
max() += smallVec_;
411 forAll(srcOverlapProcs, i)
413 label procI = srcOverlapProcs[i];
416 tgtSendCells[procI].
append(tgtCelli);
425 forAll(srcOverlapProcs, i)
427 label procI = srcOverlapProcs[i];
428 const labelList& cellIDs = tgtSendCells[procI];
431 os << UIndirectList<point>(tgtCc, cellIDs);
437 forAll(tgtOverlapProcs, i)
439 label procI = tgtOverlapProcs[i];
448 if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
450 donors[sampleI] = globalCells_.toGlobal(srcCellMap[srcCelli]);
460 forAll(srcOverlapProcs, i)
462 label procI = srcOverlapProcs[i];
463 const labelList& cellIDs = tgtSendCells[procI];
468 if (donors.size() != cellIDs.size())
476 label globalDonor = donors[donorI];
478 if (globalDonor != -1)
480 label celli = tgtCellMap[cellIDs[donorI]];
483 if (betterDonor(tgtI, allDonor[celli], srcI))
486 allStencil[celli][0] = globalDonor;
487 allDonor[celli] = srcI;
497 Foam::cellCellStencils::trackingInverseDistance::trackingInverseDistance
505 globalCells_(mesh_.nCells())
508 globalDonor_.setSize(mesh_.nCells());
523 meshParts_.setSize(
nZones);
534 (void)meshParts_[zonei].subMesh().nGeometricD();
541 <<
" mesh regions" <<
endl;
573 scalar layerRelax(dict_.getOrDefault(
"layerRelax", 1.0));
575 label
nZones = meshParts_.size();
580 pointField subPoints(mesh_.points(), meshParts_[zonei].pointMap());
582 fvMesh& subMesh = meshParts_[zonei].subMesh();
590 if (searchBoxDivisions.size() !=
nZones)
593 << searchBoxDivisions.size()
594 <<
" should equal the number of zones " <<
nZones
598 const boundBox& allBb(mesh_.bounds());
610 const fvMesh& subMesh = meshParts_[zonei].subMesh();
623 allBb.min() - 2*allBb.span(),
624 allBb.min() - allBb.span()
640 bbs[proci] = procBb[proci][zonei];
654 labelList allPatchTypes(mesh_.nCells(), OTHER);
658 if (dict_.readIfPresent(
"searchBox", globalPatchBb))
680 bool ok = markBoundaries
682 meshParts_[zonei].subMesh(),
686 patchDivisions[zonei],
689 meshParts_[zonei].cellMap(),
711 meshParts_[zonei].subMesh(),
713 patchDivisions[zonei],
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];
761 labelList allCellTypes(mesh_.nCells(), CALCULATED);
764 labelList allDonorID(mesh_.nCells(), -1);
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++)
841 createField(mesh_,
"allCellTypes", allCellTypes)
848 forAll(allPatchTypes, celli)
850 if (allCellTypes[celli] != HOLE)
852 switch (allPatchTypes[celli])
858 if (allStencil[celli].size())
860 allCellTypes[celli] = INTERPOLATED;
864 allCellTypes[celli] = HOLE;
876 createField(mesh_,
"allCellTypes_patch", allCellTypes)
887 label nCalculated = 0;
891 if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
893 if (allStencil[celli].size() == 0)
896 allCellTypes[celli] = HOLE;
897 allStencil[celli].
clear();
901 allCellTypes[celli] = INTERPOLATED;
909 Pout<<
"Detected " << nCalculated <<
" cells changing from hole"
910 <<
" to calculated. Changed to interpolated"
917 findHoles(globalCells_, mesh_,
zoneID, allStencil, allCellTypes);
924 createField(mesh_,
"allCellTypes_hole", allCellTypes)
931 walkFront(layerRelax, allStencil, allCellTypes, allWeight);
938 createField(mesh_,
"allCellTypes_front", allCellTypes)
946 cellTypes_.transfer(allCellTypes);
947 cellStencil_.setSize(mesh_.nCells());
948 cellInterpolationWeights_.setSize(mesh_.nCells());
952 if (cellTypes_[celli] == INTERPOLATED)
954 cellStencil_[celli].transfer(allStencil[celli]);
955 cellInterpolationWeights_[celli].setSize(1);
956 cellInterpolationWeights_[celli][0] = 1.0;
957 interpolationCells.
append(celli);
961 cellStencil_[celli].clear();
962 cellInterpolationWeights_[celli].clear();
965 interpolationCells_.transfer(interpolationCells);
968 cellInterpolationMap_.reset
977 cellInterpolationWeight_.
transfer(allWeight);
983 >(cellInterpolationWeight_.boundaryFieldRef(),
false);
989 mkDir(mesh_.time().timePath());
990 OBJstream str(mesh_.time().timePath()/
"injectionStencil.obj");
991 Pout<< typeName <<
" : dumping injectionStencil to "
994 cellInterpolationMap().distribute(cc);
996 forAll(cellStencil_, celli)
998 const labelList& slots = cellStencil_[celli];
1001 const point& accCc = mesh_.cellCentres()[celli];
1004 const point& donorCc = cc[slots[i]];
1015 createStencil(globalCells);
1022 cellInterpolationWeight_.instance() = mesh_.time().timeName();
1023 cellInterpolationWeight_.
write();
1031 mesh_.time().timeName(),
1039 zeroGradientFvPatchScalarField::typeName
1042 forAll(volTypes.internalField(), celli)
1044 volTypes[celli] = cellTypes_[celli];
1051 >(volTypes.boundaryFieldRef(),
false);
1055 mkDir(mesh_.time().timePath());
1056 OBJstream str(mesh_.time().timePath()/
"stencil.obj");
1057 Pout<< typeName <<
" : dumping to " << str.
name() <<
endl;
1059 cellInterpolationMap().distribute(cc);
1061 forAll(cellStencil_, celli)
1063 const labelList& slots = cellStencil_[celli];
1066 const point& accCc = mesh_.cellCentres()[celli];
1069 const point& donorCc = cc[slots[i]];
1084 forAll(interpolationCells_, i)
1086 label celli = interpolationCells_[i];
1087 const labelList& slots = cellStencil_[celli];
1089 bool hasLocal =
false;
1090 bool hasRemote =
false;
1094 if (slots[sloti] >= mesh_.nCells())
1124 Info<<
"Overset analysis : nCells : "
1127 <<
indent <<
"calculated : " << nCells[CALCULATED] <<
nl
1128 <<
indent <<
"interpolated : " << nCells[INTERPOLATED]
1129 <<
" (interpolated from local:" << nLocal
1130 <<
" mixed local/remote:" << nMixed
1131 <<
" remote:" << nRemote <<
")" <<
nl
1132 <<
indent <<
"hole : " << nCells[HOLE] <<
nl