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)
129 <<
"Did not determine new position" 130 <<
" for face " << facei
144 void Foam::polyDualMesh::getPointEdges
153 const labelList& fEdges = patch.faceEdges()[facei];
154 const face& f = patch.localFaces()[facei];
161 label edgeI = fEdges[i];
163 const edge& e = patch.edges()[edgeI];
170 if (f.nextLabel(index) == e[1])
179 if (e0 != -1 && e1 != -1)
184 else if (e[1] == pointi)
189 if (f.nextLabel(index) == e[0])
198 if (e0 != -1 && e1 != -1)
206 <<
" vertices:" << patch.localFaces()[facei]
207 <<
" that uses point:" << pointi
215 const polyPatch& patch,
216 const label patchToDualOffset,
227 label meshPointi = patch.meshPoints()[pointi];
228 const labelList& pFaces = patch.pointFaces()[pointi];
230 DynamicList<label> dualFace;
232 if (pointToDualPoint[meshPointi] >= 0)
235 dualFace.setCapacity(pFaces.size()+2+1);
237 dualFace.
append(pointToDualPoint[meshPointi]);
241 dualFace.setCapacity(pFaces.size()+2);
245 if (edgeToDualPoint[patch.meshEdges()[edgeI]] < 0)
251 dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
253 label facei = patch.edgeFaces()[edgeI][0];
259 getPointEdges(patch, facei, pointi, e0, e1);
273 dualFace.append(facei + patchToDualOffset);
277 getPointEdges(patch, facei, pointi, e0, e1);
288 if (edgeToDualPoint[patch.meshEdges()[edgeI]] >= 0)
291 dualFace.append(edgeToDualPoint[patch.meshEdges()[edgeI]]);
294 const labelList& eFaces = patch.edgeFaces()[edgeI];
296 if (eFaces.size() == 1)
303 if (eFaces[0] == facei)
327 void Foam::polyDualMesh::collectPatchInternalFace
329 const polyPatch& patch,
330 const label patchToDualOffset,
334 const label startEdgeI,
341 const labelList& meshEdges = patch.meshEdges();
342 const labelList& pFaces = patch.pointFaces()[pointi];
345 DynamicList<label> dualFace(pFaces.size());
347 DynamicList<label> featEdgeIndices(dualFace.size());
350 label edgeI = startEdgeI;
351 label facei = patch.edgeFaces()[edgeI][0];
357 getPointEdges(patch, facei, pointi, e0, e1);
371 dualFace.append(facei + patchToDualOffset);
375 getPointEdges(patch, facei, pointi, e0, e1);
386 if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
389 dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
390 featEdgeIndices.append(dualFace.size()-1);
393 if (edgeI == startEdgeI)
399 const labelList& eFaces = patch.edgeFaces()[edgeI];
401 if (eFaces[0] == facei)
413 featEdgeIndices2.transfer(featEdgeIndices);
420 forAll(featEdgeIndices2, i)
422 featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
430 void Foam::polyDualMesh::splitFace
432 const polyPatch& patch,
439 DynamicList<face>& dualFaces,
440 DynamicList<label>& dualOwner,
441 DynamicList<label>& dualNeighbour,
442 DynamicList<label>& dualRegion
448 label meshPointi = patch.meshPoints()[pointi];
450 if (pointToDualPoint[meshPointi] >= 0)
454 if (featEdgeIndices.size() < 2)
457 dualFaces.append(face(dualFace));
458 dualOwner.append(meshPointi);
459 dualNeighbour.append(-1);
460 dualRegion.append(patch.index());
467 forAll(featEdgeIndices, i)
469 label startFp = featEdgeIndices[i];
471 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
478 sz = endFp - startFp + 2;
482 sz = endFp + dualFace.size() - startFp + 2;
487 subFace[0] = pointToDualPoint[patch.meshPoints()[pointi]];
490 for (
label subFp = 1; subFp < subFace.size(); subFp++)
492 subFace[subFp] = dualFace[startFp];
494 startFp = (startFp+1) % dualFace.size();
497 dualFaces.append(face(subFace));
498 dualOwner.append(meshPointi);
499 dualNeighbour.append(-1);
500 dualRegion.append(patch.index());
507 if (featEdgeIndices.size() < 2)
510 dualFaces.append(face(dualFace));
511 dualOwner.append(meshPointi);
512 dualNeighbour.append(-1);
513 dualRegion.append(patch.index());
522 DynamicList<label> subFace(dualFace.size());
524 forAll(featEdgeIndices, featI)
526 label startFp = featEdgeIndices[featI];
527 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
533 label vertI = dualFace[fp];
535 subFace.append(vertI);
542 fp = dualFace.fcIndex(fp);
545 if (subFace.size() > 2)
550 dualFaces.append(face(subFace));
551 dualOwner.append(meshPointi);
552 dualNeighbour.append(-1);
553 dualRegion.append(patch.index());
559 if (subFace.size() > 2)
564 dualFaces.append(face(subFace));
565 dualOwner.append(meshPointi);
566 dualNeighbour.append(-1);
567 dualRegion.append(patch.index());
577 void Foam::polyDualMesh::dualPatch
579 const polyPatch& patch,
580 const label patchToDualOffset,
586 DynamicList<face>& dualFaces,
587 DynamicList<label>& dualOwner,
588 DynamicList<label>& dualNeighbour,
589 DynamicList<label>& dualRegion
592 const labelList& meshEdges = patch.meshEdges();
599 labelList doneEdgeSide(meshEdges.size(), 0);
601 boolList donePoint(patch.nPoints(),
false);
607 forAll(doneEdgeSide, patchEdgeI)
609 const labelList& eFaces = patch.edgeFaces()[patchEdgeI];
611 if (eFaces.size() == 1)
613 const edge& e = patch.edges()[patchEdgeI];
617 label bitMask = 1<<eI;
619 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
624 label pointi = e[eI];
626 label edgeI = patchEdgeI;
642 if (patch.edges()[edgeI][0] == pointi)
644 doneEdgeSide[edgeI] |= 1;
648 doneEdgeSide[edgeI] |= 2;
651 dualFaces.append(face(dualFace));
652 dualOwner.append(patch.meshPoints()[pointi]);
653 dualNeighbour.append(-1);
654 dualRegion.append(patch.index());
656 doneEdgeSide[patchEdgeI] |= bitMask;
657 donePoint[pointi] =
true;
670 if (!donePoint[pointi])
674 collectPatchInternalFace
680 patch.pointEdges()[pointi][0],
708 donePoint[pointi] =
true;
714 void Foam::polyDualMesh::calcDual
716 const polyMesh& mesh,
721 const polyBoundaryMesh& patches = mesh.boundaryMesh();
722 const label nIntFaces = mesh.nInternalFaces();
731 mesh.nFaces() - nIntFaces,
740 allBoundary.meshEdges
748 pointSet nonManifoldPoints
752 meshEdges.size() / 100
755 allBoundary.checkPointManifold(
true, &nonManifoldPoints);
757 if (nonManifoldPoints.size())
759 nonManifoldPoints.write();
762 <<
"There are " << nonManifoldPoints.size() <<
" points where" 763 <<
" the outside of the mesh is non-manifold." <<
nl 764 <<
"Such a mesh cannot be converted to a dual." <<
nl 765 <<
"Writing points at non-manifold sites to pointSet " 766 << nonManifoldPoints.name()
785 + mesh.nFaces() - nIntFaces
786 + featureEdges.size()
787 + featurePoints.size()
790 label dualPointi = 0;
794 const pointField& cellCentres = mesh.cellCentres();
796 cellPoint_.
setSize(cellCentres.size());
798 forAll(cellCentres, celli)
800 cellPoint_[celli] = dualPointi;
801 dualPoints[dualPointi++] = cellCentres[celli];
805 const pointField& faceCentres = mesh.faceCentres();
807 boundaryFacePoint_.
setSize(mesh.nFaces() - nIntFaces);
809 for (
label facei = nIntFaces; facei < mesh.nFaces(); facei++)
811 boundaryFacePoint_[facei - nIntFaces] = dualPointi;
812 dualPoints[dualPointi++] = faceCentres[facei];
819 labelList edgeToDualPoint(mesh.nEdges(), -2);
821 forAll(meshEdges, patchEdgeI)
823 label edgeI = meshEdges[patchEdgeI];
825 edgeToDualPoint[edgeI] = -1;
830 label edgeI = featureEdges[i];
832 const edge& e = mesh.edges()[edgeI];
834 edgeToDualPoint[edgeI] = dualPointi;
835 dualPoints[dualPointi++] = e.centre(mesh.points());
845 labelList pointToDualPoint(mesh.nPoints(), -3);
853 pointToDualPoint[meshPoints[i]] = -2;
864 pointToDualPoint[meshPoints[loop[j]]] = -1;
871 label pointi = featurePoints[i];
873 pointToDualPoint[pointi] = dualPointi;
874 dualPoints[dualPointi++] = mesh.points()[pointi];
881 DynamicList<face> dynDualFaces(mesh.nEdges());
882 DynamicList<label> dynDualOwner(mesh.nEdges());
883 DynamicList<label> dynDualNeighbour(mesh.nEdges());
884 DynamicList<label> dynDualRegion(mesh.nEdges());
890 forAll(meshEdges, patchEdgeI)
892 label edgeI = meshEdges[patchEdgeI];
894 const edge& e = mesh.edges()[edgeI];
897 label neighbour = -1;
911 const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
913 if (patchFaces.size() != 2)
916 <<
"Cannot handle edges with " << patchFaces.size()
917 <<
" connected boundary faces." 921 label face0 = patchFaces[0] + nIntFaces;
922 const face& f0 = mesh.faces()[face0];
924 label face1 = patchFaces[1] + nIntFaces;
929 label startFacei = -1;
934 if (f0.nextLabel(index) == owner)
947 DynamicList<label> dualFace;
949 if (edgeToDualPoint[edgeI] >= 0)
952 dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2+1);
954 dualFace.append(edgeToDualPoint[edgeI]);
958 dualFace.setCapacity(mesh.edgeCells()[edgeI].size()+2);
962 dualFace.append(mesh.nCells() + startFacei - nIntFaces);
964 label celli = mesh.faceOwner()[startFacei];
965 label facei = startFacei;
970 dualFace.append(celli);
986 if (facei == endFacei)
991 if (mesh.faceOwner()[facei] == celli)
993 celli = mesh.faceNeighbour()[facei];
997 celli = mesh.faceOwner()[facei];
1002 dualFace.append(mesh.nCells() + endFacei - nIntFaces);
1004 dynDualFaces.append(face(dualFace.shrink()));
1005 dynDualOwner.append(owner);
1006 dynDualNeighbour.append(neighbour);
1007 dynDualRegion.append(-1);
1011 const face& f = dynDualFaces.last();
1012 vector n = f.normal(dualPoints);
1013 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1016 <<
" on boundary edge:" << edgeI
1017 << mesh.points()[mesh.edges()[edgeI][0]]
1018 << mesh.points()[mesh.edges()[edgeI][1]]
1028 forAll(edgeToDualPoint, edgeI)
1030 if (edgeToDualPoint[edgeI] == -2)
1036 const edge& e = mesh.edges()[edgeI];
1039 label neighbour = -1;
1053 const labelList& eCells = mesh.edgeCells()[edgeI];
1055 label celli = eCells[0];
1063 const face& f0 = mesh.faces()[face0];
1067 bool f0OrderOk = (f0.nextLabel(index) == owner);
1069 label startFacei = -1;
1071 if (f0OrderOk == (mesh.faceOwner()[face0] == celli))
1081 DynamicList<label> dualFace(mesh.edgeCells()[edgeI].size());
1083 label facei = startFacei;
1088 dualFace.append(celli);
1104 if (facei == startFacei)
1109 if (mesh.faceOwner()[facei] == celli)
1111 celli = mesh.faceNeighbour()[facei];
1115 celli = mesh.faceOwner()[facei];
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 vector n = f.normal(dualPoints);
1128 if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
1131 <<
" on internal edge:" << edgeI
1132 << mesh.points()[mesh.edges()[edgeI][0]]
1133 << mesh.points()[mesh.edges()[edgeI][1]]
1143 dynDualFaces.shrink();
1144 dynDualOwner.shrink();
1145 dynDualNeighbour.shrink();
1146 dynDualRegion.shrink();
1148 OFstream str(
"dualInternalFaces.obj");
1149 Pout<<
"polyDualMesh::calcDual : dumping internal faces to " 1152 forAll(dualPoints, dualPointi)
1156 forAll(dynDualFaces, dualFacei)
1158 const face& f = dynDualFaces[dualFacei];
1163 str<<
' ' << f[fp]+1;
1169 const label nInternalFaces = dynDualFaces.size();
1176 const polyPatch& pp = patches[
patchi];
1181 mesh.nCells() + pp.start() - nIntFaces,
1197 faceList dualFaces(dynDualFaces.shrink(),
true);
1198 dynDualFaces.
clear();
1200 labelList dualOwner(dynDualOwner.shrink(),
true);
1201 dynDualOwner.
clear();
1203 labelList dualNeighbour(dynDualNeighbour.shrink(),
true);
1204 dynDualNeighbour.
clear();
1206 labelList dualRegion(dynDualRegion.shrink(),
true);
1207 dynDualRegion.
clear();
1214 OFstream str(
"dualFaces.obj");
1215 Pout<<
"polyDualMesh::calcDual : dumping all faces to " 1218 forAll(dualPoints, dualPointi)
1222 forAll(dualFaces, dualFacei)
1224 const face& f = dualFaces[dualFacei];
1229 str<<
' ' << f[fp]+1;
1236 cellList dualCells(mesh.nPoints());
1240 dualCells[celli].setSize(0);
1245 label celli = dualOwner[facei];
1249 label sz = cFaces.size();
1250 cFaces.setSize(sz+1);
1253 forAll(dualNeighbour, facei)
1255 label celli = dualNeighbour[facei];
1261 label sz = cFaces.size();
1262 cFaces.setSize(sz+1);
1295 labelList patchSizes(patches.size(), 0);
1297 forAll(dualRegion, facei)
1299 if (dualRegion[facei] >= 0)
1301 patchSizes[dualRegion[facei]]++;
1305 labelList patchStarts(patches.size(), 0);
1307 label facei = nInternalFaces;
1311 patchStarts[
patchi] = facei;
1313 facei += patchSizes[
patchi];
1317 Pout<<
"nFaces:" << dualFaces.size()
1318 <<
" patchSizes:" << patchSizes
1319 <<
" patchStarts:" << patchStarts
1324 List<polyPatch*> dualPatches(patches.size());
1328 const polyPatch& pp = patches[
patchi];
1330 dualPatches[
patchi] = pp.clone
1338 addPatches(dualPatches);
1364 time().findInstance(meshDir(),
"cellPoint"),
1375 "boundaryFacePoint",
1376 time().findInstance(meshDir(),
"boundaryFacePoint"),
1387 Foam::polyDualMesh::polyDualMesh
1418 "boundaryFacePoint",
1428 calcDual(mesh, featureEdges, featurePoints);
1433 Foam::polyDualMesh::polyDualMesh
1436 const scalar featureCos
1463 "boundaryFacePoint",
1475 calcFeatures(mesh, featureCos, featureEdges, featurePoints);
1476 calcDual(mesh, featureEdges, featurePoints);
1483 const scalar featureCos,
1501 labelList allRegion(allBoundary.size());
1520 const vectorField& faceNormals = allBoundary.faceNormals();
1521 const labelList& meshPoints = allBoundary.meshPoints();
1527 const labelList& eFaces = edgeFaces[edgeI];
1529 if (eFaces.
size() != 2)
1532 const edge& e = allBoundary.edges()[edgeI];
1535 << meshPoints[e[0]] <<
' ' << meshPoints[e[1]]
1536 <<
" coords:" << mesh.
points()[meshPoints[e[0]]]
1537 << mesh.
points()[meshPoints[e[1]]]
1538 <<
" has more than 2 faces connected to it:" 1541 isFeatureEdge[edgeI] =
true;
1543 else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1545 isFeatureEdge[edgeI] =
true;
1549 (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1553 isFeatureEdge[edgeI] =
true;
1565 forAll(pointEdges, pointi)
1567 const labelList& pEdges = pointEdges[pointi];
1569 label nFeatEdges = 0;
1573 if (isFeatureEdge[pEdges[i]])
1581 allFeaturePoints.append(allBoundary.meshPoints()[pointi]);
1584 featurePoints.
transfer(allFeaturePoints);
1590 Pout<<
"polyDualMesh::calcFeatures : dumping feature points to " 1595 label pointi = featurePoints[i];
1604 allBoundary.meshEdges
1618 forAll(isFeatureEdge, edgeI)
1620 if (isFeatureEdge[edgeI])
1623 allFeatureEdges.append(meshEdges[edgeI]);
1626 featureEdges.
transfer(allFeatureEdges);
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const fileName & name() const
Return the name of the stream.
List< labelList > labelListList
A List of labelList.
const labelListList & cellEdges() const
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 inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
errorManipArg< error, int > exit(error &err, const int errNo=1)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label nInternalFaces() const
const labelListList & pointEdges() const
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
void size(const label)
Override size to be inconsistent with allocated storage.
static void calcFeatures(const polyMesh &, const scalar featureCos, labelList &featureEdges, labelList &featurePoints)
Helper function to create feature edges and points based on.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Xfer< T > xferCopy(const T &)
Construct by copying the contents of the arg.
Vector< scalar > vector
A scalar version of the templated Vector.
const fileName & name() const
Return the name of the stream.
virtual const pointField & points() const
Return raw points.
List< bool > boolList
Bool container classes.
A list of faces which address into the list of points.
A List obtained as a section of another List.
vectorField pointField
pointField is a vectorField.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
void clear()
Clear the list, i.e. set size to zero.
void append(const T &)
Append an element at the end of the list.
virtual const labelList & faceOwner() const
Return face owner.
List< label > labelList
A List of labels.
virtual const faceList & faces() const
Return raw faces.
errorManip< error > abort(error &err)
void reverse(UList< T > &, const label n)
const Time & time() const
Return time.
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
~polyDualMesh()
Destructor.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void setSize(const label)
Reset size of List.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
fileName meshDir() const
Return the local mesh directory (dbDir()/meshSubDir)
label start() const
Return start label of this patch in the polyMesh face list.
Mesh consisting of general polyhedral cells.
A patch is a list of labels that address the faces in the global face list.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
List< cell > cellList
list of cells
const labelListList & edgeFaces() const