48void Foam::meshDualiser::checkPolyTopoChange(
const polyTopoChange& meshMod)
52 forAll(meshMod.points(), i)
54 points[i] = meshMod.points()[i];
72 if (newToOld[newI].size() != 1)
75 <<
"duplicate verts:" << newToOld[newI]
77 << UIndirectList<point>(
points, newToOld[newI])
86void 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;
129Foam::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];
150void 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()),
181bool 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);
201Foam::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
318Foam::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
377void 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.appendUniq(startDual);
546void 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
644void Foam::meshDualiser::createFacesAroundBoundaryPoint
647 const label patchPointi,
648 const label startFacei,
649 polyTopoChange& meshMod,
653 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
654 const polyPatch& pp =
patches[patchi];
656 const labelList& own = mesh_.faceOwner();
658 label pointi = pp.meshPoints()[patchPointi];
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),
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,
878 const labelList& featureFaces,
879 const labelList& featureEdges,
880 const labelList& singleCellFeaturePoints,
881 const labelList& multiCellFeaturePoints,
882 polyTopoChange& meshMod
892 bitSet isBoundaryEdge(mesh_.
nEdges());
897 isBoundaryEdge.
set(fEdges);
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
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
953 featureFaceSet[featureFaces[i]] =
true;
962 if (!featureFaceSet[facei])
965 <<
"Not all boundary faces marked as feature faces."
967 <<
"First conflicting face:" << 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
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
1071 pointToDualCells_[pointi].
setSize(pCells.size());
1075 pointToDualCells_[pointi][pCelli] = meshMod.addCell
1088 0.5*(mesh_.
points()[pointi]+cellCentres[pCells[pCelli]])
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
1134 label facei = featureFaces[i];
1136 faceToDualPoint_[facei] = meshMod.addPoint
1139 mesh_.
faces()[facei][0],
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
1178 label edgeI = featureEdges[i];
1180 const edge&
e = mesh_.
edges()[edgeI];
1182 edgeToDualPoint_[edgeI] = meshMod.addPoint
1199 if (edgeToDualPoint_[edgeI] == -1)
1201 const edge&
e = mesh_.
edges()[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
1228 label newDualE1 = findDualCell(eCells[i],
e[1]);
1230 if (dualE1 != newDualE1)
1232 edgeToDualPoint_[edgeI] = meshMod.addPoint
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);
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];
1349 bool foundStart =
false;
1355 edgeToDualPoint_[fEdges[fp]] != -1
1356 && !sameDualCell(facei,
f.nextLabel(fp))
1372 createFaceFromInternalFace
1385 dumpPolyTopoChange(meshMod,
"generatedFacesFromFeatFaces_");
1395 const polyPatch& pp =
patches[patchi];
1399 forAll(pointFaces, patchPointi)
1410 label startFacei = pp.start()+
pFaces[i];
1419 createFacesAroundBoundaryPoint
1434 dumpPolyTopoChange(meshMod,
"generatedFacesFromBndFaces_");
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.
void setSize(const label n)
Alias for resize()
void append(const T &val)
Append an element at the end of the list.
virtual const fileName & name() const
Get the name of the stream.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
label rcIndex(const label i) const noexcept
void size(const label n)
Older name for setAddressableSize.
label fcIndex(const label i) const noexcept
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
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.
Creates dual of polyMesh. Every point becomes a cell (or multiple cells for feature points),...
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.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
virtual const labelList & faceNeighbour() const
Return face neighbour.
const pointZoneMesh & pointZones() const noexcept
Return point zone mesh.
virtual const pointField & points() const
Return raw points.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
const labelListList & pointCells() const
label nInternalFaces() const noexcept
Number of internal faces.
const vectorField & cellCentres() const
label nFaces() const noexcept
Number of mesh faces.
const labelListList & edgeFaces() const
const labelListList & faceEdges() const
const labelListList & edgeCells() const
const cellList & cells() const
label nEdges() const
Number of mesh edges.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const labelList nFaces(UPstream::listGatherValues< label >(aMesh.nFaces()))
const labelList nEdges(UPstream::listGatherValues< label >(aMesh.nEdges()))
const labelIOList & zoneID
Geometric merging of points. See below.
List< label > labelList
A List of labels.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
label mergePoints(const PointList &points, labelList &pointToUnique, labelList &uniquePoints, const scalar mergeTol=SMALL, const bool verbose=false)
Ostream & endl(Ostream &os)
Add newline and flush stream.
List< labelList > labelListList
A List of labelList.
errorManip< error > abort(error &err)
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< bool > boolList
A List of bools.
List< edge > edgeList
A List of edges.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
constexpr char nl
The newline '\n' character (0x0a)
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]
#define forAll(list, i)
Loop across all elements in list.