65 SortableList<label> nbr(cFaces.size());
69 label faceI = cFaces[i];
71 label nbrCellI = faceNeighbour[faceI];
76 if (nbrCellI == cellI)
78 nbrCellI = faceOwner[faceI];
105 oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
110 nInternalFaces = newFaceI;
112 Pout<<
"nInternalFaces:" << nInternalFaces <<
endl;
113 Pout<<
"nFaces:" << faceOwner.size() <<
endl;
114 Pout<<
"nCells:" << cells.size() <<
endl;
117 for (label faceI = newFaceI; faceI < faceOwner.size(); faceI++)
119 oldToNew[faceI] = faceI;
126 if (oldToNew[faceI] == -1)
130 "polyDualMesh::getFaceOrder"
131 "(const labelList&, const labelList&, const label) const"
132 ) <<
"Did not determine new position"
133 <<
" for face " << faceI
147 void Foam::polyDualMesh::getPointEdges
156 const labelList& fEdges = patch.faceEdges()[faceI];
157 const face&
f = patch.localFaces()[faceI];
164 label edgeI = fEdges[i];
166 const edge&
e = patch.edges()[edgeI];
173 if (f.nextLabel(index) == e[1])
182 if (e0 != -1 && e1 != -1)
187 else if (e[1] == pointI)
192 if (f.nextLabel(index) == e[0])
201 if (e0 != -1 && e1 != -1)
208 FatalErrorIn(
"getPointEdges") <<
"Cannot find two edges on face:" << faceI
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)
250 FatalErrorIn(
"polyDualMesh::collectPatchSideFace") << edgeI
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]]);
297 const labelList& eFaces = patch.edgeFaces()[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,
344 const labelList& meshEdges = patch.meshEdges();
345 const labelList& pFaces = patch.pointFaces()[pointI];
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)
402 const labelList& eFaces = patch.edgeFaces()[edgeI];
404 if (eFaces[0] == faceI)
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)
553 dualFaces.append(face(subFace));
554 dualOwner.append(meshPointI);
555 dualNeighbour.append(-1);
556 dualRegion.append(patch.index());
562 if (subFace.size() > 2)
567 dualFaces.append(face(subFace));
568 dualOwner.append(meshPointI);
569 dualNeighbour.append(-1);
570 dualRegion.append(patch.index());
580 void Foam::polyDualMesh::dualPatch
582 const polyPatch& patch,
583 const label patchToDualOffset,
589 DynamicList<face>& dualFaces,
590 DynamicList<label>& dualOwner,
591 DynamicList<label>& dualNeighbour,
592 DynamicList<label>& dualRegion
595 const labelList& meshEdges = patch.meshEdges();
602 labelList doneEdgeSide(meshEdges.size(), 0);
604 boolList donePoint(patch.nPoints(),
false);
610 forAll(doneEdgeSide, patchEdgeI)
612 const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
614 if (eFaces.size() == 1)
616 const edge& e = patch.edges()[patchEdgeI];
620 label bitMask = 1<<eI;
622 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
627 label pointI = e[eI];
629 label edgeI = patchEdgeI;
645 if (patch.edges()[edgeI][0] == pointI)
647 doneEdgeSide[edgeI] |= 1;
651 doneEdgeSide[edgeI] |= 2;
654 dualFaces.append(face(dualFace));
655 dualOwner.append(patch.meshPoints()[pointI]);
656 dualNeighbour.append(-1);
657 dualRegion.append(patch.index());
659 doneEdgeSide[patchEdgeI] |= bitMask;
660 donePoint[pointI] =
true;
673 if (!donePoint[pointI])
677 collectPatchInternalFace
683 patch.pointEdges()[pointI][0],
711 donePoint[pointI] =
true;
717 void Foam::polyDualMesh::calcDual
719 const polyMesh&
mesh,
724 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
725 const label nIntFaces = mesh.nInternalFaces();
734 mesh.nFaces() - nIntFaces,
743 allBoundary.meshEdges
751 pointSet nonManifoldPoints
755 meshEdges.size() / 100
758 allBoundary.checkPointManifold(
true, &nonManifoldPoints);
760 if (nonManifoldPoints.size())
762 nonManifoldPoints.write();
766 "polyDualMesh::calcDual(const polyMesh&, const labelList&"
767 ", const labelList&)"
768 ) <<
"There are " << nonManifoldPoints.size() <<
" points where"
769 <<
" the outside of the mesh is non-manifold." <<
nl
770 <<
"Such a mesh cannot be converted to a dual." <<
nl
771 <<
"Writing points at non-manifold sites to pointSet "
772 << nonManifoldPoints.name()
791 + mesh.nFaces() - nIntFaces
792 + featureEdges.size()
793 + featurePoints.size()
796 label dualPointI = 0;
800 const pointField& cellCentres = mesh.cellCentres();
802 cellPoint_.
setSize(cellCentres.size());
804 forAll(cellCentres, cellI)
806 cellPoint_[cellI] = dualPointI;
807 dualPoints[dualPointI++] = cellCentres[cellI];
811 const pointField& faceCentres = mesh.faceCentres();
813 boundaryFacePoint_.
setSize(mesh.nFaces() - nIntFaces);
815 for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
817 boundaryFacePoint_[faceI - nIntFaces] = dualPointI;
818 dualPoints[dualPointI++] = faceCentres[faceI];
825 labelList edgeToDualPoint(mesh.nEdges(), -2);
827 forAll(meshEdges, patchEdgeI)
829 label edgeI = meshEdges[patchEdgeI];
831 edgeToDualPoint[edgeI] = -1;
836 label edgeI = featureEdges[i];
838 const edge& e = mesh.edges()[edgeI];
840 edgeToDualPoint[edgeI] = dualPointI;
841 dualPoints[dualPointI++] = e.centre(mesh.points());
851 labelList pointToDualPoint(mesh.nPoints(), -3);
855 const labelList& meshPoints = patches[patchI].meshPoints();
859 pointToDualPoint[meshPoints[i]] = -2;
870 pointToDualPoint[meshPoints[loop[j]]] = -1;
877 label pointI = featurePoints[i];
879 pointToDualPoint[pointI] = dualPointI;
880 dualPoints[dualPointI++] = mesh.points()[pointI];
887 DynamicList<face> dynDualFaces(mesh.nEdges());
888 DynamicList<label> dynDualOwner(mesh.nEdges());
889 DynamicList<label> dynDualNeighbour(mesh.nEdges());
890 DynamicList<label> dynDualRegion(mesh.nEdges());
896 forAll(meshEdges, patchEdgeI)
898 label edgeI = meshEdges[patchEdgeI];
900 const edge& e = mesh.edges()[edgeI];
903 label neighbour = -1;
917 const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
919 if (patchFaces.size() != 2)
922 <<
"Cannot handle edges with " << patchFaces.size()
923 <<
" connected boundary faces."
927 label face0 = patchFaces[0] + nIntFaces;
928 const face& f0 = mesh.faces()[face0];
930 label face1 = patchFaces[1] + nIntFaces;
935 label startFaceI = -1;
940 if (f0.nextLabel(index) == owner)
953 DynamicList<label> dualFace;
955 if (edgeToDualPoint[edgeI] >= 0)
958 dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2+1);
960 dualFace.append(edgeToDualPoint[edgeI]);
964 dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2);
968 dualFace.append(mesh.nCells() + startFaceI - nIntFaces);
970 label cellI = mesh.faceOwner()[startFaceI];
971 label faceI = startFaceI;
976 dualFace.append(cellI);
992 if (faceI == endFaceI)
997 if (mesh.faceOwner()[faceI] == cellI)
999 cellI = mesh.faceNeighbour()[faceI];
1003 cellI = mesh.faceOwner()[faceI];
1008 dualFace.append(mesh.nCells() + endFaceI - nIntFaces);
1010 dynDualFaces.append(face(dualFace.shrink()));
1011 dynDualOwner.append(owner);
1012 dynDualNeighbour.append(neighbour);
1013 dynDualRegion.append(-1);
1017 const face& f = dynDualFaces[dynDualFaces.size()-1];
1018 vector n = f.normal(dualPoints);
1019 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1021 WarningIn(
"calcDual") <<
"Incorrect orientation"
1022 <<
" on boundary edge:" << edgeI
1023 << mesh.points()[mesh.edges()[edgeI][0]]
1024 << mesh.points()[mesh.edges()[edgeI][1]]
1034 forAll(edgeToDualPoint, edgeI)
1036 if (edgeToDualPoint[edgeI] == -2)
1042 const edge& e = mesh.edges()[edgeI];
1045 label neighbour = -1;
1059 const labelList& eCells = mesh.edgeCells()[edgeI];
1061 label cellI = eCells[0];
1069 const face& f0 = mesh.faces()[face0];
1073 bool f0OrderOk = (f0.nextLabel(index) == owner);
1075 label startFaceI = -1;
1077 if (f0OrderOk == (mesh.faceOwner()[face0] == cellI))
1087 DynamicList<label> dualFace(mesh.edgeCells()[edgeI].size());
1089 label faceI = startFaceI;
1094 dualFace.append(cellI);
1110 if (faceI == startFaceI)
1115 if (mesh.faceOwner()[faceI] == cellI)
1117 cellI = mesh.faceNeighbour()[faceI];
1121 cellI = mesh.faceOwner()[faceI];
1125 dynDualFaces.append(face(dualFace.shrink()));
1126 dynDualOwner.append(owner);
1127 dynDualNeighbour.append(neighbour);
1128 dynDualRegion.append(-1);
1132 const face& f = dynDualFaces[dynDualFaces.size()-1];
1133 vector n = f.normal(dualPoints);
1134 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1136 WarningIn(
"calcDual") <<
"Incorrect orientation"
1137 <<
" on internal edge:" << edgeI
1138 << mesh.points()[mesh.edges()[edgeI][0]]
1139 << mesh.points()[mesh.edges()[edgeI][1]]
1149 dynDualFaces.shrink();
1150 dynDualOwner.shrink();
1151 dynDualNeighbour.shrink();
1152 dynDualRegion.shrink();
1154 OFstream str(
"dualInternalFaces.obj");
1155 Pout<<
"polyDualMesh::calcDual : dumping internal faces to "
1158 forAll(dualPoints, dualPointI)
1162 forAll(dynDualFaces, dualFaceI)
1164 const face& f = dynDualFaces[dualFaceI];
1169 str<<
' ' << f[fp]+1;
1175 const label nInternalFaces = dynDualFaces.size();
1182 const polyPatch& pp = patches[patchI];
1187 mesh.nCells() + pp.start() - nIntFaces,
1203 faceList dualFaces(dynDualFaces.shrink(),
true);
1204 dynDualFaces.
clear();
1206 labelList dualOwner(dynDualOwner.shrink(),
true);
1207 dynDualOwner.
clear();
1209 labelList dualNeighbour(dynDualNeighbour.shrink(),
true);
1210 dynDualNeighbour.
clear();
1212 labelList dualRegion(dynDualRegion.shrink(),
true);
1213 dynDualRegion.
clear();
1220 OFstream str(
"dualFaces.obj");
1221 Pout<<
"polyDualMesh::calcDual : dumping all faces to "
1224 forAll(dualPoints, dualPointI)
1228 forAll(dualFaces, dualFaceI)
1230 const face& f = dualFaces[dualFaceI];
1235 str<<
' ' << f[fp]+1;
1242 cellList dualCells(mesh.nPoints());
1246 dualCells[cellI].setSize(0);
1251 label cellI = dualOwner[faceI];
1255 label sz = cFaces.
size();
1256 cFaces.setSize(sz+1);
1259 forAll(dualNeighbour, faceI)
1261 label cellI = dualNeighbour[faceI];
1267 label sz = cFaces.
size();
1268 cFaces.setSize(sz+1);
1301 labelList patchSizes(patches.size(), 0);
1303 forAll(dualRegion, faceI)
1305 if (dualRegion[faceI] >= 0)
1307 patchSizes[dualRegion[faceI]]++;
1311 labelList patchStarts(patches.size(), 0);
1313 label faceI = nInternalFaces;
1317 patchStarts[patchI] = faceI;
1319 faceI += patchSizes[patchI];
1323 Pout<<
"nFaces:" << dualFaces.size()
1324 <<
" patchSizes:" << patchSizes
1325 <<
" patchStarts:" << patchStarts
1330 List<polyPatch*> dualPatches(patches.size());
1334 const polyPatch& pp = patches[patchI];
1336 dualPatches[patchI] = pp.clone
1344 addPatches(dualPatches);
1370 time().findInstance(meshDir(),
"cellPoint"),
1381 "boundaryFacePoint",
1382 time().findInstance(meshDir(),
"boundaryFacePoint"),
1393 Foam::polyDualMesh::polyDualMesh
1412 time().findInstance(meshDir(),
"faces"),
1424 "boundaryFacePoint",
1425 time().findInstance(meshDir(),
"faces"),
1434 calcDual(mesh, featureEdges, featurePoints);
1439 Foam::polyDualMesh::polyDualMesh
1442 const scalar featureCos
1457 time().findInstance(meshDir(),
"faces"),
1469 "boundaryFacePoint",
1470 time().findInstance(meshDir(),
"faces"),
1481 calcFeatures(mesh, featureCos, featureEdges, featurePoints);
1482 calcDual(mesh, featureEdges, featurePoints);
1489 const scalar featureCos,
1507 labelList allRegion(allBoundary.size());
1526 const vectorField& faceNormals = allBoundary.faceNormals();
1527 const labelList& meshPoints = allBoundary.meshPoints();
1533 const labelList& eFaces = edgeFaces[edgeI];
1535 if (eFaces.
size() != 2)
1538 const edge& e = allBoundary.edges()[edgeI];
1540 WarningIn(
"polyDualMesh::calcFeatures") <<
"Edge "
1541 << meshPoints[e[0]] <<
' ' << meshPoints[e[1]]
1542 <<
" coords:" << mesh.
points()[meshPoints[e[0]]]
1543 << mesh.
points()[meshPoints[e[1]]]
1544 <<
" has more than 2 faces connected to it:"
1547 isFeatureEdge[edgeI] =
true;
1549 else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1551 isFeatureEdge[edgeI] =
true;
1555 (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1559 isFeatureEdge[edgeI] =
true;
1571 forAll(pointEdges, pointI)
1573 const labelList& pEdges = pointEdges[pointI];
1575 label nFeatEdges = 0;
1579 if (isFeatureEdge[pEdges[i]])
1587 allFeaturePoints.append(allBoundary.meshPoints()[pointI]);
1590 featurePoints.
transfer(allFeaturePoints);
1596 Pout<<
"polyDualMesh::calcFeatures : dumping feature points to "
1601 label pointI = featurePoints[i];
1610 allBoundary.meshEdges
1624 forAll(isFeatureEdge, edgeI)
1626 if (isFeatureEdge[edgeI])
1629 allFeatureEdges.append(meshEdges[edgeI]);
1632 featureEdges.
transfer(allFeatureEdges);