46 namespace cellCellStencils
61 const scalar layerRelax,
67 bitSet isFront(mesh_.nFaces());
77 if (isA<oversetFvPatch>(fvm[patchI]))
82 isFront.set(fvm[patchI].start()+i);
91 const labelList& nei = mesh_.faceNeighbour();
93 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
95 label ownType = allCellTypes[own[faceI]];
96 label neiType = allCellTypes[nei[faceI]];
99 (ownType == HOLE && neiType != HOLE)
100 || (ownType != HOLE && neiType == HOLE)
110 syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
114 label faceI = mesh_.nInternalFaces();
115 faceI < mesh_.nFaces();
119 label ownType = allCellTypes[own[faceI]];
120 label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];
124 (ownType == HOLE && neiType != HOLE)
125 || (ownType != HOLE && neiType == HOLE)
138 scalar fraction = 1.0;
144 Info<<
"Front : fraction:" << fraction
148 bitSet newIsFront(mesh_.nFaces());
151 if (isFront.test(faceI))
153 label own = mesh_.faceOwner()[faceI];
154 if (allCellTypes[own] != HOLE)
156 if (allWeight[own] < fraction)
158 allWeight[own] = fraction;
166 allCellTypes[own] = INTERPOLATED;
167 newIsFront.
set(mesh_.cells()[own]);
170 if (mesh_.isInternalFace(faceI))
172 label nei = mesh_.faceNeighbour()[faceI];
173 if (allCellTypes[nei] != HOLE)
175 if (allWeight[nei] < fraction)
177 allWeight[nei] = fraction;
186 allCellTypes[nei] = INTERPOLATED;
187 newIsFront.
set(mesh_.cells()[nei]);
196 isFront.transfer(newIsFront);
198 fraction -= layerRelax;
228 for (label faceI = 0; faceI <
mesh.nInternalFaces(); faceI++)
234 (ownType == HOLE && neiType != HOLE)
235 || (ownType != HOLE && neiType == HOLE)
238 isBlockedFace[faceI] =
true;
243 syncTools::swapBoundaryCellList(
mesh,
cellTypes, nbrCellTypes);
245 for (label faceI =
mesh.nInternalFaces(); faceI <
mesh.nFaces(); faceI++)
248 label neiType = nbrCellTypes[faceI-
mesh.nInternalFaces()];
252 (ownType == HOLE && neiType != HOLE)
253 || (ownType != HOLE && neiType == HOLE)
256 isBlockedFace[faceI] =
true;
262 Info<< typeName <<
" : detected " << cellRegion.
nRegions()
263 <<
" mesh regions after overset" <<
nl <<
endl;
281 for (label faceI = 0; faceI <
mesh.nInternalFaces(); faceI++)
283 if (isBlockedFace[faceI])
285 label ownRegion = cellRegion[own[faceI]];
289 if (regionType[ownRegion] == 0)
295 regionType[ownRegion] = 1;
299 label neiRegion = cellRegion[nei[faceI]];
303 if (regionType[neiRegion] == 0)
309 regionType[neiRegion] = 1;
316 label faceI =
mesh.nInternalFaces();
317 faceI <
mesh.nFaces();
321 if (isBlockedFace[faceI])
323 label ownRegion = cellRegion[own[faceI]];
325 if (regionType[ownRegion] == 0)
331 regionType[ownRegion] = 1;
341 const fvPatch& fvp = pbm[patchI];
343 if (isA<oversetFvPatch>(fvp))
345 else if (!fvPatch::constraintType(fvp.type()))
353 label regionI = cellRegion[fc[i]];
355 if (
cellTypes[fc[i]] != HOLE && regionType[regionI] != 2)
360 regionType[regionI] = 2;
381 Pstream::listCombineScatter(regionType);
391 const fvPatch& fvp = pbm[patchI];
393 if (isA<oversetFvPatch>(fvp))
399 label regionI = cellRegion[cellI];
401 if (regionType[regionI] != 2)
403 const labelList& slots = compactStencil[cellI];
406 label otherType = cellRegionType[slots[i]];
413 regionType[regionI] = 2;
435 label
type = regionType[cellRegion[cellI]];
455 const fvPatch& fvp = pbm[patchI];
458 if (isA<oversetFvPatch>(fvp))
464 patchCellTypes[cellMap[cellI]] = OVERSET;
467 else if (!fvPatch::constraintType(fvp.type()))
474 if (patchCellTypes[cellMap[cellI]] != OVERSET)
476 patchCellTypes[cellMap[cellI]] =
PATCH;
493 const labelList& slots = addressing[cellI];
501 result[cellI] = OVERSET;
507 result[cellI] =
PATCH;
509 else if (result[cellI] == -1)
512 result[cellI] = OTHER;
527 if (result.size() != addressing.size())
541 mapPtr().distribute(work);
543 interpolatePatchTypes(addressing, work, result);
547 interpolatePatchTypes(addressing,
patchTypes, result);
554 const label subZoneID,
558 const label donorZoneID,
562 const labelList& interpolatedOtherPatchTypes,
570 forAll(subCellMap, subCellI)
572 label cellI = subCellMap[subCellI];
574 bool validDonors =
true;
575 switch (interpolatedOtherPatchTypes[subCellI])
592 allCellTypes[cellI] = HOLE;
639 label currentDiff =
mag(subZoneID-allDonorID[cellI]);
640 label thisDiff =
mag(subZoneID-donorZoneID);
644 allDonorID[cellI] == -1
645 || (thisDiff < currentDiff)
646 || (thisDiff == currentDiff && donorZoneID > allDonorID[cellI])
649 allWeights[cellI] = weights[subCellI];
652 allDonorID[cellI] = donorZoneID;
661 Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
670 overlapTolerance_(defaultOverlapTolerance_),
672 interpolationCells_(0),
673 cellInterpolationMap_(),
675 cellInterpolationWeights_(0),
676 cellInterpolationWeight_
680 "cellInterpolationWeight",
681 mesh_.facesInstance(),
689 zeroGradientFvPatchScalarField::typeName
693 nonInterpolatedFields_.insert(
"cellTypes");
694 nonInterpolatedFields_.insert(
"cellInterpolationWeight");
697 nonInterpolatedFields_.insert(
"cellDisplacement");
698 nonInterpolatedFields_.insert(
"grad(cellDisplacement)");
699 const word w(
"snGradCorr(cellDisplacement)");
700 const word d(
"((viscosity*faceDiffusivity)*magSf)");
701 nonInterpolatedFields_.insert(
"surfaceIntegrate(("+d+
"*"+w+
"))");
710 mesh_.time().timeName(),
712 IOobject::READ_IF_PRESENT,
720 Pout<<
"Reading cellTypes from time " << mesh_.time().timeName()
725 forAll(volCellTypes, celli)
728 cellTypes_[celli] = volCellTypes[celli];
749 scalar layerRelax(dict_.getOrDefault(
"layerRelax", 1.0));
763 <<
" mesh regions" <<
nl <<
endl;
771 Info<<
indent<<
"zone:" << zonei <<
" nCells:"
786 labelList allCellTypes(mesh_.nCells(), CALCULATED);
787 labelList allPatchTypes(mesh_.nCells(), OTHER);
791 labelList allDonorID(mesh_.nCells(), -1);
804 Info<<
"Marking patch-cells on zone " << partI <<
endl;
805 markPatchCells(partMesh, partCellMap, allPatchTypes);
811 <<
"After patch analysis : nCells : "
814 <<
indent <<
"other : " << nCells[OTHER] <<
nl
816 <<
indent <<
"overset: " << nCells[OVERSET] <<
nl
822 for (label srcI = 0; srcI <
meshParts.size()-1; srcI++)
827 for (label tgtI = srcI+1; tgtI <
meshParts.size(); tgtI++)
847 interpolatePatchTypes
852 interpolatedTgtPatchTypes
859 forAll(tgtCellMap, tgtCellI)
861 label cellI = tgtCellMap[tgtCellI];
862 tgtGlobalCells[tgtCellI] = globalCells.
toGlobal(cellI);
864 if (mapper.
tgtMap().valid())
866 mapper.
tgtMap()().distribute(tgtGlobalCells);
879 interpolatedTgtPatchTypes,
892 interpolatePatchTypes
897 interpolatedSrcPatchTypes
902 forAll(srcCellMap, srcCellI)
904 label cellI = srcCellMap[srcCellI];
905 srcGlobalCells[srcCellI] = globalCells.
toGlobal(cellI);
907 if (mapper.
srcMap().valid())
909 mapper.
srcMap()().distribute(srcGlobalCells);
923 interpolatedSrcPatchTypes,
940 createField(mesh_,
"allCellTypes", allCellTypes)
947 forAll(allPatchTypes, cellI)
949 if (allCellTypes[cellI] != HOLE)
951 switch (allPatchTypes[cellI])
956 scalar v = mesh_.V()[cellI];
957 scalar overlapVol =
sum(allWeights[cellI]);
958 if (overlapVol > overlapTolerance_*v)
960 allCellTypes[cellI] = INTERPOLATED;
966 allCellTypes[cellI] = HOLE;
967 allWeights[cellI].
clear();
968 allStencil[cellI].
clear();
999 createField(mesh_,
"allCellTypes_patch", allCellTypes)
1011 label nCalculated = 0;
1013 forAll(cellTypes_, celli)
1015 if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
1017 if (allStencil[celli].size() == 0)
1020 allCellTypes[celli] = HOLE;
1021 allWeights[celli].
clear();
1022 allStencil[celli].
clear();
1026 allCellTypes[celli] = INTERPOLATED;
1034 Pout<<
"Detected " << nCalculated <<
" cells changing from hole"
1035 <<
" to calculated. Changed these to interpolated"
1042 findHoles(globalCells, mesh_,
zoneID, allStencil, allCellTypes);
1048 createField(mesh_,
"allCellTypes_hole", allCellTypes)
1057 walkFront(layerRelax, allCellTypes, allWeight);
1063 createField(mesh_,
"allCellTypes_front", allCellTypes)
1070 forAll(allCellTypes, cellI)
1072 if (allCellTypes[cellI] == INTERPOLATED)
1074 if (allWeight[cellI] < SMALL || allStencil[cellI].size() == 0)
1078 allWeights[cellI].
clear();
1079 allStencil[cellI].
clear();
1080 allWeight[cellI] = 0.0;
1084 scalar
s =
sum(allWeights[cellI]);
1085 forAll(allWeights[cellI], i)
1087 allWeights[cellI][i] /=
s;
1093 allWeights[cellI].
clear();
1094 allStencil[cellI].
clear();
1107 mesh_.time().timeName(),
1115 zeroGradientFvPatchScalarField::typeName
1137 mesh_.time().timeName(),
1145 zeroGradientFvPatchScalarField::typeName
1148 forAll(volTypes.internalField(), cellI)
1150 volTypes[cellI] = allCellTypes[cellI];
1157 >(volTypes.boundaryFieldRef(),
false);
1197 cellTypes_.transfer(allCellTypes);
1198 cellStencil_.transfer(allStencil);
1199 cellInterpolationWeights_.transfer(allWeights);
1200 cellInterpolationWeight_.transfer(allWeight);
1206 >(cellInterpolationWeight_.boundaryFieldRef(),
false);
1209 forAll(cellStencil_, cellI)
1211 if (cellStencil_[cellI].size())
1213 interpolationCells.append(cellI);
1216 interpolationCells_.transfer(interpolationCells);
1220 cellInterpolationMap_.reset
1229 cellInterpolationWeight_.instance() = mesh_.time().timeName();
1230 cellInterpolationWeight_.write();
1233 mkDir(mesh_.time().timePath());
1234 OBJstream str(mesh_.time().timePath()/
"stencil2.obj");
1235 Info<< typeName <<
" : dumping to " << str.name() <<
endl;
1237 cellInterpolationMap().distribute(cc);
1239 forAll(interpolationCells_, compactI)
1241 label cellI = interpolationCells_[compactI];
1242 const labelList& slots = cellStencil_[cellI];
1244 Pout<<
"cellI:" << cellI <<
" at:"
1245 << mesh_.cellCentres()[cellI]
1246 <<
" calculated from slots:" << slots
1248 <<
" weights:" << cellInterpolationWeights_[cellI]
1253 const point& donorCc = cc[slots[i]];
1254 const point& accCc = mesh_.cellCentres()[cellI];
1264 Info<<
"Overset analysis : nCells : "
1267 <<
indent <<
"calculated : " << nCells[CALCULATED] <<
nl
1268 <<
indent <<
"interpolated : " << nCells[INTERPOLATED] <<
nl
1269 <<
indent <<
"hole : " << nCells[HOLE] <<
nl
1280 const point& sample,
1287 weights.
setSize(donorCcs.size());
1291 scalar d =
mag(sample-donorCcs[i]);