48 void Foam::meshDualiser::checkPolyTopoChange(
const polyTopoChange& meshMod)
52 forAll(meshMod.points(), i)
54 points[i] = meshMod.points()[i];
66 if (nUnique <
points.size())
72 if (newToOld[newI].size() != 1)
75 <<
"duplicate verts:" << newToOld[newI]
77 << UIndirectList<point>(
points, newToOld[newI])
86 void Foam::meshDualiser::dumpPolyTopoChange
88 const polyTopoChange& meshMod,
89 const fileName& prefix
92 OFstream str1(prefix +
"Faces.obj");
93 OFstream str2(prefix +
"Edges.obj");
95 Info<<
"Dumping current polyTopoChange. Faces to " << str1.name()
96 <<
" , points and edges to " << str2.name() <<
endl;
98 const DynamicList<point>&
points = meshMod.points();
106 const DynamicList<face>& faces = meshMod.faces();
110 const face&
f = faces[facei];
115 str1<<
' ' <<
f[fp]+1;
122 str2<<
' ' <<
f[fp]+1;
124 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 = mesh_.pointCells()[pointi].
find(celli);
145 return dualCells[index];
150 void Foam::meshDualiser::generateDualBoundaryEdges
152 const bitSet& isBoundaryEdge,
154 polyTopoChange& meshMod
157 const labelList& pEdges = mesh_.pointEdges()[pointi];
161 label edgeI = pEdges[pEdgeI];
163 if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.test(edgeI))
165 const edge&
e = mesh_.edges()[edgeI];
167 edgeToDualPoint_[edgeI] = meshMod.addPoint
169 e.centre(mesh_.points()),
181 bool Foam::meshDualiser::sameDualCell
187 if (!mesh_.isInternalFace(facei))
190 <<
"face:" << facei <<
" is not internal face." 194 label own = mesh_.faceOwner()[facei];
195 label nei = mesh_.faceNeighbour()[facei];
197 return findDualCell(own, pointi) == findDualCell(nei, pointi);
201 Foam::label Foam::meshDualiser::addInternalFace
203 const label masterPointi,
204 const label masterEdgeI,
205 const label masterFacei,
207 const bool edgeOrder,
208 const label dualCell0,
209 const label dualCell1,
210 const DynamicList<label>& verts,
211 polyTopoChange& meshMod
216 if (edgeOrder != (dualCell0 < dualCell1))
223 pointField facePoints(meshMod.points(), newFace);
234 if (nUnique < facePoints.size())
237 <<
"verts:" << verts <<
" newFace:" << newFace
238 <<
" face points:" << facePoints
245 bool zoneFlip =
false;
246 if (masterFacei != -1)
248 zoneID = mesh_.faceZones().whichZone(masterFacei);
252 const faceZone& fZone = mesh_.faceZones()[zoneID];
254 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
260 if (dualCell0 < dualCell1)
262 dualFacei = meshMod.addFace
289 dualFacei = meshMod.addFace
318 Foam::label Foam::meshDualiser::addBoundaryFace
320 const label masterPointi,
321 const label masterEdgeI,
322 const label masterFacei,
324 const label dualCelli,
326 const DynamicList<label>& verts,
327 polyTopoChange& meshMod
333 bool zoneFlip =
false;
334 if (masterFacei != -1)
336 zoneID = mesh_.faceZones().whichZone(masterFacei);
340 const faceZone& fZone = mesh_.faceZones()[zoneID];
342 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
346 label dualFacei = meshMod.addFace
377 void Foam::meshDualiser::createFacesAroundEdge
379 const bool splitFace,
380 const bitSet& isBoundaryEdge,
382 const label startFacei,
383 polyTopoChange& meshMod,
387 const edge&
e = mesh_.edges()[edgeI];
388 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
392 mesh_.faces()[startFacei],
397 edgeFaceCirculator ie
403 isBoundaryEdge.test(edgeI)
407 bool edgeOrder = ie.sameOrder(
e[0],
e[1]);
408 label startFaceLabel = ie.faceLabel();
419 DynamicList<label> verts(100);
421 if (edgeToDualPoint_[edgeI] != -1)
423 verts.append(edgeToDualPoint_[edgeI]);
425 if (faceToDualPoint_[ie.faceLabel()] != -1)
427 doneEFaces[eFaces.find(ie.faceLabel())] =
true;
428 verts.append(faceToDualPoint_[ie.faceLabel()]);
430 if (cellToDualPoint_[ie.cellLabel()] != -1)
432 verts.append(cellToDualPoint_[ie.cellLabel()]);
435 label currentDualCell0 = findDualCell(ie.cellLabel(),
e[0]);
436 label currentDualCell1 = findDualCell(ie.cellLabel(),
e[1]);
442 label facei = ie.faceLabel();
445 doneEFaces[eFaces.find(facei)] =
true;
447 if (faceToDualPoint_[facei] != -1)
449 verts.append(faceToDualPoint_[facei]);
452 label celli = ie.cellLabel();
462 label dualCell0 = findDualCell(celli,
e[0]);
463 label dualCell1 = findDualCell(celli,
e[1]);
471 dualCell0 != currentDualCell0
472 || dualCell1 != currentDualCell1
490 currentDualCell0 = dualCell0;
491 currentDualCell1 = dualCell1;
494 if (edgeToDualPoint_[edgeI] != -1)
496 verts.append(edgeToDualPoint_[edgeI]);
498 if (faceToDualPoint_[facei] != -1)
500 verts.append(faceToDualPoint_[facei]);
504 if (cellToDualPoint_[celli] != -1)
506 verts.append(cellToDualPoint_[celli]);
515 if (!isBoundaryEdge.test(edgeI))
517 label startDual = faceToDualPoint_[startFaceLabel];
521 verts.push_uniq(startDual);
546 void Foam::meshDualiser::createFaceFromInternalFace
550 polyTopoChange& meshMod
553 const face&
f = mesh_.faces()[facei];
554 const labelList& fEdges = mesh_.faceEdges()[facei];
555 label own = mesh_.faceOwner()[facei];
556 label nei = mesh_.faceNeighbour()[facei];
567 DynamicList<label> verts(100);
569 verts.append(faceToDualPoint_[facei]);
570 verts.append(edgeToDualPoint_[fEdges[fp]]);
576 label currentDualCell0 = findDualCell(own,
f[fp]);
577 label currentDualCell1 = findDualCell(nei,
f[fp]);
582 if (pointToDualPoint_[
f[fp]] != -1)
584 verts.append(pointToDualPoint_[
f[fp]]);
588 label edgeI = fEdges[fp];
590 if (edgeToDualPoint_[edgeI] != -1)
592 verts.append(edgeToDualPoint_[edgeI]);
599 label dualCell0 = findDualCell(own,
f[nextFp]);
600 label dualCell1 = findDualCell(nei,
f[nextFp]);
602 if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
605 if (edgeToDualPoint_[edgeI] == -1)
608 <<
"face:" << facei <<
" verts:" <<
f 609 <<
" points:" << UIndirectList<point>(mesh_.points(),
f)
610 <<
" no feature edge between " <<
f[fp]
611 <<
" and " <<
f[nextFp] <<
" although have different" 612 <<
" dual cells." <<
endl 613 <<
"point " <<
f[fp] <<
" has dual cells " 614 << currentDualCell0 <<
" and " << currentDualCell1
615 <<
" ; point "<<
f[nextFp] <<
" has dual cells " 616 << dualCell0 <<
" and " << dualCell1
644 void Foam::meshDualiser::createFacesAroundBoundaryPoint
647 const label patchPointi,
648 const label startFacei,
649 polyTopoChange& meshMod,
653 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
656 const labelList& own = mesh_.faceOwner();
660 if (pointToDualPoint_[pointi] == -1)
666 label facei = startFacei;
668 DynamicList<label> verts(4);
672 label index =
pFaces.find(facei-
pp.start());
675 if (donePFaces[index])
679 donePFaces[index] =
true;
682 verts.append(faceToDualPoint_[facei]);
684 label dualCelli = findDualCell(own[facei], pointi);
687 const face&
f = mesh_.faces()[facei];
688 label fp =
f.
find(pointi);
690 label edgeI = mesh_.faceEdges()[facei][prevFp];
692 if (edgeToDualPoint_[edgeI] != -1)
694 verts.append(edgeToDualPoint_[edgeI]);
698 edgeFaceCirculator circ
711 while (mesh_.isInternalFace(circ.faceLabel()));
714 facei = circ.faceLabel();
716 if (facei <
pp.start() || facei >=
pp.start()+
pp.size())
719 <<
"Walked from face on patch:" << patchi
720 <<
" to face:" << facei
721 <<
" fc:" << mesh_.faceCentres()[facei]
727 if (dualCelli != findDualCell(own[facei], pointi))
730 <<
"Different dual cells but no feature edge" 731 <<
" inbetween point:" << pointi
732 <<
" coord:" << mesh_.points()[pointi]
739 label dualCelli = findDualCell(own[facei], pointi);
759 label facei = startFacei;
762 DynamicList<label> verts(mesh_.faces()[facei].size());
765 verts.append(pointToDualPoint_[pointi]);
768 const labelList& fEdges = mesh_.faceEdges()[facei];
769 label nextEdgeI = fEdges[mesh_.faces()[facei].find(pointi)];
770 if (edgeToDualPoint_[nextEdgeI] != -1)
772 verts.
append(edgeToDualPoint_[nextEdgeI]);
777 label index =
pFaces.find(facei-
pp.start());
780 if (donePFaces[index])
784 donePFaces[index] =
true;
787 verts.append(faceToDualPoint_[facei]);
790 const labelList& fEdges = mesh_.faceEdges()[facei];
791 const face&
f = mesh_.faces()[facei];
793 label edgeI = fEdges[prevFp];
795 if (edgeToDualPoint_[edgeI] != -1)
799 verts.append(edgeToDualPoint_[edgeI]);
805 findDualCell(own[facei], pointi),
812 verts.append(pointToDualPoint_[pointi]);
813 verts.append(edgeToDualPoint_[edgeI]);
817 edgeFaceCirculator circ
830 while (mesh_.isInternalFace(circ.faceLabel()));
833 facei = circ.faceLabel();
838 && facei >=
pp.start()
839 && facei <
pp.start()+
pp.size()
842 if (verts.size() > 2)
850 findDualCell(own[facei], pointi),
862 Foam::meshDualiser::meshDualiser(
const polyMesh&
mesh)
865 pointToDualCells_(mesh_.
nPoints()),
866 pointToDualPoint_(mesh_.
nPoints(), -1),
867 cellToDualPoint_(mesh_.nCells()),
868 faceToDualPoint_(mesh_.nFaces(), -1),
869 edgeToDualPoint_(mesh_.nEdges(), -1)
877 const bool splitFace,
880 const labelList& singleCellFeaturePoints,
882 polyTopoChange& meshMod
885 const labelList& own = mesh_.faceOwner();
886 const labelList& nei = mesh_.faceNeighbour();
887 const vectorField& cellCentres = mesh_.cellCentres();
892 bitSet isBoundaryEdge(mesh_.nEdges());
893 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
895 const labelList& fEdges = mesh_.faceEdges()[facei];
897 isBoundaryEdge.
set(fEdges);
910 boolList featureFaceSet(mesh_.nFaces(),
false);
913 featureFaceSet[featureFaces[i]] =
true;
915 label facei = featureFaceSet.find(
false);
920 <<
"In split-face-mode (splitFace=true) but not all faces" 921 <<
" marked as feature faces." <<
endl 922 <<
"First conflicting face:" << facei
923 <<
" centre:" << mesh_.faceCentres()[facei]
927 boolList featureEdgeSet(mesh_.nEdges(),
false);
930 featureEdgeSet[featureEdges[i]] =
true;
932 label edgeI = featureEdgeSet.find(
false);
936 const edge&
e = mesh_.edges()[edgeI];
938 <<
"In split-face-mode (splitFace=true) but not all edges" 939 <<
" marked as feature edges." <<
endl 940 <<
"First conflicting edge:" << edgeI
942 <<
" coords:" << mesh_.points()[
e[0]] << mesh_.points()[
e[1]]
950 boolList featureFaceSet(mesh_.nFaces(),
false);
953 featureFaceSet[featureFaces[i]] =
true;
957 label facei = mesh_.nInternalFaces();
958 facei < mesh_.nFaces();
962 if (!featureFaceSet[facei])
965 <<
"Not all boundary faces marked as feature faces." 967 <<
"First conflicting face:" << facei
968 <<
" centre:" << mesh_.faceCentres()[facei]
998 autoPtr<OFstream> dualCcStr;
1001 dualCcStr.reset(
new OFstream(
"dualCc.obj"));
1002 Pout<<
"Dumping centres of dual cells to " << dualCcStr().
name()
1017 forAll(singleCellFeaturePoints, i)
1019 label pointi = singleCellFeaturePoints[i];
1021 pointToDualPoint_[pointi] = meshMod.addPoint
1023 mesh_.points()[pointi],
1025 mesh_.pointZones().whichZone(pointi),
1030 pointToDualCells_[pointi].setSize(1);
1031 pointToDualCells_[pointi][0] = meshMod.addCell
1046 forAll(multiCellFeaturePoints, i)
1048 label pointi = multiCellFeaturePoints[i];
1050 if (pointToDualCells_[pointi].size() > 0)
1053 <<
"Point " << pointi <<
" at:" << mesh_.points()[pointi]
1054 <<
" is both in singleCellFeaturePoints" 1055 <<
" and multiCellFeaturePoints." 1059 pointToDualPoint_[pointi] = meshMod.addPoint
1061 mesh_.points()[pointi],
1063 mesh_.pointZones().whichZone(pointi),
1069 const labelList& pCells = mesh_.pointCells()[pointi];
1071 pointToDualCells_[pointi].
setSize(pCells.size());
1075 pointToDualCells_[pointi][pCelli] = meshMod.addCell
1081 mesh_.cellZones().whichZone(pCells[pCelli])
1088 0.5*(mesh_.points()[pointi]+cellCentres[pCells[pCelli]])
1094 forAll(mesh_.points(), pointi)
1096 if (pointToDualCells_[pointi].empty())
1098 pointToDualCells_[pointi].setSize(1);
1099 pointToDualCells_[pointi][0] = meshMod.addCell
1119 forAll(cellToDualPoint_, celli)
1121 cellToDualPoint_[celli] = meshMod.addPoint
1124 mesh_.faces()[mesh_.cells()[celli][0]][0],
1134 label facei = featureFaces[i];
1136 faceToDualPoint_[facei] = meshMod.addPoint
1138 mesh_.faceCentres()[facei],
1139 mesh_.faces()[facei][0],
1147 for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1149 if (faceToDualPoint_[facei] == -1)
1151 const face&
f = mesh_.faces()[facei];
1155 label ownDualCell = findDualCell(own[facei],
f[fp]);
1156 label neiDualCell = findDualCell(nei[facei],
f[fp]);
1158 if (ownDualCell != neiDualCell)
1160 faceToDualPoint_[facei] = meshMod.addPoint
1162 mesh_.faceCentres()[facei],
1178 label edgeI = featureEdges[i];
1180 const edge&
e = mesh_.edges()[edgeI];
1182 edgeToDualPoint_[edgeI] = meshMod.addPoint
1184 e.centre(mesh_.points()),
1199 if (edgeToDualPoint_[edgeI] == -1)
1201 const edge&
e = mesh_.edges()[edgeI];
1206 const labelList& eCells = mesh_.edgeCells()[edgeI];
1208 label dualE0 = findDualCell(eCells[0],
e[0]);
1209 label dualE1 = findDualCell(eCells[0],
e[1]);
1211 for (label i = 1; i < eCells.size(); i++)
1213 label newDualE0 = findDualCell(eCells[i],
e[0]);
1215 if (dualE0 != newDualE0)
1217 edgeToDualPoint_[edgeI] = meshMod.addPoint
1219 e.centre(mesh_.points()),
1221 mesh_.pointZones().whichZone(
e[0]),
1228 label newDualE1 = findDualCell(eCells[i],
e[1]);
1230 if (dualE1 != newDualE1)
1232 edgeToDualPoint_[edgeI] = meshMod.addPoint
1234 e.centre(mesh_.points()),
1236 mesh_.pointZones().whichZone(
e[1]),
1248 forAll(singleCellFeaturePoints, i)
1250 generateDualBoundaryEdges
1253 singleCellFeaturePoints[i],
1257 forAll(multiCellFeaturePoints, i)
1259 generateDualBoundaryEdges
1262 multiCellFeaturePoints[i],
1271 dumpPolyTopoChange(meshMod,
"generatedPoints_");
1272 checkPolyTopoChange(meshMod);
1286 const edgeList& edges = mesh_.edges();
1290 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
1292 boolList doneEFaces(eFaces.size(),
false);
1302 label startFacei = eFaces[i];
1311 createFacesAroundEdge
1326 dumpPolyTopoChange(meshMod,
"generatedFacesFromEdges_");
1335 forAll(faceToDualPoint_, facei)
1337 if (faceToDualPoint_[facei] != -1 && mesh_.isInternalFace(facei))
1339 const face&
f = mesh_.faces()[facei];
1340 const labelList& fEdges = mesh_.faceEdges()[facei];
1349 bool foundStart =
false;
1355 edgeToDualPoint_[fEdges[fp]] != -1
1356 && !sameDualCell(facei,
f.nextLabel(fp))
1372 createFaceFromInternalFace
1385 dumpPolyTopoChange(meshMod,
"generatedFacesFromFeatFaces_");
1391 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
1399 forAll(pointFaces, patchPointi)
1410 label startFacei =
pp.start()+
pFaces[i];
1419 createFacesAroundBoundaryPoint
1434 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...
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
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())