50 namespace cellCellStencils
65 return (ijk[0]*nDivs[1] + ijk[1])*nDivs[2] + ijk[2];
75 label ij = boxI/nDivs[2];
76 label
k = boxI-ij*nDivs[2];
77 label i = ij/nDivs[1];
78 label j = ij-i*nDivs[1];
96 floor(relPt[0]/d[0]*nDivs[0]),
97 floor(relPt[1]/d[1]*nDivs[1]),
98 floor(relPt[2]/d[2]*nDivs[2])
114 const vector sz(d[0]/nDivs[0], d[1]/nDivs[1], d[2]/nDivs[2]);
116 return bb.
min()+0.5*sz+
vector(sz[0]*ids[0], sz[1]*ids[1], sz[2]*ids[2]);
126 const unsigned int val
132 for (
direction cmpt = 0; cmpt < 3; cmpt++)
134 if (maxIds[cmpt] < 0 || minIds[cmpt] > nDivs[cmpt])
141 minIds =
max(labelVector::zero, minIds);
142 maxIds =
min(maxIndex, maxIds);
144 for (label i = minIds[0]; i <= maxIds[0]; i++)
146 for (label j = minIds[1]; j <= maxIds[1]; j++)
148 for (label
k = minIds[2];
k <= maxIds[2];
k++)
180 const fvPatch& fvp = pbm[patchI];
183 if (!fvPatch::constraintType(fvp.type()))
195 faceBb.
min() -= smallVec;
196 faceBb.
max() += smallVec;
209 const fvPatch& fvp = pbm[patchI];
212 if (isA<oversetFvPatch>(fvp))
219 patchCellTypes[cellMap[fc[i]]] = patchCellType::OVERSET;
223 faceBb.
min() -= smallVec;
224 faceBb.
max() += smallVec;
228 fill(
patchTypes, bb, nDivs, faceBb, patchCellType::OVERSET);
248 vector(GREAT, GREAT, GREAT),
249 vector(-GREAT, -GREAT, -GREAT)
256 const face&
f = faces[cFaces[cFacei]];
276 const unsigned int val
284 for (
direction cmpt = 0; cmpt < 3; cmpt++)
286 if (maxIds[cmpt] < 0 || minIds[cmpt] > nDivs[cmpt])
293 minIds =
max(labelVector::zero, minIds);
294 maxIds =
min(maxIndex, maxIds);
296 for (label i = minIds[0]; i <= maxIds[0]; i++)
298 for (label j = minIds[1]; j <= maxIds[1]; j++)
300 for (label
k = minIds[2];
k <= maxIds[2];
k++)
335 const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
336 const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
338 if (srcPatchBb.
overlaps(tgtPatchBb))
341 const labelVector& zoneDivs = patchDivisions[srcI];
343 forAll(tgtCellMap, tgtCelli)
345 label celli = tgtCellMap[tgtCelli];
347 cBb.
min() -= smallVec_;
348 cBb.
max() += smallVec_;
362 allCellTypes[celli] = HOLE;
371 for (
const int procI : Pstream::allProcs())
373 if (procI != Pstream::myProcNo())
375 const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
378 if (srcPatchBb.
overlaps(tgtPatchBb))
383 os << srcPatchBb << patchDivisions[srcI] << patchParts[srcI];
388 for (
const int procI : Pstream::allProcs())
390 if (procI != Pstream::myProcNo())
394 const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
396 if (srcPatchBb.
overlaps(tgtPatchBb))
401 if (srcPatchBb != receivedBb)
405 <<
" srcPatchBb:" << srcPatchBb
406 <<
" receivedBb:" << receivedBb
413 forAll(tgtCellMap, tgtCelli)
415 label celli = tgtCellMap[tgtCelli];
417 cBb.
min() -= smallVec_;
418 cBb.
max() += smallVec_;
431 allCellTypes[celli] = HOLE;
442 const label destMesh,
443 const label currentDonorMesh,
444 const label newDonorMesh
457 if (currentDonorMesh == -1)
463 const label currentDist =
mag(currentDonorMesh-destMesh);
464 const label newDist =
mag(newDonorMesh-destMesh);
466 if (newDist < currentDist)
470 else if (newDist == currentDist && newDonorMesh > currentDonorMesh)
509 waveMethod::calculate(tgtMesh, srcMesh, tgtToSrcAddr);
510 forAll(tgtCellMap, tgtCelli)
512 label srcCelli = tgtToSrcAddr[tgtCelli];
513 if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
515 label celli = tgtCellMap[tgtCelli];
519 if (betterDonor(tgtI, allDonor[celli], srcI))
522 globalCells.
toGlobal(srcCellMap[srcCelli]);
524 allStencil[celli][0] = globalDonor;
525 allDonor[celli] = srcI;
540 for (
const int procI : Pstream::allProcs())
542 if (procI != Pstream::myProcNo())
544 if (tgtBbs[procI].overlaps(srcBbs[Pstream::myProcNo()]))
546 tgtOverlapProcs.
append(procI);
548 if (srcBbs[procI].overlaps(tgtBbs[Pstream::myProcNo()]))
550 srcOverlapProcs.
append(procI);
559 forAll(srcOverlapProcs, i)
561 label procI = srcOverlapProcs[i];
562 tgtSendCells[procI].reserve(tgtMesh.
nCells()/srcOverlapProcs.size());
566 forAll(tgtCellMap, tgtCelli)
568 label celli = tgtCellMap[tgtCelli];
569 if (srcOverlapProcs.size())
572 subBb.
min() -= smallVec_;
573 subBb.
max() += smallVec_;
575 forAll(srcOverlapProcs, i)
577 label procI = srcOverlapProcs[i];
580 tgtSendCells[procI].
append(tgtCelli);
589 forAll(srcOverlapProcs, i)
591 label procI = srcOverlapProcs[i];
592 const labelList& cellIDs = tgtSendCells[procI];
595 os << UIndirectList<point>(tgtCc, cellIDs);
601 forAll(tgtOverlapProcs, i)
603 label procI = tgtOverlapProcs[i];
613 if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
615 donors[sampleI] = globalCells.
toGlobal(srcCellMap[srcCelli]);
625 forAll(srcOverlapProcs, i)
627 label procI = srcOverlapProcs[i];
628 const labelList& cellIDs = tgtSendCells[procI];
633 if (donors.size() != cellIDs.size())
641 label globalDonor = donors[donorI];
643 if (globalDonor != -1)
645 label celli = tgtCellMap[cellIDs[donorI]];
649 if (betterDonor(tgtI, allDonor[celli], srcI))
652 allStencil[celli][0] = globalDonor;
653 allDonor[celli] = srcI;
976 for (label faceI = 0; faceI <
mesh.nInternalFaces(); faceI++)
982 (ownType == HOLE && neiType != HOLE)
983 || (ownType != HOLE && neiType == HOLE)
986 isBlockedFace[faceI] =
true;
995 syncTools::swapBoundaryCellList(
mesh,
cellTypes, nbrCellTypes);
997 for (label faceI =
mesh.nInternalFaces(); faceI <
mesh.nFaces(); faceI++)
1000 label neiType = nbrCellTypes[faceI-
mesh.nInternalFaces()];
1004 (ownType == HOLE && neiType != HOLE)
1005 || (ownType != HOLE && neiType == HOLE)
1008 isBlockedFace[faceI] =
true;
1014 << nBlocked <<
endl;
1018 const label nRegions = cellRegion.
nRegions();
1046 << nRegions <<
endl;
1067 for (label faceI = 0; faceI <
mesh.nInternalFaces(); faceI++)
1069 if (isBlockedFace[faceI])
1071 label ownRegion = cellRegion[own[faceI]];
1075 if (regionType[ownRegion] == 0)
1077 regionType[ownRegion] = 1;
1081 label neiRegion = cellRegion[nei[faceI]];
1085 if (regionType[neiRegion] == 0)
1087 regionType[neiRegion] = 1;
1094 label faceI =
mesh.nInternalFaces();
1095 faceI <
mesh.nFaces();
1099 if (isBlockedFace[faceI])
1101 label ownRegion = cellRegion[own[faceI]];
1103 if (regionType[ownRegion] == 0)
1105 regionType[ownRegion] = 1;
1115 const fvPatch& fvp = pbm[patchI];
1117 if (isA<oversetFvPatch>(fvp))
1119 else if (!fvPatch::constraintType(fvp.type()))
1124 label regionI = cellRegion[fc[i]];
1126 if (
cellTypes[fc[i]] != HOLE && regionType[regionI] != 2)
1128 regionType[regionI] = 2;
1155 Pstream::listCombineScatter(regionType);
1170 const fvPatch& fvp = pbm[patchI];
1172 if (isA<oversetFvPatch>(fvp))
1177 label cellI = fc[i];
1178 label regionI = cellRegion[cellI];
1180 if (regionType[regionI] != 2)
1182 const labelList& slots = compactStencil[cellI];
1185 label otherType = cellRegionType[slots[i]];
1192 regionType[regionI] = 2;
1204 << nChanged <<
endl;
1214 forAll(cellRegion, cellI)
1216 label
type = regionType[cellRegion[cellI]];
1229 const scalar wantedFraction,
1234 const cell& cFaces = mesh_.cells()[cellI];
1237 label nbrFacei = cFaces[i];
1238 if (fraction[nbrFacei] < wantedFraction)
1240 fraction[nbrFacei] = wantedFraction;
1241 isFront.
set(nbrFacei);
1249 const scalar layerRelax,
1256 bitSet isFront(mesh_.nFaces());
1265 if (isA<oversetFvPatch>(fvm[patchI]))
1267 const labelList& fc = fvm[patchI].faceCells();
1270 label cellI = fc[i];
1271 if (allCellTypes[cellI] == INTERPOLATED)
1276 isFront.
set(fvm[patchI].start()+i);
1285 const labelList& own = mesh_.faceOwner();
1286 const labelList& nei = mesh_.faceNeighbour();
1288 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1290 label ownType = allCellTypes[own[faceI]];
1291 label neiType = allCellTypes[nei[faceI]];
1294 (ownType == HOLE && neiType != HOLE)
1295 || (ownType != HOLE && neiType == HOLE)
1305 syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
1309 label faceI = mesh_.nInternalFaces();
1310 faceI < mesh_.nFaces();
1314 label ownType = allCellTypes[own[faceI]];
1315 label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];
1319 (ownType == HOLE && neiType != HOLE)
1320 || (ownType != HOLE && neiType == HOLE)
1336 if (isFront.
test(faceI))
1338 fraction[faceI] = 1.0;
1346 bitSet newIsFront(mesh_.nFaces());
1350 if (isFront.
test(faceI))
1352 label own = mesh_.faceOwner()[faceI];
1353 if (allCellTypes[own] != HOLE)
1355 if (allWeight[own] < fraction[faceI])
1359 if (allStencil[own].size())
1361 allWeight[own] = fraction[faceI];
1362 allCellTypes[own] = INTERPOLATED;
1368 fraction[faceI]-layerRelax,
1375 allWeight[own] = 0.0;
1376 allCellTypes[own] = HOLE;
1388 if (mesh_.isInternalFace(faceI))
1390 label nei = mesh_.faceNeighbour()[faceI];
1391 if (allCellTypes[nei] != HOLE)
1393 if (allWeight[nei] < fraction[faceI])
1395 if (allStencil[nei].size())
1397 allWeight[nei] = fraction[faceI];
1398 allCellTypes[nei] = INTERPOLATED;
1402 fraction[faceI]-layerRelax,
1409 allWeight[nei] = 0.0;
1410 allCellTypes[nei] = HOLE;
1429 fraction.transfer(newFraction);
1443 weights.
setSize(donorCcs.size());
1489 const vector greatPoint(GREAT, GREAT, GREAT);
1491 boolList isValidDonor(mesh_.nCells(),
true);
1492 forAll(cellTypes_, celli)
1494 if (cellTypes_[celli] == HOLE)
1496 isValidDonor[celli] =
false;
1502 bitSet doneAcceptor(interpolationCells_.size());
1511 forAll(interpolationCells_, i)
1513 if (!doneAcceptor[i])
1515 label cellI = interpolationCells_[i];
1516 const point& cc = mesh_.cellCentres()[cellI];
1517 const labelList& slots = cellStencil_[cellI];
1519 if (slots.size() != 1)
1527 label elemI = slots[slotI];
1545 mapDistributeBase::distribute<point, minMagSqrEqOp<point>,
flipOp>
1547 Pstream::commsTypes::nonBlocking,
1550 cellInterpolationMap().constructMap(),
1552 cellInterpolationMap().subMap(),
1558 UPstream::msgType(),
1559 cellInterpolationMap().comm()
1568 if (
samples[cellI] != greatPoint)
1570 donorCells.append(cellI);
1592 label cellI = donorCells[i];
1593 const pointList& donorCentres = donorCellCentres[cellI];
1605 cellInterpolationMap().distribute(donorCellCells);
1606 cellInterpolationMap().distribute(donorWeights);
1607 cellInterpolationMap().distribute(
samples);
1610 forAll(interpolationCells_, i)
1612 if (!doneAcceptor[i])
1614 label cellI = interpolationCells_[i];
1615 const labelList& slots = cellStencil_[cellI];
1617 if (slots.size() != 1)
1623 label slotI = slots[0];
1626 if (
samples[slotI] == mesh_.cellCentres()[cellI])
1628 cellStencil_[cellI].
transfer(donorCellCells[slotI]);
1629 cellInterpolationWeights_[cellI].transfer
1635 doneAcceptor.set(i);
1643 cellInterpolationMap_.reset
1657 Foam::cellCellStencils::inverseDistance::inverseDistance
1668 interpolationCells_(0),
1669 cellInterpolationMap_(),
1671 cellInterpolationWeights_(0),
1672 cellInterpolationWeight_
1676 "cellInterpolationWeight",
1677 mesh_.facesInstance(),
1685 zeroGradientFvPatchScalarField::typeName
1689 nonInterpolatedFields_.insert(
"cellInterpolationWeight");
1690 nonInterpolatedFields_.insert(
"cellTypes");
1691 nonInterpolatedFields_.insert(
"maxMagWeight");
1694 nonInterpolatedFields_.insert(
"cellDisplacement");
1695 nonInterpolatedFields_.insert(
"grad(cellDisplacement)");
1696 const word w(
"snGradCorr(cellDisplacement)");
1697 const word d(
"((viscosity*faceDiffusivity)*magSf)");
1698 nonInterpolatedFields_.insert(
"surfaceIntegrate(("+d+
"*"+w+
"))");
1707 mesh_.time().timeName(),
1709 IOobject::READ_IF_PRESENT,
1717 Pout<<
"Reading cellTypes from time " << mesh_.time().timeName()
1722 forAll(volCellTypes, celli)
1725 cellTypes_[celli] = volCellTypes[celli];
1746 scalar layerRelax(dict_.getOrDefault(
"layerRelax", 1.0));
1748 scalar tol = dict_.getOrDefault(
"tolerance", 1
e-10);
1749 smallVec_ = mesh_.bounds().span()*tol;
1763 const boundBox& allBb(mesh_.bounds());
1799 allBb.min() - 2*allBb.span(),
1800 allBb.min() - allBb.span()
1816 bbs[procI] = procBb[procI][zoneI];
1828 labelList allPatchTypes(mesh_.nCells(), OTHER);
1832 if (dict_.readIfPresent(
"searchBox", globalPatchBb))
1846 if (dict_.readIfPresent(
"searchBoxDivisions", globalDivs))
1848 patchDivisions = globalDivs;
1854 if (mesh_.nGeometricD() == 1)
1856 nDivs = mesh_.nCells();
1858 else if (mesh_.nGeometricD() == 2)
1860 nDivs = label(
Foam::sqrt(scalar(mesh_.nCells())));
1864 nDivs = label(
Foam::cbrt(scalar(mesh_.nCells())));
1879 forAll(patchParts, zoneI)
1886 patchDivisions[zoneI][0]
1887 *patchDivisions[zoneI][1]
1888 *patchDivisions[zoneI][2]
1897 patchDivisions[zoneI],
1909 <<
" mesh regions" <<
endl;
1915 <<
" voxels:" << patchDivisions[zoneI]
1925 labelList allCellTypes(mesh_.nCells(), CALCULATED);
1928 labelList allDonorID(mesh_.nCells(), -1);
1935 for (label srcI = 0; srcI <
meshParts.size()-1; srcI++)
1937 for (label tgtI = srcI+1; tgtI <
meshParts.size(); tgtI++)
1971 for (label srcI = 0; srcI <
meshParts.size()-1; srcI++)
1973 for (label tgtI = srcI+1; tgtI <
meshParts.size(); tgtI++)
2008 createField(mesh_,
"allCellTypes", allCellTypes)
2014 forAll(allPatchTypes, cellI)
2016 if (allCellTypes[cellI] != HOLE)
2018 switch (allPatchTypes[cellI])
2024 if (allStencil[cellI].size())
2026 allCellTypes[cellI] = INTERPOLATED;
2040 allCellTypes[cellI] = HOLE;
2051 createField(mesh_,
"allCellTypes_patch", allCellTypes)
2062 label nCalculated = 0;
2064 forAll(cellTypes_, celli)
2066 if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
2068 if (allStencil[celli].size() == 0)
2071 allCellTypes[celli] = HOLE;
2072 allStencil[celli].
clear();
2076 allCellTypes[celli] = INTERPOLATED;
2084 Pout<<
"Detected " << nCalculated <<
" cells changing from hole"
2085 <<
" to calculated. Changed to interpolated"
2092 findHoles(globalCells, mesh_,
zoneID, allStencil, allCellTypes);
2098 createField(mesh_,
"allCellTypes_hole", allCellTypes)
2105 forAll(allStencil, celli)
2107 stencilSize[celli] = allStencil[celli].size();
2111 createField(mesh_,
"allStencil_hole", stencilSize)
2119 walkFront(layerRelax, allStencil, allCellTypes, allWeight);
2125 createField(mesh_,
"allCellTypes_front", allCellTypes)
2133 cellTypes_.transfer(allCellTypes);
2134 cellStencil_.setSize(mesh_.nCells());
2135 cellInterpolationWeights_.setSize(mesh_.nCells());
2137 forAll(cellTypes_, cellI)
2139 if (cellTypes_[cellI] == INTERPOLATED)
2141 cellStencil_[cellI].transfer(allStencil[cellI]);
2142 cellInterpolationWeights_[cellI].setSize(1);
2143 cellInterpolationWeights_[cellI][0] = 1.0;
2144 interpolationCells.
append(cellI);
2148 cellStencil_[cellI].clear();
2149 cellInterpolationWeights_[cellI].clear();
2152 interpolationCells_.transfer(interpolationCells);
2155 cellInterpolationMap_.reset
2159 cellInterpolationWeight_.
transfer(allWeight);
2164 >(cellInterpolationWeight_.boundaryFieldRef(),
false);
2170 mesh_.time().
write();
2173 mkDir(mesh_.time().timePath());
2174 OBJstream str(mesh_.time().timePath()/
"injectionStencil.obj");
2175 Pout<<
type() <<
" : dumping injectionStencil to "
2176 << str.name() <<
endl;
2178 cellInterpolationMap().distribute(cc);
2180 forAll(cellStencil_, celli)
2182 const labelList& slots = cellStencil_[celli];
2185 const point& accCc = mesh_.cellCentres()[celli];
2188 const point& donorCc = cc[slots[i]];
2197 createStencil(globalCells);
2203 cellInterpolationWeight_.instance() = mesh_.time().timeName();
2204 cellInterpolationWeight_.write();
2209 forAll(cellStencil_, celli)
2211 const scalarList& wghts = cellInterpolationWeights_[celli];
2214 if (
mag(wghts[i]) >
mag(maxMagWeight[celli]))
2216 maxMagWeight[celli] = wghts[i];
2219 if (
mag(maxMagWeight[celli]) > 1)
2222 Pout<<
"cell:" << celli
2223 <<
" at:" << cc[celli]
2224 <<
" zone:" <<
zoneID[celli]
2225 <<
" donors:" << cellStencil_[celli]
2226 <<
" weights:" << wghts
2236 createField(mesh_,
"maxMagWeight", maxMagWeight)
2242 >(tfld.
ref().boundaryFieldRef(),
false);
2250 createField(mesh_,
"cellTypes", cellTypes_)
2257 >(tfld.
ref().boundaryFieldRef(),
false);
2263 mkDir(mesh_.time().timePath());
2265 cellInterpolationMap().distribute(cc);
2270 mesh_.time().timePath()
2271 +
"/stencil_" +
name(zonei) +
".obj"
2273 Pout<<
type() <<
" : dumping to " << str.name() <<
endl;
2277 forAll(subMeshCellMap, subcelli)
2279 const label celli = subMeshCellMap[subcelli];
2280 const labelList& slots = cellStencil_[celli];
2281 const point& accCc = mesh_.cellCentres()[celli];
2284 const point& donorCc = cc[slots[i]];
2298 forAll(interpolationCells_, i)
2300 label celli = interpolationCells_[i];
2301 const labelList& slots = cellStencil_[celli];
2303 bool hasLocal =
false;
2304 bool hasRemote =
false;
2308 if (slots[sloti] >= mesh_.nCells())
2338 Info<<
"Overset analysis : nCells : "
2341 <<
indent <<
"calculated : " << nCells[CALCULATED] <<
nl
2342 <<
indent <<
"interpolated : " << nCells[INTERPOLATED]
2343 <<
" (interpolated from local:" << nLocal
2344 <<
" mixed local/remote:" << nMixed
2345 <<
" remote:" << nRemote <<
")" <<
nl
2346 <<
indent <<
"hole : " << nCells[HOLE] <<
nl