41void Foam::conformalVoronoiMesh::calcDualMesh
44 labelList& boundaryPts,
50 pointField& cellCentres,
51 labelList& cellToDelaunayVertex,
52 labelListList& patchToDelaunayVertex,
53 bitSet& boundaryFacesToRemove
58 setVertexSizeAndAlignment();
60 timeCheck(
"After setVertexSizeAndAlignment");
62 indexDualVertices(
points, boundaryPts);
65 Info<<
nl <<
"Merging identical points" <<
endl;
68 mergeIdenticalDualVertices(
points, boundaryPts);
73 timeCheck(
"Before createFacesOwnerNeighbourAndPatches");
75 createFacesOwnerNeighbourAndPatches
83 patchToDelaunayVertex,
84 boundaryFacesToRemove,
92 cellToDelaunayVertex = removeUnusedCells(owner, neighbour);
94 cellCentres =
pointField(cellCentres, cellToDelaunayVertex);
96 removeUnusedPoints(faces,
points, boundaryPts);
102void Foam::conformalVoronoiMesh::calcTetMesh
113 labelList vertexMap(number_of_vertices());
118 pointToDelaunayVertex.setSize(number_of_vertices());
122 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
123 vit != finite_vertices_end();
127 if (vit->internalPoint() || vit->boundaryPoint())
129 vertexMap[vit->index()] = vertI;
131 pointToDelaunayVertex[vertI] = vit->index();
137 pointToDelaunayVertex.setSize(vertI);
143 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
144 cit != finite_cells_end();
148 if (cit->internalOrBoundaryDualVertex())
150 cit->cellIndex() = celli++;
158 patchNames = geometryToConformTo_.patchNames();
166 List<DynamicList<face>> patchFaces(
nPatches);
167 List<DynamicList<label>> patchOwners(
nPatches);
169 faces.setSize(number_of_finite_facets());
171 owner.setSize(number_of_finite_facets());
173 neighbour.setSize(number_of_finite_facets());
177 labelList verticesOnTriFace(3, label(-1));
179 face newFace(verticesOnTriFace);
183 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
184 fit != finite_facets_end();
188 const Cell_handle
c1(fit->first);
189 const label oppositeVertex = fit->second;
190 const Cell_handle
c2(
c1->neighbor(oppositeVertex));
192 if (
c1->hasFarPoint() &&
c2->hasFarPoint())
198 label c1I =
c1->cellIndex();
199 label c2I =
c2->cellIndex();
201 label ownerCell = -1;
202 label neighbourCell = -1;
204 for (label i = 0; i < 3; i++)
206 verticesOnTriFace[i] = vertexMap
208 c1->vertex(vertex_triple_index(oppositeVertex, i))->index()
212 newFace = face(verticesOnTriFace);
214 if (
c1->hasFarPoint() ||
c2->hasFarPoint())
217 if (
c1->hasFarPoint())
230 label patchIndex = geometryToConformTo_.findPatch
235 if (patchIndex == -1)
241 <<
"did not find a surface patch. Adding to "
246 patchFaces[patchIndex].append(newFace);
247 patchOwners[patchIndex].append(ownerCell);
267 faces[facei] = newFace;
268 owner[facei] = ownerCell;
269 neighbour[facei] = neighbourCell;
274 label nInternalFaces = facei;
276 faces.setSize(nInternalFaces);
277 owner.setSize(nInternalFaces);
278 neighbour.setSize(nInternalFaces);
280 sortFaces(faces, owner, neighbour);
299void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
307 label nPtsMerged = 0;
308 label nPtsMergedSum = 0;
312 Map<label> dualPtIndexMap;
314 nPtsMerged = mergeIdenticalDualVertices
320 reindexDualVertices(dualPtIndexMap, boundaryPts);
322 reduce(nPtsMerged, sumOp<label>());
324 nPtsMergedSum += nPtsMerged;
326 }
while (nPtsMerged > 0);
328 if (nPtsMergedSum > 0)
330 Info<<
" Merged " << nPtsMergedSum <<
" points " <<
endl;
335Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
338 Map<label>& dualPtIndexMap
341 label nPtsMerged = 0;
345 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
346 fit != finite_facets_end();
350 const Cell_handle
c1(fit->first);
351 const label oppositeVertex = fit->second;
352 const Cell_handle
c2(
c1->neighbor(oppositeVertex));
354 if (is_infinite(c1) || is_infinite(c2))
359 label& c1I =
c1->cellIndex();
360 label& c2I =
c2->cellIndex();
362 if ((c1I != c2I) && !
c1->hasFarPoint() && !
c2->hasFarPoint())
384 dualPtIndexMap.insert(c1I, c1I);
385 dualPtIndexMap.insert(c2I, c1I);
389 dualPtIndexMap.insert(c1I, c2I);
390 dualPtIndexMap.insert(c2I, c2I);
400 Info<<
"mergeIdenticalDualVertices:" <<
endl
401 <<
" zero-length edges : "
674void Foam::conformalVoronoiMesh::deferredCollapseFaceSet
681 DynamicList<label> faceLabels;
685 if (deferredCollapseFaces.found(Pair<label>(owner[nI], neighbour[nI])))
687 faceLabels.append(nI);
691 Pout<<
"facesToCollapse" <<
nl << faceLabels <<
endl;
696Foam::conformalVoronoiMesh::createPolyMeshFromPoints
708 bitSet boundaryFacesToRemove;
710 timeCheck(
"Start of checkPolyMeshQuality");
712 Info<<
nl <<
"Creating polyMesh to assess quality" <<
endl;
714 createFacesOwnerNeighbourAndPatches
722 patchToDelaunayVertex,
723 boundaryFacesToRemove,
729 labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
730 cellCentres =
pointField(cellCentres, cellToDelaunayVertex);
736 "foamyHexMesh_temporary",
752 label nValidPatches = 0;
756 label totalPatchSize =
patchDicts[
p].get<label>(
"nFaces");
765 if (totalPatchSize > 0)
767 patchDicts[
p].set(
"transform",
"coincidentFullMatch");
769 patches[nValidPatches] =
new processorPolyPatch
774 pMesh.boundaryMesh(),
784 reduce(totalPatchSize, sumOp<label>());
786 if (totalPatchSize > 0)
809void Foam::conformalVoronoiMesh::checkCellSizing()
811 Info<<
"Checking cell sizes..."<<
endl;
813 timeCheck(
"Start of Cell Sizing");
815 labelList boundaryPts(number_of_finite_cells(), internal);
818 indexDualVertices(ptsField, boundaryPts);
821 mergeIdenticalDualVertices(ptsField, boundaryPts);
823 autoPtr<polyMesh>
meshPtr = createPolyMeshFromPoints(ptsField);
824 const polyMesh& pMesh =
meshPtr();
829 DynamicList<label> checkFaces(
identity(pMesh.nFaces()));
832 Info<<
"Running checkMesh on mesh with " << pMesh.nCells()
835 const dictionary&
dict
836 = foamyHexMeshControls().foamyHexMeshDict();
838 const dictionary& meshQualityDict
841 const scalar maxNonOrtho =
844 label nWrongFaces = 0;
846 if (maxNonOrtho < 180.0 - SMALL)
860 label nNonOrthogonal =
returnReduce(wrongFaces.size(), sumOp<label>());
862 Info<<
" non-orthogonality > " << maxNonOrtho
863 <<
" degrees : " << nNonOrthogonal <<
endl;
865 nWrongFaces += nNonOrthogonal;
868 labelHashSet protrudingCells = findOffsetPatchFaces(pMesh, 0.25);
870 label nProtrudingCells = protrudingCells.size();
872 Info<<
" protruding/intruding cells : " << nProtrudingCells <<
endl;
874 nWrongFaces += nProtrudingCells;
885 Info<<
" Found total of " << nWrongFaces <<
" bad faces" <<
endl;
891 for (
const label facei : wrongFaces)
893 const label faceOwner = pMesh.faceOwner()[facei];
894 const label faceNeighbour = pMesh.faceNeighbour()[facei];
896 cellsToResizeMap.insert(faceOwner);
897 cellsToResizeMap.insert(faceNeighbour);
900 cellsToResizeMap += protrudingCells;
902 pointField cellsToResize(cellsToResizeMap.size());
905 for (label celli = 0; celli < pMesh.nCells(); ++celli)
907 if (cellsToResizeMap.found(celli))
909 cellsToResize[
count++] = pMesh.cellCentres()[celli];
913 Info<<
" DISABLED: Automatically re-sizing " << cellsToResize.size()
914 <<
" cells that are attached to the bad faces: " <<
endl;
919 timeCheck(
"End of Cell Sizing");
921 Info<<
"Finished checking cell sizes"<<
endl;
927 const polyMesh&
mesh,
928 const scalar allowedOffset
931 timeCheck(
"Start findRemainingProtrusionSet");
935 cellSet offsetBoundaryCells
938 "foamyHexMesh_protrudingCells",
953 const face&
f = localFaces[pLFI];
957 const scalar targetSize = targetCellSize(faceCentre);
962 geometryToConformTo_.findSurfaceNearest
973 && (
mag(pHit.hitPoint() - faceCentre) > allowedOffset*targetSize)
976 offsetBoundaryCells.insert(fCell[pLFI]);
981 if (foamyHexMeshControls().objOutput())
983 offsetBoundaryCells.write();
986 return std::move(offsetBoundaryCells);
995 autoPtr<polyMesh>
meshPtr = createPolyMeshFromPoints(pts);
998 timeCheck(
"polyMesh created, checking quality");
1002 DynamicList<label> checkFaces(pMesh.nFaces());
1006 scalar faceAreaLimit = SMALL;
1010 if (
mag(fAreas[fI]) > faceAreaLimit)
1012 checkFaces.append(fI);
1016 Info<<
nl <<
"Excluding "
1017 <<
returnReduce(fAreas.size() - checkFaces.size(), sumOp<label>())
1018 <<
" faces from check, < " << faceAreaLimit <<
" area" <<
endl;
1020 const dictionary&
dict
1021 = foamyHexMeshControls().foamyHexMeshDict();
1023 const dictionary& meshQualityDict
1037 label nInvalidPolyhedra = 0;
1043 if (
cells[cI].size() < 4 &&
cells[cI].size() > 0)
1049 nInvalidPolyhedra++;
1051 wrongFaces.insert(
cells[cI]);
1055 Info<<
" cells with more than 1 but fewer than 4 faces : "
1063 for (label fI = 0; fI < pMesh.nInternalFaces(); fI++)
1065 nInternalFaces[pMesh.faceOwner()[fI]]++;
1066 nInternalFaces[pMesh.faceNeighbour()[fI]]++;
1069 const polyBoundaryMesh&
patches = pMesh.boundaryMesh();
1079 nInternalFaces[owners[i]]++;
1084 label oneInternalFaceCells = 0;
1086 forAll(nInternalFaces, cI)
1088 if (nInternalFaces[cI] <= 1)
1090 oneInternalFaceCells++;
1091 wrongFaces.insert(
cells[cI]);
1095 Info<<
" cells with with zero or one non-boundary face : "
1101 bitSet ptToBeLimited(pts.size(),
false);
1103 for (
const label facei : wrongFaces)
1105 const face
f = pMesh.faces()[facei];
1107 ptToBeLimited.
set(
f);
1141 label maxFilterCount = 0;
1145 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1146 cit != finite_cells_end();
1150 label cI = cit->cellIndex();
1154 if (ptToBeLimited.test(cI))
1156 cit->filterCount()++;
1159 if (cit->filterCount() > maxFilterCount)
1161 maxFilterCount = cit->filterCount();
1166 Info<<
nl <<
"Maximum number of filter limits applied: "
1173Foam::label Foam::conformalVoronoiMesh::classifyBoundaryPoint
1178 if (cit->boundaryDualVertex())
1180 if (cit->featurePointDualVertex())
1182 return featurePoint;
1184 else if (cit->featureEdgeDualVertex())
1193 else if (cit->baffleSurfaceDualVertex())
1197 else if (cit->baffleEdgeDualVertex())
1208void Foam::conformalVoronoiMesh::indexDualVertices
1216 this->resetCellCount();
1218 label nConstrainedVertices = 0;
1219 if (foamyHexMeshControls().guardFeaturePoints())
1223 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1224 vit != finite_vertices_end();
1228 if (vit->constrained())
1230 vit->index() = number_of_finite_cells() + nConstrainedVertices;
1231 nConstrainedVertices++;
1236 pts.setSize(number_of_finite_cells() + nConstrainedVertices);
1239 number_of_finite_cells() + nConstrainedVertices,
1243 if (foamyHexMeshControls().guardFeaturePoints())
1245 nConstrainedVertices = 0;
1248 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1249 vit != finite_vertices_end();
1253 if (vit->constrained())
1255 pts[number_of_finite_cells() + nConstrainedVertices] =
1258 boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
1261 nConstrainedVertices++;
1272 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1273 cit != finite_cells_end();
1284 if (!cit->hasFarPoint())
1286 cit->cellIndex() = getNewCellIndex();
1294 typedef CGAL::Exact_predicates_exact_constructions_kernel Exact;
1295 typedef CGAL::Point_3<Exact> ExactPoint;
1297 List<labelPair> cellVerticesPair(4);
1298 List<ExactPoint> cellVertices(4);
1300 for (label vI = 0; vI < 4; ++vI)
1304 cit->vertex(vI)->procIndex(),
1305 cit->vertex(vI)->index()
1308 cellVertices[vI] = ExactPoint
1310 cit->vertex(vI)->point().x(),
1311 cit->vertex(vI)->point().y(),
1312 cit->vertex(vI)->point().z()
1319 oldToNew =
invert(oldToNew.size(), oldToNew);
1322 ExactPoint synchronisedDual = CGAL::circumcenter
1332 CGAL::to_double(synchronisedDual.x()),
1333 CGAL::to_double(synchronisedDual.y()),
1334 CGAL::to_double(synchronisedDual.z())
1339 pts[cit->cellIndex()] = cit->dual();
1343 if (foamyHexMeshControls().snapFeaturePoints())
1345 if (cit->featurePointDualVertex())
1353 geometryToConformTo_.findFeaturePointNearest
1356 sqr(targetCellSize(dual)),
1365 Info<<
"Dual = " << dual <<
nl
1366 <<
" Nearest = " << fpHit.hitPoint() <<
endl;
1369 pts[cit->cellIndex()] = fpHit.hitPoint();
1456 boundaryPts[cit->cellIndex()] = classifyBoundaryPoint(cit);
1470void Foam::conformalVoronoiMesh::reindexDualVertices
1472 const Map<label>& dualPtIndexMap,
1478 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1479 cit != finite_cells_end();
1483 if (dualPtIndexMap.found(cit->cellIndex()))
1485 cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
1486 boundaryPts[cit->cellIndex()] =
1489 boundaryPts[cit->cellIndex()],
1490 boundaryPts[dualPtIndexMap[cit->cellIndex()]]
1497Foam::label Foam::conformalVoronoiMesh::createPatchInfo
1503 patchNames = geometryToConformTo_.patchNames();
1507 const PtrList<dictionary>& patchInfo = geometryToConformTo_.patchInfo();
1511 if (patchInfo.set(patchi))
1513 patchDicts.set(patchi,
new dictionary(patchInfo[patchi]));
1528 patchNames[defaultPatchIndex] =
"foamyHexMesh_defaultPatch";
1529 patchDicts.set(defaultPatchIndex,
new dictionary());
1540 List<boolList> procUsedList
1551 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1552 vit != finite_vertices_end();
1559 if (vit->referred())
1561 procUsed[vit->procIndex()] =
true;
1568 forAll(procUsedList, proci)
1574 procUsed[proci] =
true;
1592 for (label pI = nNonProcPatches; pI < nTotalPatches; ++pI)
1618 patchDicts[nNonProcPatches + procAddI].set(
"neighbProcNo", pUI);
1625 return defaultPatchIndex;
1629Foam::vector Foam::conformalVoronoiMesh::calcSharedPatchNormal
1638 for (label cI = 0; cI < 4; ++cI)
1640 if (
c1->neighbor(cI) != c2 && !
c1->vertex(cI)->constrained())
1642 if (
c1->vertex(cI)->internalBoundaryPoint())
1644 patchEdge[0] =
topoint(
c1->vertex(cI)->point());
1648 patchEdge[1] =
topoint(
c1->vertex(cI)->point());
1655 return vector(patchEdge[1] - patchEdge[0]);
1659bool Foam::conformalVoronoiMesh::boundaryDualFace
1665 label nInternal = 0;
1666 label nExternal = 0;
1668 for (label cI = 0; cI < 4; ++cI)
1670 if (
c1->neighbor(cI) != c2 && !
c1->vertex(cI)->constrained())
1672 if (
c1->vertex(cI)->internalBoundaryPoint())
1676 else if (
c1->vertex(cI)->externalBoundaryPoint())
1683 Info<<
"in = " << nInternal <<
" out = " << nExternal <<
endl;
1685 return (nInternal == 1 && nExternal == 1);
1689void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
1698 bitSet& boundaryFacesToRemove,
1699 bool includeEmptyPatches
1707 forAll(procNeighbours, patchi)
1709 procNeighbours[patchi] =
1710 patchDicts[patchi].getOrDefault<label>(
"neighbProcNo", -1);
1713 List<DynamicList<face>> patchFaces(
nPatches);
1714 List<DynamicList<label>> patchOwners(
nPatches);
1716 List<DynamicList<label>> patchPPSlaves(
nPatches);
1717 List<DynamicList<bool>> indirectPatchFace(
nPatches);
1720 faces.setSize(number_of_finite_edges());
1721 owner.setSize(number_of_finite_edges());
1722 neighbour.setSize(number_of_finite_edges());
1723 boundaryFacesToRemove.setSize(number_of_finite_edges(),
false);
1725 labelPairPairDynListList procPatchSortingIndex(
nPatches);
1727 label dualFacei = 0;
1729 if (foamyHexMeshControls().guardFeaturePoints())
1731 OBJstream startCellStr(
"startingCell.obj");
1732 OBJstream featurePointFacesStr(
"ftPtFaces.obj");
1733 OBJstream featurePointDualsStr(
"ftPtDuals.obj");
1734 OFstream cellStr(
"vertexCells.obj");
1740 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1741 vit != finite_vertices_end();
1745 if (vit->constrained())
1748 std::list<Cell_handle> vertexCells;
1749 finite_incident_cells(vit, std::back_inserter(vertexCells));
1751 Cell_handle startCell;
1755 std::list<Cell_handle>::iterator vcit = vertexCells.begin();
1756 vcit != vertexCells.end();
1760 if ((*vcit)->featurePointExternalCell())
1765 if ((*vcit)->real())
1767 featurePointDualsStr.write
1775 if (startCell ==
nullptr)
1777 Pout<<
"Start cell is null!" <<
endl;
1781 Cell_handle vc1 = startCell;
1784 Info<<
"c1 index = " << vc1->cellIndex() <<
" "
1785 << vc1->dual() <<
endl;
1787 for (label cI = 0; cI < 4; ++cI)
1789 Info<<
"c1 = " << cI <<
" "
1790 << vc1->neighbor(cI)->cellIndex() <<
" v = "
1791 << vc1->neighbor(cI)->dual() <<
endl;
1793 Info<< vc1->vertex(cI)->info();
1796 Cell_handle nextCell;
1798 for (label cI = 0; cI < 4; ++cI)
1800 if (vc1->vertex(cI)->externalBoundaryPoint())
1802 vc2 = vc1->neighbor(cI);
1804 Info<<
" c2 is neighbor "
1806 <<
" of c1" <<
endl;
1808 for (label cI = 0; cI < 4; ++cI)
1810 Info<<
" c2 = " << cI <<
" "
1811 << vc2->neighbor(cI)->cellIndex() <<
" v = "
1812 << vc2->vertex(cI)->index() <<
endl;
1816 f[0] = vit->index();
1817 f[1] = vc1->cellIndex();
1818 f[2] = vc2->cellIndex();
1826 vector correctNormal = calcSharedPatchNormal(vc1, vc2);
1827 correctNormal.normalise();
1829 Info<<
" cN " << correctNormal <<
endl;
1831 vector fN =
f.areaNormal(pts);
1833 if (
mag(fN) < SMALL)
1842 if ((fN & correctNormal) > 0)
1852 label own = vit->index();
1856 Info<<
"Start walk from " << vc1->cellIndex()
1857 <<
" to " << vc2->cellIndex() <<
endl;
1864 Info<<
" Walk from " << vc1->cellIndex()
1865 <<
" " << vc1->dual()
1866 <<
" to " << vc2->cellIndex()
1867 <<
" " << vc2->dual()
1875 geometryToConformTo_.findPatch
1880 f[1] = vc1->cellIndex();
1881 f[2] = vc2->cellIndex();
1883 patchFaces[patchIndex].
append(
f);
1884 patchOwners[patchIndex].append(own);
1885 patchPPSlaves[patchIndex].append(own);
1888 Cell_handle nextCell;
1890 Info<<
" c1 vertices " << vc2->dual() <<
endl;
1891 for (label cI = 0; cI < 4; ++cI)
1893 Info<<
" " << vc2->vertex(cI)->info();
1895 Info<<
" c1 neighbour vertices " <<
endl;
1896 for (label cI = 0; cI < 4; ++cI)
1900 !vc2->vertex(cI)->constrained()
1901 && vc2->neighbor(cI) != vc1
1902 && !is_infinite(vc2->neighbor(cI))
1905 vc2->neighbor(cI)->featurePointExternalCell()
1906 || vc2->neighbor(cI)->featurePointInternalCell()
1908 && vc2->neighbor(cI)->hasConstrainedPoint()
1918 Info<<
" neighbour " << cI <<
" "
1919 << vc2->neighbor(cI)->dual() <<
endl;
1920 for (label
I = 0;
I < 4; ++
I)
1923 << vc2->neighbor(cI)->vertex(
I)->info();
1928 for (label cI = 0; cI < 4; ++cI)
1932 !vc2->vertex(cI)->constrained()
1933 && vc2->neighbor(cI) != vc1
1934 && !is_infinite(vc2->neighbor(cI))
1937 vc2->neighbor(cI)->featurePointExternalCell()
1938 || vc2->neighbor(cI)->featurePointInternalCell()
1940 && vc2->neighbor(cI)->hasConstrainedPoint()
1944 if (boundaryDualFace(vc2, vc2->neighbor(cI)))
1946 nextCell = vc2->neighbor(cI);
1956 }
while (vc1 != startCell && iter < 100);
1963 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1964 eit != finite_edges_end();
1968 Cell_handle
c = eit->first;
1969 Vertex_handle vA =
c->vertex(eit->second);
1970 Vertex_handle vB =
c->vertex(eit->third);
1972 if (vA->constrained() && vB->constrained())
1979 (vA->constrained() && vB->internalOrBoundaryPoint())
1980 || (vB->constrained() && vA->internalOrBoundaryPoint())
1983 face newDualFace = buildDualFace(eit);
1988 if (ownerAndNeighbour(vA, vB, own, nei))
1994 faces[dualFacei] = newDualFace;
1995 owner[dualFacei] = own;
1996 neighbour[dualFacei] = nei;
2002 (vA->internalOrBoundaryPoint() && !vA->referred())
2003 || (vB->internalOrBoundaryPoint() && !vB->referred())
2008 (vA->internalPoint() && vB->externalBoundaryPoint())
2009 || (vB->internalPoint() && vA->externalBoundaryPoint())
2012 Cell_circulator ccStart = incident_cells(*eit);
2013 Cell_circulator cc1 = ccStart;
2014 Cell_circulator cc2 = cc1;
2018 bool skipEdge =
false;
2024 cc1->hasFarPoint() || cc2->hasFarPoint()
2025 || is_infinite(cc1) || is_infinite(cc2)
2028 Pout<<
"Ignoring edge between internal and external: "
2039 }
while (cc1 != ccStart);
2054 face newDualFace = buildDualFace(eit);
2056 if (newDualFace.size() >= 3)
2061 if (ownerAndNeighbour(vA, vB, own, nei))
2066 label patchIndex = -1;
2075 if (isProcBoundaryEdge(eit))
2080 label procIndex =
max(vA->procIndex(), vB->procIndex());
2084 procNeighbours.find(vA->procIndex()),
2085 procNeighbours.find(vB->procIndex())
2096 DynamicList<labelPairPair>& sortingIndex =
2097 procPatchSortingIndex[patchIndex];
2099 if (vB->internalOrBoundaryPoint() && vB->referred())
2105 labelPair(vA->index(), vA->procIndex()),
2116 labelPair(vB->index(), vB->procIndex()),
2127 DynamicList<labelPairPair>& sortingIndex =
2128 procPatchSortingIndex[patchIndex];
2130 if (vA->internalOrBoundaryPoint() && vA->referred())
2136 labelPair(vA->index(), vA->procIndex()),
2147 labelPair(vB->index(), vB->procIndex()),
2166 patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
2169 if (patchIndex == -1)
2178 patchIndex = geometryToConformTo_.findPatch
2184 patchFaces[patchIndex].append(newDualFace);
2185 patchOwners[patchIndex].append(own);
2191 vA->boundaryPoint() && vB->boundaryPoint()
2192 && !ptPairs_.isPointPair(vA, vB)
2193 && !ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
2196 indirectPatchFace[patchIndex].append(
true);
2200 indirectPatchFace[patchIndex].append(
false);
2204 if (vA->internalOrBoundaryPoint())
2206 patchPPSlaves[patchIndex].append(vB->index());
2210 patchPPSlaves[patchIndex].append(vA->index());
2217 !vA->boundaryPoint()
2218 || !vB->boundaryPoint()
2219 || ptPairs_.isPointPair(vA, vB)
2220 || ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
2223 patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
2229 && geometryToConformTo_.patchInfo().set(patchIndex)
2234 patchFaces[patchIndex].append(newDualFace);
2235 patchOwners[patchIndex].append(own);
2236 indirectPatchFace[patchIndex].append(
false);
2240 patchFaces[patchIndex].append(newDualFace);
2241 patchOwners[patchIndex].append(nei);
2242 indirectPatchFace[patchIndex].append(
false);
2247 <
labelPair(vA->index(), vA->procIndex())
2250 patchPPSlaves[patchIndex].append(vB->index());
2251 patchPPSlaves[patchIndex].append(vB->index());
2255 patchPPSlaves[patchIndex].append(vA->index());
2256 patchPPSlaves[patchIndex].append(vA->index());
2263 faces[dualFacei] = newDualFace;
2264 owner[dualFacei] = own;
2265 neighbour[dualFacei] = nei;
2274 if (!patchFaces[defaultPatchIndex].empty())
2276 Pout<<
nl << patchFaces[defaultPatchIndex].size()
2277 <<
" faces were not able to have their patch determined from "
2279 <<
nl <<
"Adding to patch " <<
patchNames[defaultPatchIndex]
2283 label nInternalFaces = dualFacei;
2285 faces.setSize(nInternalFaces);
2286 owner.setSize(nInternalFaces);
2287 neighbour.setSize(nInternalFaces);
2289 timeCheck(
"polyMesh quality checked");
2291 sortFaces(faces, owner, neighbour);
2298 procPatchSortingIndex
2301 timeCheck(
"faces, owner, neighbour sorted");
2309 boundaryFacesToRemove,
2316 patchPointPairSlaves.setSize(
nPatches);
2317 forAll(patchPPSlaves, patchi)
2319 patchPointPairSlaves[patchi].transfer(patchPPSlaves[patchi]);
2322 if (foamyHexMeshControls().objOutput())
2324 Info<<
"Writing processor interfaces" <<
endl;
2328 if (patchFaces[nbI].size() > 0)
2330 const label neighbour =
2331 patchDicts[nbI].getOrDefault<label>(
"neighbProcNo", -1);
2333 faceList procPatchFaces = patchFaces[nbI];
2339 forAll(procPatchFaces, fI)
2341 procPatchFaces[fI].flip();
2345 if (neighbour != -1)
2356 time().
path()/fName,
2367void Foam::conformalVoronoiMesh::sortFaces
2388 List<labelPair> ownerNeighbourPair(owner.size());
2390 forAll(ownerNeighbourPair, oNI)
2392 ownerNeighbourPair[oNI] =
labelPair(owner[oNI], neighbour[oNI]);
2396 <<
"Sorting faces, owner and neighbour into upper triangular order"
2400 oldToNew =
invert(oldToNew.size(), oldToNew);
2408void Foam::conformalVoronoiMesh::sortProcPatches
2410 List<DynamicList<face>>& patchFaces,
2411 List<DynamicList<label>>& patchOwners,
2412 List<DynamicList<label>>& patchPointPairSlaves,
2413 labelPairPairDynListList& patchSortingIndices
2421 forAll(patchSortingIndices, patchi)
2423 faceList& faces = patchFaces[patchi];
2425 DynamicList<label>& slaves = patchPointPairSlaves[patchi];
2426 DynamicList<labelPairPair>& sortingIndices
2427 = patchSortingIndices[patchi];
2429 if (!sortingIndices.empty())
2433 faces.size() != sortingIndices.size()
2434 || owner.size() != sortingIndices.size()
2435 || slaves.size() != sortingIndices.size()
2439 <<
"patch size and size of sorting indices is inconsistent "
2440 <<
" for patch " << patchi <<
nl
2441 <<
" faces.size() " << faces.size() <<
nl
2442 <<
" owner.size() " << owner.size() <<
nl
2443 <<
" slaves.size() " << slaves.size() <<
nl
2444 <<
" sortingIndices.size() "
2445 << sortingIndices.size()
2450 oldToNew =
invert(oldToNew.size(), oldToNew);
2461void Foam::conformalVoronoiMesh::addPatches
2463 const label nInternalFaces,
2467 bitSet& boundaryFacesToRemove,
2468 const List<DynamicList<face>>& patchFaces,
2469 const List<DynamicList<label>>& patchOwners,
2470 const List<DynamicList<bool>>& indirectPatchFace
2473 label nBoundaryFaces = 0;
2478 patchDicts[
p].set(
"startFace", nInternalFaces + nBoundaryFaces);
2480 nBoundaryFaces += patchFaces[
p].size();
2483 faces.setSize(nInternalFaces + nBoundaryFaces);
2484 owner.setSize(nInternalFaces + nBoundaryFaces);
2485 boundaryFacesToRemove.setSize(nInternalFaces + nBoundaryFaces);
2487 label facei = nInternalFaces;
2493 faces[facei] = patchFaces[
p][
f];
2494 owner[facei] = patchOwners[
p][
f];
2495 boundaryFacesToRemove[facei] = indirectPatchFace[
p][
f];
2503void Foam::conformalVoronoiMesh::removeUnusedPoints
2510 Info<<
nl <<
"Removing unused points" <<
endl;
2512 bitSet ptUsed(pts.size(),
false);
2518 const face&
f = faces[fI];
2525 labelList oldToNew(pts.size(), label(-1));
2532 if (ptUsed.test(ptUI))
2534 oldToNew[ptUI] = pointi++;
2546 pts.setSize(pointi);
2547 boundaryPts.setSize(pointi);
2564 Info<<
nl <<
"Removing unused cells" <<
endl;
2566 bitSet cellUsed(vertexCount(),
false);
2570 cellUsed.set(owner);
2571 cellUsed.set(neighbour);
2575 labelList oldToNew(cellUsed.size(), label(-1));
2582 if (cellUsed.test(cellUI))
2584 oldToNew[cellUI] = celli++;
2594 DynamicList<label> unusedCells;
2598 if (!cellUsed.test(cUI))
2600 unusedCells.append(cUI);
2604 if (unusedCells.size() > 0)
2608 <<
" unused cell labels" <<
endl;
2612 label& o = owner[oI];
2619 label&
n = neighbour[nI];
Various functions to operate on Lists.
InfoProxy< IOstream > info() const
Return info proxy.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
void setSize(const label n)
Alias for resize()
void append(const T &val)
Append an element at the end of the list.
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
void setSize(const label newLen)
Same as resize()
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
static bool & parRun() noexcept
Test if this a parallel run.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
void checkMesh() const
Debug: Check coupled mesh for correctness.
static const complex max
complex (VGREAT,VGREAT)
UPtrList< const labelUList > faceCells() const
Return a list of faceCells for each patch.
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
label nCells() const noexcept
Number of mesh cells.
int myProcNo() const noexcept
Return processor number.
static word newName(const label myProcNo, const label neighbProcNo)
Return the name of a processorPolyPatch.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
const polyBoundaryMesh & patches
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
Foam::autoPtr< Foam::dynamicFvMesh > meshPtr
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define WarningInFunction
Report a warning using Foam::Warning.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
const dimensionedScalar c
Speed of light in a vacuum.
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
const wordList internal
Standard dimensioned field types (scalar, vector, tensor, etc)
const std::string patch
OpenFOAM patch number as a std::string.
Pair< label > labelPair
A pair of labels.
List< word > wordList
A List of words.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
List< cell > cellList
A List of cells.
line< point, const point & > linePointRef
A line using referred points.
pointFromPoint topoint(const Point &P)
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
static const Identity< scalar > I
label findLower(const ListType &input, const T &val, const label start, const ComparePredicate &comp)
vector point
Point is a vector.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Pair< labelPair > labelPairPair
A pair of labelPairs.
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< bool > boolList
A List of bools.
static constexpr const zero Zero
Global zero (0)
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
UList< label > labelUList
A UList of labels.
HashSet< labelPair, Foam::Hash< labelPair > > labelPairHashSet
A HashSet for a labelPair. The hashing is based on labelPair (FixedList) and is thus non-commutative.
constexpr char nl
The newline '\n' character (0x0a)
wordList patchNames(nPatches)
PtrList< dictionary > patchDicts
#define forAll(list, i)
Loop across all elements in list.
static const char *const typeName
The type name used in ensight case files.