38 void Foam::conformalVoronoiMesh::calcDualMesh
46 PtrList<dictionary>& patchDicts,
50 PackedBoolList& boundaryFacesToRemove
55 setVertexSizeAndAlignment();
57 timeCheck(
"After setVertexSizeAndAlignment");
59 indexDualVertices(points, boundaryPts);
62 Info<<
nl <<
"Merging identical points" <<
endl;
65 mergeIdenticalDualVertices(points, boundaryPts);
70 timeCheck(
"Before createFacesOwnerNeighbourAndPatches");
72 createFacesOwnerNeighbourAndPatches
80 patchToDelaunayVertex,
81 boundaryFacesToRemove,
89 cellToDelaunayVertex = removeUnusedCells(owner, neighbour);
91 cellCentres =
pointField(cellCentres, cellToDelaunayVertex);
93 removeUnusedPoints(faces, points, boundaryPts);
99 void Foam::conformalVoronoiMesh::calcTetMesh
107 PtrList<dictionary>& patchDicts
110 labelList vertexMap(number_of_vertices());
114 points.setSize(number_of_vertices());
115 pointToDelaunayVertex.setSize(number_of_vertices());
119 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
120 vit != finite_vertices_end();
124 if (vit->internalPoint() || vit->boundaryPoint())
126 vertexMap[vit->index()] = vertI;
127 points[vertI] =
topoint(vit->point());
128 pointToDelaunayVertex[vertI] = vit->index();
133 points.setSize(vertI);
134 pointToDelaunayVertex.setSize(vertI);
140 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
141 cit != finite_cells_end();
145 if (cit->internalOrBoundaryDualVertex())
147 cit->cellIndex() = celli++;
155 patchNames = geometryToConformTo_.
patchNames();
157 patchNames.
setSize(patchNames.size() + 1);
159 patchNames[patchNames.size() - 1] =
"foamyHexMesh_defaultPatch";
161 label nPatches = patchNames.size();
163 List<DynamicList<face>> patchFaces(nPatches, DynamicList<face>(0));
165 List<DynamicList<label>> patchOwners(nPatches, DynamicList<label>(0));
167 faces.setSize(number_of_finite_facets());
169 owner.setSize(number_of_finite_facets());
171 neighbour.setSize(number_of_finite_facets());
177 face newFace(verticesOnTriFace);
181 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
182 fit != finite_facets_end();
187 const label oppositeVertex = fit->second;
190 if (c1->hasFarPoint() && c2->hasFarPoint())
196 label c1I = c1->cellIndex();
197 label c2I = c2->cellIndex();
199 label ownerCell = -1;
200 label neighbourCell = -1;
202 for (
label i = 0; i < 3; i++)
204 verticesOnTriFace[i] = vertexMap
206 c1->vertex(vertex_triple_index(oppositeVertex, i))->index()
210 newFace = face(verticesOnTriFace);
212 if (c1->hasFarPoint() || c2->hasFarPoint())
215 if (c1->hasFarPoint())
230 newFace.centre(points)
233 if (patchIndex == -1)
235 patchIndex = patchNames.size() - 1;
238 << newFace.centre(points) <<
nl 239 <<
"did not find a surface patch. Adding to " 240 << patchNames[patchIndex]
244 patchFaces[patchIndex].append(newFace);
245 patchOwners[patchIndex].append(ownerCell);
265 faces[facei] = newFace;
266 owner[facei] = ownerCell;
267 neighbour[facei] = neighbourCell;
272 label nInternalFaces = facei;
274 faces.setSize(nInternalFaces);
275 owner.setSize(nInternalFaces);
276 neighbour.setSize(nInternalFaces);
278 sortFaces(faces, owner, neighbour);
297 void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
305 label nPtsMerged = 0;
306 label nPtsMergedSum = 0;
310 Map<label> dualPtIndexMap;
312 nPtsMerged = mergeIdenticalDualVertices
318 reindexDualVertices(dualPtIndexMap, boundaryPts);
320 reduce(nPtsMerged, sumOp<label>());
322 nPtsMergedSum += nPtsMerged;
324 }
while (nPtsMerged > 0);
326 if (nPtsMergedSum > 0)
328 Info<<
" Merged " << nPtsMergedSum <<
" points " <<
endl;
333 Foam::label Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
336 Map<label>& dualPtIndexMap
339 label nPtsMerged = 0;
343 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
344 fit != finite_facets_end();
349 const label oppositeVertex = fit->second;
352 if (is_infinite(c1) || is_infinite(c2))
357 label& c1I = c1->cellIndex();
358 label& c2I = c2->cellIndex();
360 if ((c1I != c2I) && !c1->hasFarPoint() && !c2->hasFarPoint())
382 dualPtIndexMap.insert(c1I, c1I);
383 dualPtIndexMap.insert(c2I, c1I);
387 dualPtIndexMap.insert(c1I, c2I);
388 dualPtIndexMap.insert(c2I, c2I);
398 Info<<
"mergeIdenticalDualVertices:" << endl
399 <<
" zero-length edges : " 672 void Foam::conformalVoronoiMesh::deferredCollapseFaceSet
676 const HashSet<
labelPair, labelPair::Hash<>>& deferredCollapseFaces
679 DynamicList<label> faceLabels;
683 if (deferredCollapseFaces.found(Pair<label>(owner[nI], neighbour[nI])))
685 faceLabels.append(nI);
689 Pout<<
"facesToCollapse" <<
nl << faceLabels <<
endl;
694 Foam::conformalVoronoiMesh::createPolyMeshFromPoints
706 PackedBoolList boundaryFacesToRemove;
708 timeCheck(
"Start of checkPolyMeshQuality");
710 Info<<
nl <<
"Creating polyMesh to assess quality" <<
endl;
712 createFacesOwnerNeighbourAndPatches
720 patchToDelaunayVertex,
721 boundaryFacesToRemove,
727 labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
728 cellCentres =
pointField(cellCentres, cellToDelaunayVertex);
736 "foamyHexMesh_temporary",
751 List<polyPatch*>
patches(patchNames.size());
753 label nValidPatches = 0;
762 && word(patchDicts[
p].
lookup(
"type")) == processorPolyPatch::typeName
766 if (totalPatchSize > 0)
768 patchDicts[
p].set(
"transform",
"coincidentFullMatch");
770 patches[nValidPatches] =
new processorPolyPatch
775 pMesh.boundaryMesh(),
776 processorPolyPatch::typeName
785 reduce(totalPatchSize, sumOp<label>());
787 if (totalPatchSize > 0)
802 patches.setSize(nValidPatches);
804 pMesh.addPatches(patches);
810 void Foam::conformalVoronoiMesh::checkCellSizing()
812 Info<<
"Checking cell sizes..."<<
endl;
816 labelList boundaryPts(number_of_finite_cells(),
internal);
819 indexDualVertices(ptsField, boundaryPts);
822 mergeIdenticalDualVertices(ptsField, boundaryPts);
824 autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(ptsField);
825 const polyMesh& pMesh =
meshPtr();
830 DynamicList<label> checkFaces(
identity(pMesh.nFaces()));
833 Info<<
"Running checkMesh on mesh with " << pMesh.nCells()
836 const dictionary& dict
839 const dictionary& meshQualityDict
840 = dict.
subDict(
"meshQualityControls");
842 const scalar maxNonOrtho =
843 readScalar(meshQualityDict.lookup(
"maxNonOrtho",
true));
845 label nWrongFaces = 0;
847 if (maxNonOrtho < 180.0 - small)
863 Info<<
" non-orthogonality > " << maxNonOrtho
864 <<
" degrees : " << nNonOrthogonal <<
endl;
866 nWrongFaces += nNonOrthogonal;
869 labelHashSet protrudingCells = findOffsetPatchFaces(pMesh, 0.25);
871 label nProtrudingCells = protrudingCells.size();
873 Info<<
" protruding/intruding cells : " << nProtrudingCells <<
endl;
875 nWrongFaces += nProtrudingCells;
886 Info<<
" Found total of " << nWrongFaces <<
" bad faces" <<
endl;
894 const label faceOwner = pMesh.faceOwner()[iter.key()];
895 const label faceNeighbour = pMesh.faceNeighbour()[iter.key()];
897 if (!cellsToResizeMap.found(faceOwner))
899 cellsToResizeMap.insert(faceOwner);
902 if (!cellsToResizeMap.found(faceNeighbour))
904 cellsToResizeMap.insert(faceNeighbour);
908 cellsToResizeMap += protrudingCells;
910 pointField cellsToResize(cellsToResizeMap.size());
913 for (
label celli = 0; celli < pMesh.nCells(); ++celli)
915 if (cellsToResizeMap.found(celli))
917 cellsToResize[count++] = pMesh.cellCentres()[celli];
921 Info<<
" DISABLED: Automatically re-sizing " << cellsToResize.size()
922 <<
" cells that are attached to the bad faces: " <<
endl;
929 Info<<
"Finished checking cell sizes"<<
endl;
935 const polyMesh& mesh,
936 const scalar allowedOffset
939 timeCheck(
"Start findRemainingProtrusionSet");
941 const polyBoundaryMesh& patches = mesh.boundaryMesh();
943 cellSet offsetBoundaryCells
946 "foamyHexMesh_protrudingCells",
952 const polyPatch& patch = patches[
patchi];
954 const faceList& localFaces = patch.localFaces();
955 const pointField& localPoints = patch.localPoints();
957 const labelList& fCell = patch.faceCells();
961 const face& f = localFaces[pLFI];
965 const scalar targetSize = targetCellSize(faceCentre);
981 && (
mag(pHit.hitPoint() - faceCentre) > allowedOffset*targetSize)
984 offsetBoundaryCells.insert(fCell[pLFI]);
991 offsetBoundaryCells.write();
994 return offsetBoundaryCells;
1003 autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(pts);
1006 timeCheck(
"polyMesh created, checking quality");
1010 DynamicList<label> checkFaces(pMesh.nFaces());
1014 scalar faceAreaLimit = small;
1018 if (
mag(fAreas[fI]) > faceAreaLimit)
1020 checkFaces.append(fI);
1024 Info<<
nl <<
"Excluding " 1025 <<
returnReduce(fAreas.size() - checkFaces.size(), sumOp<label>())
1026 <<
" faces from check, < " << faceAreaLimit <<
" area" <<
endl;
1028 const dictionary& dict
1031 const dictionary& meshQualityDict
1032 = dict.
subDict(
"meshQualityControls");
1045 label nInvalidPolyhedra = 0;
1047 const cellList& cells = pMesh.cells();
1051 if (cells[cI].size() < 4 && cells[cI].size() > 0)
1057 nInvalidPolyhedra++;
1061 wrongFaces.insert(cells[cI][cFI]);
1066 Info<<
" cells with more than 1 but fewer than 4 faces : " 1074 for (
label fI = 0; fI < pMesh.nInternalFaces(); fI++)
1076 nInternalFaces[pMesh.faceOwner()[fI]]++;
1077 nInternalFaces[pMesh.faceNeighbour()[fI]]++;
1080 const polyBoundaryMesh& patches = pMesh.boundaryMesh();
1084 if (patches[
patchi].coupled())
1090 nInternalFaces[owners[i]]++;
1095 label oneInternalFaceCells = 0;
1097 forAll(nInternalFaces, cI)
1099 if (nInternalFaces[cI] <= 1)
1101 oneInternalFaceCells++;
1105 wrongFaces.insert(cells[cI][cFI]);
1110 Info<<
" cells with with zero or one non-boundary face : " 1116 PackedBoolList ptToBeLimited(pts.size(),
false);
1120 const face f = pMesh.faces()[iter.key()];
1124 ptToBeLimited[f[fPtI]] =
true;
1169 label maxFilterCount = 0;
1173 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1174 cit != finite_cells_end();
1178 label cI = cit->cellIndex();
1182 if (ptToBeLimited[cI] ==
true)
1184 cit->filterCount()++;
1187 if (cit->filterCount() > maxFilterCount)
1189 maxFilterCount = cit->filterCount();
1194 Info<<
nl <<
"Maximum number of filter limits applied: " 1201 Foam::label Foam::conformalVoronoiMesh::classifyBoundaryPoint
1206 if (cit->boundaryDualVertex())
1208 if (cit->featurePointDualVertex())
1212 else if (cit->featureEdgeDualVertex())
1221 else if (cit->baffleSurfaceDualVertex())
1225 else if (cit->baffleEdgeDualVertex())
1236 void Foam::conformalVoronoiMesh::indexDualVertices
1246 label nConstrainedVertices = 0;
1251 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1252 vit != finite_vertices_end();
1256 if (vit->constrained())
1258 vit->index() = number_of_finite_cells() + nConstrainedVertices;
1259 nConstrainedVertices++;
1264 pts.setSize(number_of_finite_cells() + nConstrainedVertices);
1267 number_of_finite_cells() + nConstrainedVertices,
1273 nConstrainedVertices = 0;
1276 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1277 vit != finite_vertices_end();
1281 if (vit->constrained())
1283 pts[number_of_finite_cells() + nConstrainedVertices] =
1286 boundaryPts[number_of_finite_cells() + nConstrainedVertices] =
1289 nConstrainedVertices++;
1300 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1301 cit != finite_cells_end();
1312 if (!cit->hasFarPoint())
1322 typedef CGAL::Exact_predicates_exact_constructions_kernel Exact;
1323 typedef CGAL::Point_3<Exact> ExactPoint;
1325 List<labelPair> cellVerticesPair(4);
1326 List<ExactPoint> cellVertices(4);
1328 for (
label vI = 0; vI < 4; ++vI)
1332 cit->vertex(vI)->procIndex(),
1333 cit->vertex(vI)->index()
1336 cellVertices[vI] = ExactPoint
1338 cit->vertex(vI)->point().x(),
1339 cit->vertex(vI)->point().y(),
1340 cit->vertex(vI)->point().z()
1348 oldToNew =
invert(oldToNew.size(), oldToNew);
1351 ExactPoint synchronisedDual = CGAL::circumcenter
1361 CGAL::to_double(synchronisedDual.x()),
1362 CGAL::to_double(synchronisedDual.y()),
1363 CGAL::to_double(synchronisedDual.z())
1368 pts[cit->cellIndex()] = cit->dual();
1374 if (cit->featurePointDualVertex())
1385 sqr(targetCellSize(dual)),
1394 Info<<
"Dual = " << dual <<
nl 1395 <<
" Nearest = " << fpHit.hitPoint() <<
endl;
1398 pts[cit->cellIndex()] = fpHit.hitPoint();
1485 boundaryPts[cit->cellIndex()] = classifyBoundaryPoint(cit);
1499 void Foam::conformalVoronoiMesh::reindexDualVertices
1501 const Map<label>& dualPtIndexMap,
1507 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1508 cit != finite_cells_end();
1512 if (dualPtIndexMap.found(cit->cellIndex()))
1514 cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
1515 boundaryPts[cit->cellIndex()] =
1518 boundaryPts[cit->cellIndex()],
1519 boundaryPts[dualPtIndexMap[cit->cellIndex()]]
1526 Foam::label Foam::conformalVoronoiMesh::createPatchInfo
1529 PtrList<dictionary>& patchDicts
1532 patchNames = geometryToConformTo_.
patchNames();
1534 patchDicts.
setSize(patchNames.size() + 1);
1536 const PtrList<dictionary>& patchInfo = geometryToConformTo_.
patchInfo();
1540 if (patchInfo.set(
patchi))
1542 patchDicts.set(
patchi,
new dictionary(patchInfo[
patchi]));
1546 patchDicts.set(
patchi,
new dictionary());
1550 wallPolyPatch::typeName
1555 patchNames.setSize(patchNames.size() + 1);
1556 label defaultPatchIndex = patchNames.size() - 1;
1557 patchNames[defaultPatchIndex] =
"foamyHexMesh_defaultPatch";
1558 patchDicts.set(defaultPatchIndex,
new dictionary());
1559 patchDicts[defaultPatchIndex].set
1562 wallPolyPatch::typeName
1565 label nProcPatches = 0;
1569 List<boolList> procUsedList
1580 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1581 vit != finite_vertices_end();
1588 if (vit->referred())
1590 procUsed[vit->procIndex()] =
true;
1598 forAll(procUsedList, proci)
1604 procUsed[proci] =
true;
1617 label nNonProcPatches = patchNames.size();
1618 label nTotalPatches = nNonProcPatches + nProcPatches;
1620 patchNames.setSize(nTotalPatches);
1621 patchDicts.setSize(nTotalPatches);
1622 for (
label pI = nNonProcPatches; pI < nTotalPatches; ++pI)
1624 patchDicts.set(pI,
new dictionary());
1633 patchNames[nNonProcPatches + procAddI] =
1636 patchDicts[nNonProcPatches + procAddI].set
1639 processorPolyPatch::typeName
1642 patchDicts[nNonProcPatches + procAddI].set
1648 patchDicts[nNonProcPatches + procAddI].set(
"neighbProcNo", pUI);
1655 return defaultPatchIndex;
1659 Foam::vector Foam::conformalVoronoiMesh::calcSharedPatchNormal
1668 for (
label cI = 0; cI < 4; ++cI)
1670 if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
1672 if (c1->vertex(cI)->internalBoundaryPoint())
1674 patchEdge[0] =
topoint(c1->vertex(cI)->point());
1678 patchEdge[1] =
topoint(c1->vertex(cI)->point());
1685 return vector(patchEdge[1] - patchEdge[0]);
1689 bool Foam::conformalVoronoiMesh::boundaryDualFace
1695 label nInternal = 0;
1696 label nExternal = 0;
1698 for (
label cI = 0; cI < 4; ++cI)
1700 if (c1->neighbor(cI) != c2 && !c1->vertex(cI)->constrained())
1702 if (c1->vertex(cI)->internalBoundaryPoint())
1706 else if (c1->vertex(cI)->externalBoundaryPoint())
1713 Info<<
"in = " << nInternal <<
" out = " << nExternal <<
endl;
1715 return (nInternal == 1 && nExternal == 1);
1719 void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
1726 PtrList<dictionary>& patchDicts,
1728 PackedBoolList& boundaryFacesToRemove,
1729 bool includeEmptyPatches
1732 const label defaultPatchIndex = createPatchInfo(patchNames, patchDicts);
1734 const label nPatches = patchNames.size();
1743 patchDicts[
patchi].found(
"neighbProcNo")
1750 List<DynamicList<face>> patchFaces(nPatches, DynamicList<face>(0));
1751 List<DynamicList<label>> patchOwners(nPatches, DynamicList<label>(0));
1753 List<DynamicList<label>> patchPPSlaves(nPatches, DynamicList<label>(0));
1755 List<DynamicList<bool>> indirectPatchFace(nPatches, DynamicList<bool>(0));
1758 faces.setSize(number_of_finite_edges());
1759 owner.setSize(number_of_finite_edges());
1760 neighbour.setSize(number_of_finite_edges());
1761 boundaryFacesToRemove.setSize(number_of_finite_edges(),
false);
1765 label dualFacei = 0;
1769 OBJstream startCellStr(
"startingCell.obj");
1770 OBJstream featurePointFacesStr(
"ftPtFaces.obj");
1771 OBJstream featurePointDualsStr(
"ftPtDuals.obj");
1772 OFstream cellStr(
"vertexCells.obj");
1778 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1779 vit != finite_vertices_end();
1783 if (vit->constrained())
1786 std::list<Cell_handle> vertexCells;
1787 finite_incident_cells(vit, std::back_inserter(vertexCells));
1793 std::list<Cell_handle>::iterator vcit = vertexCells.begin();
1794 vcit != vertexCells.end();
1798 if ((*vcit)->featurePointExternalCell())
1803 if ((*vcit)->real())
1805 featurePointDualsStr.write
1813 if (startCell ==
nullptr)
1815 Pout<<
"Start cell is null!" <<
endl;
1822 Info<<
"c1 index = " << vc1->cellIndex() <<
" " 1823 << vc1->dual() <<
endl;
1825 for (
label cI = 0; cI < 4; ++cI)
1827 Info<<
"c1 = " << cI <<
" " 1828 << vc1->neighbor(cI)->cellIndex() <<
" v = " 1829 << vc1->neighbor(cI)->dual() <<
endl;
1831 Info<< vc1->vertex(cI)->info();
1836 for (
label cI = 0; cI < 4; ++cI)
1838 if (vc1->vertex(cI)->externalBoundaryPoint())
1840 vc2 = vc1->neighbor(cI);
1842 Info<<
" c2 is neighbor " 1844 <<
" of c1" <<
endl;
1846 for (
label cI = 0; cI < 4; ++cI)
1848 Info<<
" c2 = " << cI <<
" " 1849 << vc2->neighbor(cI)->cellIndex() <<
" v = " 1850 << vc2->vertex(cI)->index() <<
endl;
1854 f[0] = vit->index();
1855 f[1] = vc1->cellIndex();
1856 f[2] = vc2->cellIndex();
1864 vector correctNormal = calcSharedPatchNormal(vc1, vc2);
1865 correctNormal /=
mag(correctNormal);
1867 Info<<
" cN " << correctNormal <<
endl;
1871 if (
mag(fN) < small)
1880 if ((fN & correctNormal) > 0)
1890 label own = vit->index();
1894 Info<<
"Start walk from " << vc1->cellIndex()
1895 <<
" to " << vc2->cellIndex() <<
endl;
1902 Info<<
" Walk from " << vc1->cellIndex()
1903 <<
" " << vc1->dual()
1904 <<
" to " << vc2->cellIndex()
1905 <<
" " << vc2->dual()
1918 f[1] = vc1->cellIndex();
1919 f[2] = vc2->cellIndex();
1921 patchFaces[patchIndex].append(f);
1922 patchOwners[patchIndex].append(own);
1923 patchPPSlaves[patchIndex].append(own);
1928 Info<<
" c1 vertices " << vc2->dual() <<
endl;
1929 for (
label cI = 0; cI < 4; ++cI)
1931 Info<<
" " << vc2->vertex(cI)->info();
1933 Info<<
" c1 neighbour vertices " <<
endl;
1934 for (
label cI = 0; cI < 4; ++cI)
1938 !vc2->vertex(cI)->constrained()
1939 && vc2->neighbor(cI) != vc1
1940 && !is_infinite(vc2->neighbor(cI))
1943 vc2->neighbor(cI)->featurePointExternalCell()
1944 || vc2->neighbor(cI)->featurePointInternalCell()
1946 && vc2->neighbor(cI)->hasConstrainedPoint()
1956 Info<<
" neighbour " << cI <<
" " 1957 << vc2->neighbor(cI)->dual() <<
endl;
1961 << vc2->neighbor(cI)->vertex(
I)->info();
1966 for (
label cI = 0; cI < 4; ++cI)
1970 !vc2->vertex(cI)->constrained()
1971 && vc2->neighbor(cI) != vc1
1972 && !is_infinite(vc2->neighbor(cI))
1975 vc2->neighbor(cI)->featurePointExternalCell()
1976 || vc2->neighbor(cI)->featurePointInternalCell()
1978 && vc2->neighbor(cI)->hasConstrainedPoint()
1982 if (boundaryDualFace(vc2, vc2->neighbor(cI)))
1984 nextCell = vc2->neighbor(cI);
1994 }
while (vc1 != startCell && iter < 100);
2001 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
2002 eit != finite_edges_end();
2010 if (vA->constrained() && vB->constrained())
2017 (vA->constrained() && vB->internalOrBoundaryPoint())
2018 || (vB->constrained() && vA->internalOrBoundaryPoint())
2021 face newDualFace = buildDualFace(eit);
2026 if (ownerAndNeighbour(vA, vB, own, nei))
2032 faces[dualFacei] = newDualFace;
2033 owner[dualFacei] = own;
2034 neighbour[dualFacei] = nei;
2040 (vA->internalOrBoundaryPoint() && !vA->referred())
2041 || (vB->internalOrBoundaryPoint() && !vB->referred())
2046 (vA->internalPoint() && vB->externalBoundaryPoint())
2047 || (vB->internalPoint() && vA->externalBoundaryPoint())
2050 Cell_circulator ccStart = incident_cells(*eit);
2051 Cell_circulator cc1 = ccStart;
2052 Cell_circulator cc2 = cc1;
2056 bool skipEdge =
false;
2062 cc1->hasFarPoint() || cc2->hasFarPoint()
2063 || is_infinite(cc1) || is_infinite(cc2)
2066 Pout<<
"Ignoring edge between internal and external: " 2077 }
while (cc1 != ccStart);
2092 face newDualFace = buildDualFace(eit);
2094 if (newDualFace.size() >= 3)
2099 if (ownerAndNeighbour(vA, vB, own, nei))
2104 label patchIndex = -1;
2113 if (isProcBoundaryEdge(eit))
2118 label procIndex =
max(vA->procIndex(), vB->procIndex());
2122 findIndex(procNeighbours, vA->procIndex()),
2123 findIndex(procNeighbours, vB->procIndex())
2134 DynamicList<Pair<labelPair>>& sortingIndex =
2135 procPatchSortingIndex[patchIndex];
2137 if (vB->internalOrBoundaryPoint() && vB->referred())
2143 labelPair(vA->index(), vA->procIndex()),
2154 labelPair(vB->index(), vB->procIndex()),
2165 DynamicList<Pair<labelPair>>& sortingIndex =
2166 procPatchSortingIndex[patchIndex];
2168 if (vA->internalOrBoundaryPoint() && vA->referred())
2174 labelPair(vA->index(), vA->procIndex()),
2185 labelPair(vB->index(), vB->procIndex()),
2204 patchIndex = geometryToConformTo_.
findPatch(ptA, ptB);
2207 if (patchIndex == -1)
2216 patchIndex = geometryToConformTo_.
findPatch 2222 patchFaces[patchIndex].append(newDualFace);
2223 patchOwners[patchIndex].append(own);
2229 vA->boundaryPoint() && vB->boundaryPoint()
2234 indirectPatchFace[patchIndex].append(
true);
2238 indirectPatchFace[patchIndex].append(
false);
2242 if (vA->internalOrBoundaryPoint())
2244 patchPPSlaves[patchIndex].append(vB->index());
2248 patchPPSlaves[patchIndex].append(vA->index());
2255 !vA->boundaryPoint()
2256 || !vB->boundaryPoint()
2261 patchIndex = geometryToConformTo_.
findPatch(ptA, ptB);
2267 && geometryToConformTo_.
patchInfo().set(patchIndex)
2272 patchFaces[patchIndex].append(newDualFace);
2273 patchOwners[patchIndex].append(own);
2274 indirectPatchFace[patchIndex].append(
false);
2278 patchFaces[patchIndex].append(newDualFace);
2279 patchOwners[patchIndex].append(nei);
2280 indirectPatchFace[patchIndex].append(
false);
2285 <
labelPair(vA->index(), vA->procIndex())
2288 patchPPSlaves[patchIndex].append(vB->index());
2289 patchPPSlaves[patchIndex].append(vB->index());
2293 patchPPSlaves[patchIndex].append(vA->index());
2294 patchPPSlaves[patchIndex].append(vA->index());
2301 faces[dualFacei] = newDualFace;
2302 owner[dualFacei] = own;
2303 neighbour[dualFacei] = nei;
2312 if (!patchFaces[defaultPatchIndex].empty())
2314 Pout<<
nl << patchFaces[defaultPatchIndex].size()
2315 <<
" faces were not able to have their patch determined from " 2317 <<
nl <<
"Adding to patch " << patchNames[defaultPatchIndex]
2321 label nInternalFaces = dualFacei;
2323 faces.setSize(nInternalFaces);
2324 owner.setSize(nInternalFaces);
2325 neighbour.setSize(nInternalFaces);
2329 sortFaces(faces, owner, neighbour);
2336 procPatchSortingIndex
2339 timeCheck(
"faces, owner, neighbour sorted");
2347 boundaryFacesToRemove,
2354 patchPointPairSlaves.setSize(nPatches);
2357 patchPointPairSlaves[
patchi].transfer(patchPPSlaves[
patchi]);
2362 Info<<
"Writing processor interfaces" <<
endl;
2366 if (patchFaces[nbI].size() > 0)
2368 const label neighbour =
2370 patchDicts[nbI].found(
"neighbProcNo")
2375 faceList procPatchFaces = patchFaces[nbI];
2381 forAll(procPatchFaces, fI)
2383 procPatchFaces[fI] = procPatchFaces[fI].reverseFace();
2387 if (neighbour != -1)
2409 void Foam::conformalVoronoiMesh::sortFaces
2430 List<labelPair> ownerNeighbourPair(owner.size());
2432 forAll(ownerNeighbourPair, oNI)
2434 ownerNeighbourPair[oNI] =
labelPair(owner[oNI], neighbour[oNI]);
2438 <<
"Sorting faces, owner and neighbour into upper triangular order" 2445 oldToNew =
invert(oldToNew.size(), oldToNew);
2453 void Foam::conformalVoronoiMesh::sortProcPatches
2455 List<DynamicList<face>>& patchFaces,
2456 List<DynamicList<label>>& patchOwners,
2457 List<DynamicList<label>>& patchPointPairSlaves,
2470 DynamicList<label>& slaves = patchPointPairSlaves[
patchi];
2471 DynamicList<Pair<labelPair>>& sortingIndices
2472 = patchSortingIndices[
patchi];
2474 if (!sortingIndices.empty())
2478 faces.size() != sortingIndices.size()
2479 || owner.size() != sortingIndices.size()
2480 || slaves.size() != sortingIndices.size()
2484 <<
"patch size and size of sorting indices is inconsistent " 2486 <<
" faces.size() " << faces.size() <<
nl 2487 <<
" owner.size() " << owner.size() <<
nl 2488 <<
" slaves.size() " << slaves.size() <<
nl 2489 <<
" sortingIndices.size() " 2490 << sortingIndices.size()
2498 oldToNew =
invert(oldToNew.size(), oldToNew);
2509 void Foam::conformalVoronoiMesh::addPatches
2511 const label nInternalFaces,
2514 PtrList<dictionary>& patchDicts,
2515 PackedBoolList& boundaryFacesToRemove,
2516 const List<DynamicList<face>>& patchFaces,
2517 const List<DynamicList<label>>& patchOwners,
2518 const List<DynamicList<bool>>& indirectPatchFace
2521 label nBoundaryFaces = 0;
2525 patchDicts[
p].set(
"nFaces", patchFaces[
p].size());
2526 patchDicts[
p].set(
"startFace", nInternalFaces + nBoundaryFaces);
2528 nBoundaryFaces += patchFaces[
p].size();
2531 faces.setSize(nInternalFaces + nBoundaryFaces);
2532 owner.setSize(nInternalFaces + nBoundaryFaces);
2533 boundaryFacesToRemove.setSize(nInternalFaces + nBoundaryFaces);
2535 label facei = nInternalFaces;
2541 faces[facei] = patchFaces[
p][
f];
2542 owner[facei] = patchOwners[
p][
f];
2543 boundaryFacesToRemove[facei] = indirectPatchFace[
p][
f];
2551 void Foam::conformalVoronoiMesh::removeUnusedPoints
2558 Info<<
nl <<
"Removing unused points" <<
endl;
2560 PackedBoolList ptUsed(pts.size(),
false);
2566 const face& f = faces[fI];
2570 ptUsed[f[fPtI]] =
true;
2583 if (ptUsed[ptUI] ==
true)
2585 oldToNew[ptUI] = pointi++;
2597 pts.setSize(pointi);
2598 boundaryPts.setSize(pointi);
2615 Info<<
nl <<
"Removing unused cells" <<
endl;
2623 cellUsed[owner[oI]] =
true;
2628 cellUsed[neighbour[nI]] =
true;
2640 if (cellUsed[cellUI] ==
true)
2642 oldToNew[cellUI] = celli++;
2652 DynamicList<label> unusedCells;
2656 if (cellUsed[cUI] ==
false)
2658 unusedCells.append(cUI);
2662 if (unusedCells.size() > 0)
2666 <<
" unused cell labels" <<
endl;
2670 label& o = owner[oI];
2677 label&
n = neighbour[nI];
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
List< labelList > labelListList
A List of labelList.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
A HashTable with keys but without contents.
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
pointFromPoint topoint(const Point &P)
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
static const Vector< scalar > max
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
static word newName(const label myProcNo, const label neighbProcNo)
Return the name of a processorPolyPatch.
Ostream & endl(Ostream &os)
Add newline and flush stream.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Xfer< T > xferCopy(const T &)
Construct by copying the contents of the arg.
PointIndexHit< point > pointIndexHit
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Vector< scalar > vector
A scalar version of the templated Vector.
const dictionary & foamyHexMeshDict() const
Return the foamyHexMeshDict.
UList< label > labelUList
Various functions to operate on Lists.
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
List< bool > boolList
Bool container classes.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
vectorField pointField
pointField is a vectorField.
stressControl lookup("compactNormalStress") >> compactNormalStress
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
line< point, const point & > linePointRef
Line using referred points.
static const Identity< scalar > I
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
wordList patchNames(nPatches)
Pair< label > labelPair
Label pair.
List< label > labelList
A List of labels.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
label readLabel(Istream &is)
void reverse(UList< T > &, const label n)
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
bool isPointPair(const Vertex_handle &vA, const Vertex_handle &vB) const
List< word > wordList
A List of words.
void setSize(const label)
Reset size of List.
void resetCellCount()
Set the cell count to zero.
static bool & parRun()
Is this a parallel run?
static label nProcs(const label communicator=0)
Number of processes in parallel run.
label vertexCount() const
Return the vertex count (the next unique vertex index)
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
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.
label getNewCellIndex() const
Create a new unique cell index and return.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static autoPtr< polyPatch > New(const word &patchType, const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm)
Return a pointer to a new patch created on freestore from.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
virtual Ostream & write(const token &)=0
Write next token to stream.
PtrList< dictionary > patchDicts
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
List< cell > cellList
list of cells
InfoProxy< IOstream > info() const
Return info proxy.