48 void Foam::meshDualiser::checkPolyTopoChange(
const polyTopoChange& meshMod)
52 forAll(meshMod.points(), i)
54 points[i] = meshMod.points()[i];
66 if (nUnique <
points.size())
68 CompactListList<label> newToOld
75 if (newToOld[newi].size() != 1)
78 <<
"duplicate verts:" << newToOld[newi]
80 << UIndirectList<point>(
points, newToOld[newi])
89 void Foam::meshDualiser::dumpPolyTopoChange
91 const polyTopoChange& meshMod,
92 const fileName& prefix
95 OFstream str1(prefix +
"Faces.obj");
96 OFstream str2(prefix +
"Edges.obj");
98 Info<<
"Dumping current polyTopoChange. Faces to " << str1.name()
99 <<
" , points and edges to " << str2.name() <<
endl;
101 const DynamicList<point>&
points = meshMod.points();
109 const DynamicList<face>& faces = meshMod.faces();
113 const face&
f = faces[facei];
118 str1<<
' ' <<
f[fp]+1;
125 str2<<
' ' <<
f[fp]+1;
127 str2<<
' ' <<
f[0]+1 <<
nl;
132 Foam::label Foam::meshDualiser::findDualCell
138 const labelList& dualCells = pointToDualCells_[pointi];
140 if (dualCells.size() == 1)
146 label index = mesh_.pointCells()[pointi].
find(celli);
148 return dualCells[index];
153 void Foam::meshDualiser::generateDualBoundaryEdges
155 const bitSet& isBoundaryEdge,
157 polyTopoChange& meshMod
160 const labelList& pEdges = mesh_.pointEdges()[pointi];
164 label edgeI = pEdges[pEdgeI];
166 if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.test(edgeI))
168 const edge&
e = mesh_.edges()[edgeI];
170 edgeToDualPoint_[edgeI] = meshMod.addPoint
172 e.centre(mesh_.points()),
184 bool Foam::meshDualiser::sameDualCell
190 if (!mesh_.isInternalFace(facei))
193 <<
"face:" << facei <<
" is not internal face." 197 label own = mesh_.faceOwner()[facei];
198 label nei = mesh_.faceNeighbour()[facei];
200 return findDualCell(own, pointi) == findDualCell(nei, pointi);
204 Foam::label Foam::meshDualiser::addInternalFace
206 const label masterPointi,
207 const label masterEdgeI,
208 const label masterFacei,
210 const bool edgeOrder,
211 const label dualCell0,
212 const label dualCell1,
213 const DynamicList<label>& verts,
214 polyTopoChange& meshMod
219 if (edgeOrder != (dualCell0 < dualCell1))
226 pointField facePoints(meshMod.points(), newFace);
237 if (nUnique < facePoints.size())
240 <<
"verts:" << verts <<
" newFace:" << newFace
241 <<
" face points:" << facePoints
248 bool zoneFlip =
false;
249 if (masterFacei != -1)
251 zoneID = mesh_.faceZones().whichZone(masterFacei);
255 const faceZone& fZone = mesh_.faceZones()[zoneID];
257 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
263 if (dualCell0 < dualCell1)
265 dualFacei = meshMod.addFace
292 dualFacei = meshMod.addFace
321 Foam::label Foam::meshDualiser::addBoundaryFace
323 const label masterPointi,
324 const label masterEdgeI,
325 const label masterFacei,
327 const label dualCelli,
329 const DynamicList<label>& verts,
330 polyTopoChange& meshMod
336 bool zoneFlip =
false;
337 if (masterFacei != -1)
339 zoneID = mesh_.faceZones().whichZone(masterFacei);
343 const faceZone& fZone = mesh_.faceZones()[zoneID];
345 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
349 label dualFacei = meshMod.addFace
380 void Foam::meshDualiser::createFacesAroundEdge
382 const bool splitFace,
383 const bitSet& isBoundaryEdge,
385 const label startFacei,
386 polyTopoChange& meshMod,
390 const edge&
e = mesh_.edges()[edgeI];
391 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
395 mesh_.faces()[startFacei],
400 edgeFaceCirculator ie
406 isBoundaryEdge.test(edgeI)
410 bool edgeOrder = ie.sameOrder(
e[0],
e[1]);
411 label startFaceLabel = ie.faceLabel();
422 DynamicList<label> verts(100);
424 if (edgeToDualPoint_[edgeI] != -1)
426 verts.append(edgeToDualPoint_[edgeI]);
428 if (faceToDualPoint_[ie.faceLabel()] != -1)
430 doneEFaces[eFaces.find(ie.faceLabel())] =
true;
431 verts.append(faceToDualPoint_[ie.faceLabel()]);
433 if (cellToDualPoint_[ie.cellLabel()] != -1)
435 verts.append(cellToDualPoint_[ie.cellLabel()]);
438 label currentDualCell0 = findDualCell(ie.cellLabel(),
e[0]);
439 label currentDualCell1 = findDualCell(ie.cellLabel(),
e[1]);
445 label facei = ie.faceLabel();
448 doneEFaces[eFaces.find(facei)] =
true;
450 if (faceToDualPoint_[facei] != -1)
452 verts.append(faceToDualPoint_[facei]);
455 label celli = ie.cellLabel();
465 label dualCell0 = findDualCell(celli,
e[0]);
466 label dualCell1 = findDualCell(celli,
e[1]);
474 dualCell0 != currentDualCell0
475 || dualCell1 != currentDualCell1
493 currentDualCell0 = dualCell0;
494 currentDualCell1 = dualCell1;
497 if (edgeToDualPoint_[edgeI] != -1)
499 verts.append(edgeToDualPoint_[edgeI]);
501 if (faceToDualPoint_[facei] != -1)
503 verts.append(faceToDualPoint_[facei]);
507 if (cellToDualPoint_[celli] != -1)
509 verts.append(cellToDualPoint_[celli]);
518 if (!isBoundaryEdge.test(edgeI))
520 label startDual = faceToDualPoint_[startFaceLabel];
524 verts.push_uniq(startDual);
549 void Foam::meshDualiser::createFaceFromInternalFace
553 polyTopoChange& meshMod
556 const face&
f = mesh_.faces()[facei];
557 const labelList& fEdges = mesh_.faceEdges()[facei];
558 label own = mesh_.faceOwner()[facei];
559 label nei = mesh_.faceNeighbour()[facei];
570 DynamicList<label> verts(100);
572 verts.append(faceToDualPoint_[facei]);
573 verts.append(edgeToDualPoint_[fEdges[fp]]);
579 label currentDualCell0 = findDualCell(own,
f[fp]);
580 label currentDualCell1 = findDualCell(nei,
f[fp]);
585 if (pointToDualPoint_[
f[fp]] != -1)
587 verts.append(pointToDualPoint_[
f[fp]]);
591 label edgeI = fEdges[fp];
593 if (edgeToDualPoint_[edgeI] != -1)
595 verts.append(edgeToDualPoint_[edgeI]);
602 label dualCell0 = findDualCell(own,
f[nextFp]);
603 label dualCell1 = findDualCell(nei,
f[nextFp]);
605 if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
608 if (edgeToDualPoint_[edgeI] == -1)
611 <<
"face:" << facei <<
" verts:" <<
f 612 <<
" points:" << UIndirectList<point>(mesh_.points(),
f)
613 <<
" no feature edge between " <<
f[fp]
614 <<
" and " <<
f[nextFp] <<
" although have different" 615 <<
" dual cells." <<
endl 616 <<
"point " <<
f[fp] <<
" has dual cells " 617 << currentDualCell0 <<
" and " << currentDualCell1
618 <<
" ; point "<<
f[nextFp] <<
" has dual cells " 619 << dualCell0 <<
" and " << dualCell1
647 void Foam::meshDualiser::createFacesAroundBoundaryPoint
650 const label patchPointi,
651 const label startFacei,
652 polyTopoChange& meshMod,
656 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
659 const labelList& own = mesh_.faceOwner();
663 if (pointToDualPoint_[pointi] == -1)
669 label facei = startFacei;
671 DynamicList<label> verts(4);
675 label index =
pFaces.find(facei-
pp.start());
678 if (donePFaces[index])
682 donePFaces[index] =
true;
685 verts.append(faceToDualPoint_[facei]);
687 label dualCelli = findDualCell(own[facei], pointi);
690 const face&
f = mesh_.faces()[facei];
691 label fp =
f.
find(pointi);
693 label edgeI = mesh_.faceEdges()[facei][prevFp];
695 if (edgeToDualPoint_[edgeI] != -1)
697 verts.append(edgeToDualPoint_[edgeI]);
701 edgeFaceCirculator circ
714 while (mesh_.isInternalFace(circ.faceLabel()));
717 facei = circ.faceLabel();
719 if (facei <
pp.start() || facei >=
pp.start()+
pp.size())
722 <<
"Walked from face on patch:" << patchi
723 <<
" to face:" << facei
724 <<
" fc:" << mesh_.faceCentres()[facei]
730 if (dualCelli != findDualCell(own[facei], pointi))
733 <<
"Different dual cells but no feature edge" 734 <<
" inbetween point:" << pointi
735 <<
" coord:" << mesh_.points()[pointi]
742 label dualCelli = findDualCell(own[facei], pointi);
762 label facei = startFacei;
765 DynamicList<label> verts(mesh_.faces()[facei].size());
768 verts.append(pointToDualPoint_[pointi]);
771 const labelList& fEdges = mesh_.faceEdges()[facei];
772 label nextEdgeI = fEdges[mesh_.faces()[facei].find(pointi)];
773 if (edgeToDualPoint_[nextEdgeI] != -1)
775 verts.
append(edgeToDualPoint_[nextEdgeI]);
780 label index =
pFaces.find(facei-
pp.start());
783 if (donePFaces[index])
787 donePFaces[index] =
true;
790 verts.append(faceToDualPoint_[facei]);
793 const labelList& fEdges = mesh_.faceEdges()[facei];
794 const face&
f = mesh_.faces()[facei];
796 label edgeI = fEdges[prevFp];
798 if (edgeToDualPoint_[edgeI] != -1)
802 verts.append(edgeToDualPoint_[edgeI]);
808 findDualCell(own[facei], pointi),
815 verts.append(pointToDualPoint_[pointi]);
816 verts.append(edgeToDualPoint_[edgeI]);
820 edgeFaceCirculator circ
833 while (mesh_.isInternalFace(circ.faceLabel()));
836 facei = circ.faceLabel();
841 && facei >=
pp.start()
842 && facei <
pp.start()+
pp.size()
845 if (verts.size() > 2)
853 findDualCell(own[facei], pointi),
865 Foam::meshDualiser::meshDualiser(
const polyMesh&
mesh)
868 pointToDualCells_(mesh_.
nPoints()),
869 pointToDualPoint_(mesh_.
nPoints(), -1),
870 cellToDualPoint_(mesh_.nCells()),
871 faceToDualPoint_(mesh_.nFaces(), -1),
872 edgeToDualPoint_(mesh_.nEdges(), -1)
880 const bool splitFace,
883 const labelList& singleCellFeaturePoints,
885 polyTopoChange& meshMod
888 const labelList& own = mesh_.faceOwner();
889 const labelList& nei = mesh_.faceNeighbour();
890 const vectorField& cellCentres = mesh_.cellCentres();
895 bitSet isBoundaryEdge(mesh_.nEdges());
896 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
898 const labelList& fEdges = mesh_.faceEdges()[facei];
900 isBoundaryEdge.
set(fEdges);
913 boolList featureFaceSet(mesh_.nFaces(),
false);
916 featureFaceSet[featureFaces[i]] =
true;
918 label facei = featureFaceSet.find(
false);
923 <<
"In split-face-mode (splitFace=true) but not all faces" 924 <<
" marked as feature faces." <<
endl 925 <<
"First conflicting face:" << facei
926 <<
" centre:" << mesh_.faceCentres()[facei]
930 boolList featureEdgeSet(mesh_.nEdges(),
false);
933 featureEdgeSet[featureEdges[i]] =
true;
935 label edgeI = featureEdgeSet.find(
false);
939 const edge&
e = mesh_.edges()[edgeI];
941 <<
"In split-face-mode (splitFace=true) but not all edges" 942 <<
" marked as feature edges." <<
endl 943 <<
"First conflicting edge:" << edgeI
945 <<
" coords:" << mesh_.points()[
e[0]] << mesh_.points()[
e[1]]
953 boolList featureFaceSet(mesh_.nFaces(),
false);
956 featureFaceSet[featureFaces[i]] =
true;
960 label facei = mesh_.nInternalFaces();
961 facei < mesh_.nFaces();
965 if (!featureFaceSet[facei])
968 <<
"Not all boundary faces marked as feature faces." 970 <<
"First conflicting face:" << facei
971 <<
" centre:" << mesh_.faceCentres()[facei]
1001 autoPtr<OFstream> dualCcStr;
1004 dualCcStr.reset(
new OFstream(
"dualCc.obj"));
1005 Pout<<
"Dumping centres of dual cells to " << dualCcStr().
name()
1020 forAll(singleCellFeaturePoints, i)
1022 label pointi = singleCellFeaturePoints[i];
1024 pointToDualPoint_[pointi] = meshMod.addPoint
1026 mesh_.points()[pointi],
1028 mesh_.pointZones().whichZone(pointi),
1033 pointToDualCells_[pointi].setSize(1);
1034 pointToDualCells_[pointi][0] = meshMod.addCell
1049 forAll(multiCellFeaturePoints, i)
1051 label pointi = multiCellFeaturePoints[i];
1053 if (pointToDualCells_[pointi].size() > 0)
1056 <<
"Point " << pointi <<
" at:" << mesh_.points()[pointi]
1057 <<
" is both in singleCellFeaturePoints" 1058 <<
" and multiCellFeaturePoints." 1062 pointToDualPoint_[pointi] = meshMod.addPoint
1064 mesh_.points()[pointi],
1066 mesh_.pointZones().whichZone(pointi),
1072 const labelList& pCells = mesh_.pointCells()[pointi];
1074 pointToDualCells_[pointi].
setSize(pCells.size());
1078 pointToDualCells_[pointi][pCelli] = meshMod.addCell
1084 mesh_.cellZones().whichZone(pCells[pCelli])
1091 0.5*(mesh_.points()[pointi]+cellCentres[pCells[pCelli]])
1097 forAll(mesh_.points(), pointi)
1099 if (pointToDualCells_[pointi].empty())
1101 pointToDualCells_[pointi].setSize(1);
1102 pointToDualCells_[pointi][0] = meshMod.addCell
1122 forAll(cellToDualPoint_, celli)
1124 cellToDualPoint_[celli] = meshMod.addPoint
1127 mesh_.faces()[mesh_.cells()[celli][0]][0],
1137 label facei = featureFaces[i];
1139 faceToDualPoint_[facei] = meshMod.addPoint
1141 mesh_.faceCentres()[facei],
1142 mesh_.faces()[facei][0],
1150 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1152 if (faceToDualPoint_[facei] == -1)
1154 const face&
f = mesh_.faces()[facei];
1158 label ownDualCell = findDualCell(own[facei],
f[fp]);
1159 label neiDualCell = findDualCell(nei[facei],
f[fp]);
1161 if (ownDualCell != neiDualCell)
1163 faceToDualPoint_[facei] = meshMod.addPoint
1165 mesh_.faceCentres()[facei],
1181 label edgeI = featureEdges[i];
1183 const edge&
e = mesh_.edges()[edgeI];
1185 edgeToDualPoint_[edgeI] = meshMod.addPoint
1187 e.centre(mesh_.points()),
1202 if (edgeToDualPoint_[edgeI] == -1)
1204 const edge&
e = mesh_.edges()[edgeI];
1209 const labelList& eCells = mesh_.edgeCells()[edgeI];
1211 label dualE0 = findDualCell(eCells[0],
e[0]);
1212 label dualE1 = findDualCell(eCells[0],
e[1]);
1214 for (label i = 1; i < eCells.size(); i++)
1216 label newDualE0 = findDualCell(eCells[i],
e[0]);
1218 if (dualE0 != newDualE0)
1220 edgeToDualPoint_[edgeI] = meshMod.addPoint
1222 e.centre(mesh_.points()),
1224 mesh_.pointZones().whichZone(
e[0]),
1231 label newDualE1 = findDualCell(eCells[i],
e[1]);
1233 if (dualE1 != newDualE1)
1235 edgeToDualPoint_[edgeI] = meshMod.addPoint
1237 e.centre(mesh_.points()),
1239 mesh_.pointZones().whichZone(
e[1]),
1251 forAll(singleCellFeaturePoints, i)
1253 generateDualBoundaryEdges
1256 singleCellFeaturePoints[i],
1260 forAll(multiCellFeaturePoints, i)
1262 generateDualBoundaryEdges
1265 multiCellFeaturePoints[i],
1274 dumpPolyTopoChange(meshMod,
"generatedPoints_");
1275 checkPolyTopoChange(meshMod);
1289 const edgeList& edges = mesh_.edges();
1293 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
1295 boolList doneEFaces(eFaces.size(),
false);
1305 label startFacei = eFaces[i];
1314 createFacesAroundEdge
1329 dumpPolyTopoChange(meshMod,
"generatedFacesFromEdges_");
1338 forAll(faceToDualPoint_, facei)
1340 if (faceToDualPoint_[facei] != -1 && mesh_.isInternalFace(facei))
1342 const face&
f = mesh_.faces()[facei];
1343 const labelList& fEdges = mesh_.faceEdges()[facei];
1352 bool foundStart =
false;
1358 edgeToDualPoint_[fEdges[fp]] != -1
1359 && !sameDualCell(facei,
f.nextLabel(fp))
1375 createFaceFromInternalFace
1388 dumpPolyTopoChange(meshMod,
"generatedFacesFromFeatFaces_");
1394 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
1402 forAll(pointFaces, patchPointi)
1413 label startFacei =
pp.start()+
pFaces[i];
1422 createFacesAroundBoundaryPoint
1437 dumpPolyTopoChange(meshMod,
"generatedFacesFromBndFaces_");
void setRefinement(const bool splitFace, const labelList &featureFaces, const labelList &featureEdges, const labelList &singleCellFeaturePoints, const labelList &multiCellFeaturePoints, polyTopoChange &meshMod)
Insert all changes into meshMod to convert the polyMesh into.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void append(const T &val)
Append an element at the end of the list.
List< edge > edgeList
List of edge.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
static label getMinIndex(const face &f, const label v0, const label v1)
Helper: find index in face of edge or -1. Index is such that edge is.
virtual const fileName & name() const override
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
List of labelList.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
label mergePoints(const PointList &points, labelList &pointToUnique, labelList &uniquePoints, const scalar mergeTol=SMALL, const bool verbose=false)
Calculate merge mapping, preserving the original point order. All points closer/equal mergeTol are to...
CompactListList< label > invertOneToManyCompact(const label len, const labelUList &map)
Invert one-to-many compact map. Unmapped elements will be size 0.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
errorManip< error > abort(error &err)
label find(const T &val) const
Find index of the first occurrence of the value.
const labelListList & pointFaces() const
Return point-face addressing.
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
Geometric merging of points. See below.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
const polyBoundaryMesh & patches
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< label > labelList
A List of labels.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())