68 SortableList<label> nbr(cFaces.size());
72 label facei = cFaces[i];
79 if (nbrCelli == celli)
108 oldToNew[cFaces[nbr.indices()[i]]] = newFacei++;
120 for (label facei = newFacei; facei <
faceOwner.size(); facei++)
122 oldToNew[facei] = facei;
129 if (oldToNew[facei] == -1)
132 <<
"Did not determine new position"
133 <<
" for face " << facei
147 void Foam::polyDualMesh::getPointEdges
157 const face&
f =
patch.localFaces()[facei];
164 label edgeI = fEdges[i];
166 const edge&
e =
patch.edges()[edgeI];
171 label index =
f.find(pointi);
173 if (
f.nextLabel(index) ==
e[1])
182 if (e0 != -1 && e1 != -1)
187 else if (
e[1] == pointi)
190 label index =
f.find(pointi);
192 if (
f.nextLabel(index) ==
e[0])
201 if (e0 != -1 && e1 != -1)
209 <<
" vertices:" <<
patch.localFaces()[facei]
210 <<
" that uses point:" << pointi
218 const polyPatch&
patch,
219 const label patchToDualOffset,
230 label meshPointi =
patch.meshPoints()[pointi];
233 DynamicList<label> dualFace;
235 if (pointToDualPoint[meshPointi] >= 0)
238 dualFace.setCapacity(
pFaces.size()+2+1);
240 dualFace.append(pointToDualPoint[meshPointi]);
244 dualFace.setCapacity(
pFaces.size()+2);
248 if (edgeToDualPoint[
patch.meshEdges()[edgeI]] < 0)
254 dualFace.append(edgeToDualPoint[
patch.meshEdges()[edgeI]]);
256 label facei =
patch.edgeFaces()[edgeI][0];
262 getPointEdges(
patch, facei, pointi, e0, e1);
276 dualFace.append(facei + patchToDualOffset);
280 getPointEdges(
patch, facei, pointi, e0, e1);
291 if (edgeToDualPoint[
patch.meshEdges()[edgeI]] >= 0)
294 dualFace.append(edgeToDualPoint[
patch.meshEdges()[edgeI]]);
299 if (eFaces.size() == 1)
306 if (eFaces[0] == facei)
330 void Foam::polyDualMesh::collectPatchInternalFace
332 const polyPatch&
patch,
333 const label patchToDualOffset,
337 const label startEdgeI,
348 DynamicList<label> dualFace(
pFaces.size());
350 DynamicList<label> featEdgeIndices(dualFace.size());
353 label edgeI = startEdgeI;
354 label facei =
patch.edgeFaces()[edgeI][0];
360 getPointEdges(
patch, facei, pointi, e0, e1);
374 dualFace.append(facei + patchToDualOffset);
378 getPointEdges(
patch, facei, pointi, e0, e1);
389 if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
392 dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
393 featEdgeIndices.append(dualFace.size()-1);
396 if (edgeI == startEdgeI)
404 if (eFaces[0] == facei)
414 dualFace2.transfer(dualFace);
416 featEdgeIndices2.transfer(featEdgeIndices);
423 forAll(featEdgeIndices2, i)
425 featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
433 void Foam::polyDualMesh::splitFace
435 const polyPatch&
patch,
442 DynamicList<face>& dualFaces,
443 DynamicList<label>& dualOwner,
444 DynamicList<label>& dualNeighbour,
445 DynamicList<label>& dualRegion
451 label meshPointi =
patch.meshPoints()[pointi];
453 if (pointToDualPoint[meshPointi] >= 0)
457 if (featEdgeIndices.size() < 2)
460 dualFaces.append(face(dualFace));
461 dualOwner.append(meshPointi);
462 dualNeighbour.append(-1);
463 dualRegion.append(
patch.index());
470 forAll(featEdgeIndices, i)
472 label startFp = featEdgeIndices[i];
474 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
481 sz = endFp - startFp + 2;
485 sz = endFp + dualFace.size() - startFp + 2;
490 subFace[0] = pointToDualPoint[
patch.meshPoints()[pointi]];
493 for (label subFp = 1; subFp < subFace.size(); subFp++)
495 subFace[subFp] = dualFace[startFp];
497 startFp = (startFp+1) % dualFace.size();
500 dualFaces.append(face(subFace));
501 dualOwner.append(meshPointi);
502 dualNeighbour.append(-1);
503 dualRegion.append(
patch.index());
510 if (featEdgeIndices.size() < 2)
513 dualFaces.append(face(dualFace));
514 dualOwner.append(meshPointi);
515 dualNeighbour.append(-1);
516 dualRegion.append(
patch.index());
525 DynamicList<label> subFace(dualFace.size());
527 forAll(featEdgeIndices, featI)
529 label startFp = featEdgeIndices[featI];
530 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
536 label vertI = dualFace[fp];
538 subFace.append(vertI);
545 fp = dualFace.fcIndex(fp);
548 if (subFace.size() > 2)
551 dualFaces.append(face(subFace));
552 dualOwner.append(meshPointi);
553 dualNeighbour.append(-1);
554 dualRegion.append(
patch.index());
560 if (subFace.size() > 2)
563 dualFaces.append(face(subFace));
564 dualOwner.append(meshPointi);
565 dualNeighbour.append(-1);
566 dualRegion.append(
patch.index());
576 void Foam::polyDualMesh::dualPatch
578 const polyPatch&
patch,
579 const label patchToDualOffset,
585 DynamicList<face>& dualFaces,
586 DynamicList<label>& dualOwner,
587 DynamicList<label>& dualNeighbour,
588 DynamicList<label>& dualRegion
600 bitSet donePoint(
patch.nPoints(),
false);
606 forAll(doneEdgeSide, patchEdgeI)
610 if (eFaces.size() == 1)
612 const edge&
e =
patch.edges()[patchEdgeI];
616 label bitMask = 1<<eI;
618 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
623 label pointi =
e[eI];
625 label edgeI = patchEdgeI;
641 if (
patch.edges()[edgeI][0] == pointi)
643 doneEdgeSide[edgeI] |= 1;
647 doneEdgeSide[edgeI] |= 2;
650 dualFaces.append(face(dualFace));
651 dualOwner.append(
patch.meshPoints()[pointi]);
652 dualNeighbour.append(-1);
653 dualRegion.append(
patch.index());
655 doneEdgeSide[patchEdgeI] |= bitMask;
656 donePoint.set(pointi);
669 if (!donePoint.test(pointi))
673 collectPatchInternalFace
679 patch.pointEdges()[pointi][0],
707 donePoint[pointi] =
true;
713 void Foam::polyDualMesh::calcDual
715 const polyMesh&
mesh,
739 allBoundary.meshEdges
747 pointSet nonManifoldPoints
751 meshEdges.
size() / 100
754 allBoundary.checkPointManifold(
true, &nonManifoldPoints);
756 if (nonManifoldPoints.size())
758 nonManifoldPoints.write();
761 <<
"There are " << nonManifoldPoints.size() <<
" points where"
762 <<
" the outside of the mesh is non-manifold." <<
nl
763 <<
"Such a mesh cannot be converted to a dual." <<
nl
764 <<
"Writing points at non-manifold sites to pointSet "
765 << nonManifoldPoints.name()
785 + featureEdges.size()
786 + featurePoints.size()
789 label dualPointi = 0;
795 cellPoint_.setSize(cellCentres.size());
797 forAll(cellCentres, celli)
799 cellPoint_[celli] = dualPointi;
800 dualPoints[dualPointi++] = cellCentres[celli];
806 boundaryFacePoint_.setSize(
mesh.
nFaces() - nIntFaces);
808 for (label facei = nIntFaces; facei <
mesh.
nFaces(); facei++)
810 boundaryFacePoint_[facei - nIntFaces] = dualPointi;
811 dualPoints[dualPointi++] = faceCentres[facei];
820 forAll(meshEdges, patchEdgeI)
822 label edgeI = meshEdges[patchEdgeI];
824 edgeToDualPoint[edgeI] = -1;
829 label edgeI = featureEdges[i];
833 edgeToDualPoint[edgeI] = dualPointi;
834 dualPoints[dualPointi++] =
e.centre(
mesh.
points());
852 pointToDualPoint[meshPoints[i]] = -2;
863 pointToDualPoint[meshPoints[loop[j]]] = -1;
870 label pointi = featurePoints[i];
872 pointToDualPoint[pointi] = dualPointi;
873 dualPoints[dualPointi++] =
mesh.
points()[pointi];
880 DynamicList<face> dynDualFaces(
mesh.
nEdges());
881 DynamicList<label> dynDualOwner(
mesh.
nEdges());
882 DynamicList<label> dynDualNeighbour(
mesh.
nEdges());
883 DynamicList<label> dynDualRegion(
mesh.
nEdges());
889 forAll(meshEdges, patchEdgeI)
891 label edgeI = meshEdges[patchEdgeI];
896 label neighbour = -1;
910 const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
912 if (patchFaces.size() != 2)
915 <<
"Cannot handle edges with " << patchFaces.size()
916 <<
" connected boundary faces."
920 label face0 = patchFaces[0] + nIntFaces;
923 label face1 = patchFaces[1] + nIntFaces;
928 label startFacei = -1;
931 label index = f0.find(neighbour);
933 if (f0.nextLabel(index) == owner)
946 DynamicList<label> dualFace;
948 if (edgeToDualPoint[edgeI] >= 0)
953 dualFace.
append(edgeToDualPoint[edgeI]);
961 dualFace.append(
mesh.
nCells() + startFacei - nIntFaces);
964 label facei = startFacei;
985 if (facei == endFacei)
1003 dynDualFaces.append(face(dualFace.shrink()));
1004 dynDualOwner.append(owner);
1005 dynDualNeighbour.append(neighbour);
1006 dynDualRegion.append(-1);
1010 const face&
f = dynDualFaces.last();
1011 const vector areaNorm =
f.areaNormal(dualPoints);
1013 if (((
mesh.
points()[owner] - dualPoints[
f[0]]) & areaNorm) > 0)
1016 <<
" on boundary edge:" << edgeI
1028 forAll(edgeToDualPoint, edgeI)
1030 if (edgeToDualPoint[edgeI] == -2)
1039 label neighbour = -1;
1055 label celli = eCells[0];
1065 label index = f0.find(neighbour);
1067 bool f0OrderOk = (f0.nextLabel(index) == owner);
1069 label startFacei = -1;
1081 DynamicList<label> dualFace(
mesh.
edgeCells()[edgeI].size());
1083 label facei = startFacei;
1104 if (facei == startFacei)
1119 dynDualFaces.
append(face(dualFace.shrink()));
1120 dynDualOwner.append(owner);
1121 dynDualNeighbour.append(neighbour);
1122 dynDualRegion.append(-1);
1126 const face&
f = dynDualFaces.last();
1127 const vector areaNorm =
f.areaNormal(dualPoints);
1129 if (((
mesh.
points()[owner] - dualPoints[
f[0]]) & areaNorm) > 0)
1132 <<
" on internal edge:" << edgeI
1144 dynDualFaces.shrink();
1145 dynDualOwner.shrink();
1146 dynDualNeighbour.shrink();
1147 dynDualRegion.shrink();
1149 OFstream str(
"dualInternalFaces.obj");
1150 Pout<<
"polyDualMesh::calcDual : dumping internal faces to "
1153 forAll(dualPoints, dualPointi)
1157 forAll(dynDualFaces, dualFacei)
1159 const face&
f = dynDualFaces[dualFacei];
1164 str<<
' ' <<
f[fp]+1;
1170 const label nInternalFaces = dynDualFaces.size();
1177 const polyPatch& pp =
patches[patchi];
1198 faceList dualFaces(std::move(dynDualFaces));
1199 labelList dualOwner(std::move(dynDualOwner));
1200 labelList dualNeighbour(std::move(dynDualNeighbour));
1201 labelList dualRegion(std::move(dynDualRegion));
1207 OFstream str(
"dualFaces.obj");
1208 Pout<<
"polyDualMesh::calcDual : dumping all faces to "
1211 forAll(dualPoints, dualPointi)
1215 forAll(dualFaces, dualFacei)
1217 const face&
f = dualFaces[dualFacei];
1222 str<<
' ' <<
f[fp]+1;
1233 dualCells[celli].setSize(0);
1238 label celli = dualOwner[facei];
1242 label sz = cFaces.size();
1246 forAll(dualNeighbour, facei)
1248 label celli = dualNeighbour[facei];
1254 label sz = cFaces.size();
1290 forAll(dualRegion, facei)
1292 if (dualRegion[facei] >= 0)
1294 patchSizes[dualRegion[facei]]++;
1300 label facei = nInternalFaces;
1304 patchStarts[patchi] = facei;
1306 facei += patchSizes[patchi];
1310 Pout<<
"nFaces:" << dualFaces.size()
1311 <<
" patchSizes:" << patchSizes
1312 <<
" patchStarts:" << patchStarts
1317 List<polyPatch*> dualPatches(
patches.size());
1321 const polyPatch& pp =
patches[patchi];
1323 dualPatches[patchi] = pp.
clone
1331 addPatches(dualPatches);
1357 time().findInstance(meshDir(),
"cellPoint"),
1368 "boundaryFacePoint",
1369 time().findInstance(meshDir(),
"boundaryFacePoint"),
1380 Foam::polyDualMesh::polyDualMesh
1393 time().findInstance(meshDir(),
"faces"),
1405 "boundaryFacePoint",
1406 time().findInstance(meshDir(),
"faces"),
1415 calcDual(
mesh, featureEdges, featurePoints);
1420 Foam::polyDualMesh::polyDualMesh
1423 const scalar featureCos
1432 time().findInstance(meshDir(),
"faces"),
1444 "boundaryFacePoint",
1445 time().findInstance(meshDir(),
"faces"),
1456 calcFeatures(
mesh, featureCos, featureEdges, featurePoints);
1457 calcDual(
mesh, featureEdges, featurePoints);
1464 const scalar featureCos,
1482 labelList allRegion(allBoundary.size());
1502 const labelList& meshPoints = allBoundary.meshPoints();
1504 bitSet isFeatureEdge(edgeFaces.size(),
false);
1508 const labelList& eFaces = edgeFaces[edgeI];
1510 if (eFaces.size() != 2)
1513 const edge&
e = allBoundary.edges()[edgeI];
1516 << meshPoints[
e[0]] <<
' ' << meshPoints[
e[1]]
1519 <<
" has more than 2 faces connected to it:"
1520 << eFaces.size() <<
endl;
1522 isFeatureEdge.set(edgeI);
1524 else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1526 isFeatureEdge.set(edgeI);
1534 isFeatureEdge.set(edgeI);
1546 forAll(pointEdges, pointi)
1548 const labelList& pEdges = pointEdges[pointi];
1550 label nFeatEdges = 0;
1554 if (isFeatureEdge.test(pEdges[i]))
1562 allFeaturePoints.append(allBoundary.meshPoints()[pointi]);
1565 featurePoints.
transfer(allFeaturePoints);
1571 Pout<<
"polyDualMesh::calcFeatures : dumping feature points to "
1576 label pointi = featurePoints[i];
1585 allBoundary.meshEdges
1599 forAll(isFeatureEdge, edgeI)
1601 if (isFeatureEdge.test(edgeI))
1604 allFeatureEdges.append(meshEdges[edgeI]);
1607 featureEdges.
transfer(allFeatureEdges);