42 void Foam::meshDualiser::checkPolyTopoChange(
const polyTopoChange& meshMod)
46 forAll(meshMod.points(), i)
48 points[i] = meshMod.points()[i];
68 if (newToOld[newI].size() != 1)
72 "meshDualiser::checkPolyTopoChange(const polyTopoChange&)"
73 ) <<
"duplicate verts:" << newToOld[newI]
75 << UIndirectList<point>(
points, newToOld[newI])()
84 void Foam::meshDualiser::dumpPolyTopoChange
86 const polyTopoChange& meshMod,
87 const fileName& prefix
90 OFstream str1(prefix +
"Faces.obj");
91 OFstream str2(prefix +
"Edges.obj");
93 Info<<
"Dumping current polyTopoChange. Faces to " << str1.name()
94 <<
" , points and edges to " << str2.name() <<
endl;
96 const DynamicList<point>&
points = meshMod.points();
104 const DynamicList<face>& faces = meshMod.faces();
108 const face&
f = faces[faceI];
113 str1<<
' ' << f[fp]+1;
120 str2<<
' ' << f[fp]+1;
122 str2<<
' ' << f[0]+1 <<
nl;
129 Foam::label Foam::meshDualiser::findDualCell
135 const labelList& dualCells = pointToDualCells_[pointI];
137 if (dualCells.size() == 1)
143 label index =
findIndex(mesh_.pointCells()[pointI], cellI);
145 return dualCells[index];
152 void Foam::meshDualiser::generateDualBoundaryEdges
156 polyTopoChange& meshMod
159 const labelList& pEdges = mesh_.pointEdges()[pointI];
163 label edgeI = pEdges[pEdgeI];
165 if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.get(edgeI) == 1)
167 const edge&
e = mesh_.edges()[edgeI];
169 edgeToDualPoint_[edgeI] = meshMod.addPoint
171 e.centre(mesh_.points()),
183 bool Foam::meshDualiser::sameDualCell
189 if (!mesh_.isInternalFace(faceI))
192 <<
"face:" << faceI <<
" is not internal face."
196 label own = mesh_.faceOwner()[faceI];
197 label nei = mesh_.faceNeighbour()[faceI];
199 return findDualCell(own, pointI) == findDualCell(nei, pointI);
203 Foam::label Foam::meshDualiser::addInternalFace
205 const label masterPointI,
206 const label masterEdgeI,
207 const label masterFaceI,
209 const bool edgeOrder,
210 const label dualCell0,
211 const label dualCell1,
212 const DynamicList<label>& verts,
213 polyTopoChange& meshMod
218 if (edgeOrder != (dualCell0 < dualCell1))
225 pointField facePoints(meshMod.points(), newFace);
241 <<
"verts:" << verts <<
" newFace:" << newFace
242 <<
" face points:" << facePoints
249 bool zoneFlip =
false;
250 if (masterFaceI != -1)
252 zoneID = mesh_.faceZones().whichZone(masterFaceI);
256 const faceZone& fZone = mesh_.faceZones()[zoneID];
258 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
264 if (dualCell0 < dualCell1)
266 dualFaceI = meshMod.addFace
293 dualFaceI = meshMod.addFace
322 Foam::label Foam::meshDualiser::addBoundaryFace
324 const label masterPointI,
325 const label masterEdgeI,
326 const label masterFaceI,
328 const label dualCellI,
330 const DynamicList<label>& verts,
331 polyTopoChange& meshMod
337 bool zoneFlip =
false;
338 if (masterFaceI != -1)
340 zoneID = mesh_.faceZones().whichZone(masterFaceI);
344 const faceZone& fZone = mesh_.faceZones()[zoneID];
346 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFaceI)];
350 label dualFaceI = meshMod.addFace
381 void Foam::meshDualiser::createFacesAroundEdge
383 const bool splitFace,
386 const label startFaceI,
387 polyTopoChange& meshMod,
391 const edge&
e = mesh_.edges()[edgeI];
392 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
396 mesh_.faces()[startFaceI],
401 edgeFaceCirculator ie
407 isBoundaryEdge.get(edgeI) == 1
411 bool edgeOrder = ie.sameOrder(e[0], e[1]);
412 label startFaceLabel = ie.faceLabel();
423 DynamicList<label> verts(100);
425 if (edgeToDualPoint_[edgeI] != -1)
427 verts.append(edgeToDualPoint_[edgeI]);
429 if (faceToDualPoint_[ie.faceLabel()] != -1)
431 doneEFaces[
findIndex(eFaces, ie.faceLabel())] =
true;
432 verts.append(faceToDualPoint_[ie.faceLabel()]);
434 if (cellToDualPoint_[ie.cellLabel()] != -1)
436 verts.append(cellToDualPoint_[ie.cellLabel()]);
439 label currentDualCell0 = findDualCell(ie.cellLabel(), e[0]);
440 label currentDualCell1 = findDualCell(ie.cellLabel(), e[1]);
446 label faceI = ie.faceLabel();
449 doneEFaces[
findIndex(eFaces, faceI)] =
true;
451 if (faceToDualPoint_[faceI] != -1)
453 verts.append(faceToDualPoint_[faceI]);
456 label cellI = ie.cellLabel();
466 label dualCell0 = findDualCell(cellI, e[0]);
467 label dualCell1 = findDualCell(cellI, e[1]);
475 dualCell0 != currentDualCell0
476 || dualCell1 != currentDualCell1
494 currentDualCell0 = dualCell0;
495 currentDualCell1 = dualCell1;
498 if (edgeToDualPoint_[edgeI] != -1)
500 verts.append(edgeToDualPoint_[edgeI]);
502 if (faceToDualPoint_[faceI] != -1)
504 verts.append(faceToDualPoint_[faceI]);
508 if (cellToDualPoint_[cellI] != -1)
510 verts.append(cellToDualPoint_[cellI]);
519 if (isBoundaryEdge.get(edgeI) == 0)
521 label startDual = faceToDualPoint_[startFaceLabel];
523 if (startDual != -1 &&
findIndex(verts, startDual) == -1)
525 verts.append(startDual);
550 void Foam::meshDualiser::createFaceFromInternalFace
554 polyTopoChange& meshMod
557 const face&
f = mesh_.faces()[faceI];
558 const labelList& fEdges = mesh_.faceEdges()[faceI];
559 label own = mesh_.faceOwner()[faceI];
560 label nei = mesh_.faceNeighbour()[faceI];
571 DynamicList<label> verts(100);
573 verts.append(faceToDualPoint_[faceI]);
574 verts.append(edgeToDualPoint_[fEdges[fp]]);
580 label currentDualCell0 = findDualCell(own, f[fp]);
581 label currentDualCell1 = findDualCell(nei, f[fp]);
586 if (pointToDualPoint_[f[fp]] != -1)
588 verts.append(pointToDualPoint_[f[fp]]);
592 label edgeI = fEdges[fp];
594 if (edgeToDualPoint_[edgeI] != -1)
596 verts.append(edgeToDualPoint_[edgeI]);
600 label nextFp = f.fcIndex(fp);
603 label dualCell0 = findDualCell(own, f[nextFp]);
604 label dualCell1 = findDualCell(nei, f[nextFp]);
606 if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
609 if (edgeToDualPoint_[edgeI] == -1)
612 <<
"face:" << faceI <<
" verts:" << f
613 <<
" points:" << UIndirectList<point>(mesh_.points(),
f)()
614 <<
" no feature edge between " << f[fp]
615 <<
" and " << f[nextFp] <<
" although have different"
616 <<
" dual cells." <<
endl
617 <<
"point " << f[fp] <<
" has dual cells "
618 << currentDualCell0 <<
" and " << currentDualCell1
619 <<
" ; point "<< f[nextFp] <<
" has dual cells "
620 << dualCell0 <<
" and " << dualCell1
648 void Foam::meshDualiser::createFacesAroundBoundaryPoint
651 const label patchPointI,
652 const label startFaceI,
653 polyTopoChange& meshMod,
657 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
658 const polyPatch& pp = patches[patchI];
660 const labelList& own = mesh_.faceOwner();
662 label pointI = pp.meshPoints()[patchPointI];
664 if (pointToDualPoint_[pointI] == -1)
670 label faceI = startFaceI;
672 DynamicList<label> verts(4);
676 label index =
findIndex(pFaces, faceI-pp.start());
679 if (donePFaces[index])
683 donePFaces[index] =
true;
686 verts.append(faceToDualPoint_[faceI]);
688 label dualCellI = findDualCell(own[faceI], pointI);
691 const face& f = mesh_.faces()[faceI];
693 label prevFp = f.rcIndex(fp);
694 label edgeI = mesh_.faceEdges()[faceI][prevFp];
696 if (edgeToDualPoint_[edgeI] != -1)
698 verts.append(edgeToDualPoint_[edgeI]);
702 edgeFaceCirculator circ
715 while (mesh_.isInternalFace(circ.faceLabel()));
718 faceI = circ.faceLabel();
720 if (faceI < pp.start() || faceI >= pp.start()+pp.size())
724 "meshDualiser::createFacesAroundBoundaryPoint(..)"
725 ) <<
"Walked from face on patch:" << patchI
726 <<
" to face:" << faceI
727 <<
" fc:" << mesh_.faceCentres()[faceI]
728 <<
" on patch:" << patches.whichPatch(faceI)
733 if (dualCellI != findDualCell(own[faceI], pointI))
737 "meshDualiser::createFacesAroundBoundaryPoint(..)"
738 ) <<
"Different dual cells but no feature edge"
739 <<
" inbetween point:" << pointI
740 <<
" coord:" << mesh_.points()[pointI]
747 label dualCellI = findDualCell(own[faceI], pointI);
767 label faceI = startFaceI;
770 DynamicList<label> verts(mesh_.faces()[faceI].size());
773 verts.append(pointToDualPoint_[pointI]);
776 const labelList& fEdges = mesh_.faceEdges()[faceI];
777 label nextEdgeI = fEdges[
findIndex(mesh_.faces()[faceI], pointI)];
778 if (edgeToDualPoint_[nextEdgeI] != -1)
780 verts.
append(edgeToDualPoint_[nextEdgeI]);
785 label index =
findIndex(pFaces, faceI-pp.start());
788 if (donePFaces[index])
792 donePFaces[index] =
true;
795 verts.append(faceToDualPoint_[faceI]);
798 const labelList& fEdges = mesh_.faceEdges()[faceI];
799 const face& f = mesh_.faces()[faceI];
801 label edgeI = fEdges[prevFp];
803 if (edgeToDualPoint_[edgeI] != -1)
807 verts.append(edgeToDualPoint_[edgeI]);
813 findDualCell(own[faceI], pointI),
820 verts.append(pointToDualPoint_[pointI]);
821 verts.append(edgeToDualPoint_[edgeI]);
825 edgeFaceCirculator circ
838 while (mesh_.isInternalFace(circ.faceLabel()));
841 faceI = circ.faceLabel();
846 && faceI >= pp.start()
847 && faceI < pp.start()+pp.size()
850 if (verts.size() > 2)
858 findDualCell(own[faceI], pointI),
870 Foam::meshDualiser::meshDualiser(
const polyMesh&
mesh)
873 pointToDualCells_(mesh_.
nPoints()),
874 pointToDualPoint_(mesh_.
nPoints(), -1),
875 cellToDualPoint_(mesh_.nCells()),
876 faceToDualPoint_(mesh_.nFaces(), -1),
877 edgeToDualPoint_(mesh_.nEdges(), -1)
888 const bool splitFace,
891 const labelList& singleCellFeaturePoints,
893 polyTopoChange& meshMod
910 isBoundaryEdge.set(fEdges[i], 1);
927 featureFaceSet[featureFaces[i]] =
true;
929 label faceI =
findIndex(featureFaceSet,
false);
935 "meshDualiser::setRefinement"
936 "(const labelList&, const labelList&, const labelList&"
937 ", const labelList, polyTopoChange&)"
938 ) <<
"In split-face-mode (splitFace=true) but not all faces"
939 <<
" marked as feature faces." <<
endl
940 <<
"First conflicting face:" << faceI
948 featureEdgeSet[featureEdges[i]] =
true;
950 label edgeI =
findIndex(featureEdgeSet,
false);
954 const edge& e = mesh_.
edges()[edgeI];
957 "meshDualiser::setRefinement"
958 "(const labelList&, const labelList&, const labelList&"
959 ", const labelList, polyTopoChange&)"
960 ) <<
"In split-face-mode (splitFace=true) but not all edges"
961 <<
" marked as feature edges." <<
endl
962 <<
"First conflicting edge:" << edgeI
964 <<
" coords:" << mesh_.
points()[e[0]] << mesh_.
points()[e[1]]
975 featureFaceSet[featureFaces[i]] =
true;
984 if (!featureFaceSet[faceI])
988 "meshDualiser::setRefinement"
989 "(const labelList&, const labelList&, const labelList&"
990 ", const labelList, polyTopoChange&)"
991 ) <<
"Not all boundary faces marked as feature faces."
993 <<
"First conflicting face:" << faceI
1024 autoPtr<OFstream> dualCcStr;
1027 dualCcStr.reset(
new OFstream(
"dualCc.obj"));
1028 Pout<<
"Dumping centres of dual cells to " << dualCcStr().
name()
1043 forAll(singleCellFeaturePoints, i)
1045 label pointI = singleCellFeaturePoints[i];
1047 pointToDualPoint_[pointI] = meshMod.addPoint
1056 pointToDualCells_[pointI].
setSize(1);
1057 pointToDualCells_[pointI][0] = meshMod.addCell
1065 if (dualCcStr.valid())
1072 forAll(multiCellFeaturePoints, i)
1074 label pointI = multiCellFeaturePoints[i];
1076 if (pointToDualCells_[pointI].size() > 0)
1080 "meshDualiser::setRefinement"
1081 "(const labelList&, const labelList&, const labelList&"
1082 ", const labelList, polyTopoChange&)"
1083 ) <<
"Point " << pointI <<
" at:" << mesh_.
points()[pointI]
1084 <<
" is both in singleCellFeaturePoints"
1085 <<
" and multiCellFeaturePoints."
1089 pointToDualPoint_[pointI] = meshMod.addPoint
1101 pointToDualCells_[pointI].
setSize(pCells.size());
1105 pointToDualCells_[pointI][pCellI] = meshMod.addCell
1113 if (dualCcStr.valid())
1118 0.5*(mesh_.
points()[pointI]+cellCentres[pCells[pCellI]])
1126 if (pointToDualCells_[pointI].empty())
1128 pointToDualCells_[pointI].
setSize(1);
1129 pointToDualCells_[pointI][0] = meshMod.addCell
1138 if (dualCcStr.valid())
1149 forAll(cellToDualPoint_, cellI)
1151 cellToDualPoint_[cellI] = meshMod.addPoint
1164 label faceI = featureFaces[i];
1166 faceToDualPoint_[faceI] = meshMod.addPoint
1169 mesh_.
faces()[faceI][0],
1179 if (faceToDualPoint_[faceI] == -1)
1181 const face& f = mesh_.
faces()[faceI];
1185 label ownDualCell = findDualCell(own[faceI], f[fp]);
1186 label neiDualCell = findDualCell(nei[faceI], f[fp]);
1188 if (ownDualCell != neiDualCell)
1190 faceToDualPoint_[faceI] = meshMod.addPoint
1208 label edgeI = featureEdges[i];
1210 const edge& e = mesh_.
edges()[edgeI];
1212 edgeToDualPoint_[edgeI] = meshMod.addPoint
1214 e.centre(mesh_.
points()),
1229 if (edgeToDualPoint_[edgeI] == -1)
1231 const edge& e = mesh_.
edges()[edgeI];
1238 label dualE0 = findDualCell(eCells[0], e[0]);
1239 label dualE1 = findDualCell(eCells[0], e[1]);
1241 for (label i = 1; i < eCells.size(); i++)
1243 label newDualE0 = findDualCell(eCells[i], e[0]);
1245 if (dualE0 != newDualE0)
1247 edgeToDualPoint_[edgeI] = meshMod.addPoint
1249 e.centre(mesh_.
points()),
1258 label newDualE1 = findDualCell(eCells[i], e[1]);
1260 if (dualE1 != newDualE1)
1262 edgeToDualPoint_[edgeI] = meshMod.addPoint
1264 e.centre(mesh_.
points()),
1278 forAll(singleCellFeaturePoints, i)
1280 generateDualBoundaryEdges
1283 singleCellFeaturePoints[i],
1287 forAll(multiCellFeaturePoints, i)
1289 generateDualBoundaryEdges
1292 multiCellFeaturePoints[i],
1301 dumpPolyTopoChange(meshMod,
"generatedPoints_");
1302 checkPolyTopoChange(meshMod);
1322 boolList doneEFaces(eFaces.size(),
false);
1332 label startFaceI = eFaces[i];
1341 createFacesAroundEdge
1356 dumpPolyTopoChange(meshMod,
"generatedFacesFromEdges_");
1365 forAll(faceToDualPoint_, faceI)
1367 if (faceToDualPoint_[faceI] != -1 && mesh_.
isInternalFace(faceI))
1369 const face& f = mesh_.
faces()[faceI];
1379 bool foundStart =
false;
1385 edgeToDualPoint_[fEdges[fp]] != -1
1386 && !sameDualCell(faceI, f.nextLabel(fp))
1402 createFaceFromInternalFace
1415 dumpPolyTopoChange(meshMod,
"generatedFacesFromFeatFaces_");
1421 const polyBoundaryMesh& patches = mesh_.
boundaryMesh();
1425 const polyPatch& pp = patches[patchI];
1429 forAll(pointFaces, patchPointI)
1431 const labelList& pFaces = pointFaces[patchPointI];
1433 boolList donePFaces(pFaces.size(),
false);
1440 label startFaceI = pp.start()+pFaces[i];
1449 createFacesAroundBoundaryPoint
1464 dumpPolyTopoChange(meshMod,
"generatedFacesFromBndFaces_");