45 void Foam::meshDualiser::dumpPolyTopoChange
47 const polyTopoChange& meshMod,
48 const fileName& prefix
51 OFstream str1(prefix +
"Faces.obj");
52 OFstream str2(prefix +
"Edges.obj");
54 Info<<
"Dumping current polyTopoChange. Faces to " << str1.name()
55 <<
" , points and edges to " << str2.name() <<
endl;
57 const DynamicList<point>&
points = meshMod.points();
65 const DynamicList<face>& faces = meshMod.faces();
69 const face&
f = faces[facei];
74 str1<<
' ' <<
f[fp]+1;
81 str2<<
' ' <<
f[fp]+1;
83 str2<<
' ' <<
f[0]+1 <<
nl;
94 const labelList& dualCells = pointToDualCells_[pointi];
96 if (dualCells.size() == 1)
104 return dualCells[index];
109 void Foam::meshDualiser::generateDualBoundaryEdges
111 const PackedBoolList& isBoundaryEdge,
113 polyTopoChange& meshMod
116 const labelList& pEdges = mesh_.pointEdges()[pointi];
120 label edgeI = pEdges[pEdgeI];
122 if (edgeToDualPoint_[edgeI] == -1 && isBoundaryEdge.get(edgeI) == 1)
124 const edge&
e = mesh_.edges()[edgeI];
126 edgeToDualPoint_[edgeI] = meshMod.addPoint
128 e.centre(mesh_.points()),
139 bool Foam::meshDualiser::sameDualCell
145 if (!mesh_.isInternalFace(facei))
148 <<
"face:" << facei <<
" is not internal face."
152 label own = mesh_.faceOwner()[facei];
153 label nei = mesh_.faceNeighbour()[facei];
155 return findDualCell(own, pointi) == findDualCell(nei, pointi);
161 const label masterPointi,
162 const label masterEdgeI,
163 const label masterFacei,
165 const bool edgeOrder,
166 const label dualCell0,
167 const label dualCell1,
168 const DynamicList<label>& verts,
169 polyTopoChange& meshMod
174 if (edgeOrder != (dualCell0 < dualCell1))
181 if (dualCell0 < dualCell1)
183 dualFacei = meshMod.addFace
206 dualFacei = meshMod.addFace
233 const label masterPointi,
234 const label masterEdgeI,
235 const label masterFacei,
237 const label dualCelli,
239 const DynamicList<label>& verts,
240 polyTopoChange& meshMod
245 label dualFacei = meshMod.addFace
272 void Foam::meshDualiser::createFacesAroundEdge
274 const bool splitFace,
275 const PackedBoolList& isBoundaryEdge,
277 const label startFacei,
278 polyTopoChange& meshMod,
282 const edge&
e = mesh_.edges()[edgeI];
283 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
287 mesh_.faces()[startFacei],
292 edgeFaceCirculator ie
298 isBoundaryEdge.get(edgeI) == 1
302 bool edgeOrder = ie.sameOrder(
e[0],
e[1]);
303 label startFaceLabel = ie.faceLabel();
314 DynamicList<label> verts(100);
316 if (edgeToDualPoint_[edgeI] != -1)
318 verts.append(edgeToDualPoint_[edgeI]);
320 if (faceToDualPoint_[ie.faceLabel()] != -1)
322 doneEFaces[
findIndex(eFaces, ie.faceLabel())] =
true;
323 verts.append(faceToDualPoint_[ie.faceLabel()]);
325 if (cellToDualPoint_[ie.cellLabel()] != -1)
327 verts.append(cellToDualPoint_[ie.cellLabel()]);
330 label currentDualCell0 = findDualCell(ie.cellLabel(),
e[0]);
331 label currentDualCell1 = findDualCell(ie.cellLabel(),
e[1]);
337 label facei = ie.faceLabel();
340 doneEFaces[
findIndex(eFaces, facei)] =
true;
342 if (faceToDualPoint_[facei] != -1)
344 verts.append(faceToDualPoint_[facei]);
347 label celli = ie.cellLabel();
357 label dualCell0 = findDualCell(celli,
e[0]);
358 label dualCell1 = findDualCell(celli,
e[1]);
366 dualCell0 != currentDualCell0
367 || dualCell1 != currentDualCell1
385 currentDualCell0 = dualCell0;
386 currentDualCell1 = dualCell1;
389 if (edgeToDualPoint_[edgeI] != -1)
391 verts.append(edgeToDualPoint_[edgeI]);
393 if (faceToDualPoint_[facei] != -1)
395 verts.append(faceToDualPoint_[facei]);
399 if (cellToDualPoint_[celli] != -1)
401 verts.append(cellToDualPoint_[celli]);
410 if (isBoundaryEdge.get(edgeI) == 0)
412 label startDual = faceToDualPoint_[startFaceLabel];
414 if (startDual != -1 &&
findIndex(verts, startDual) == -1)
416 verts.append(startDual);
441 void Foam::meshDualiser::createFaceFromInternalFace
445 polyTopoChange& meshMod
448 const face&
f = mesh_.faces()[facei];
449 const labelList& fEdges = mesh_.faceEdges()[facei];
450 label own = mesh_.faceOwner()[facei];
451 label nei = mesh_.faceNeighbour()[facei];
462 DynamicList<label> verts(100);
464 verts.append(faceToDualPoint_[facei]);
465 verts.append(edgeToDualPoint_[fEdges[fp]]);
471 label currentDualCell0 = findDualCell(own,
f[fp]);
472 label currentDualCell1 = findDualCell(nei,
f[fp]);
477 if (pointToDualPoint_[
f[fp]] != -1)
479 verts.append(pointToDualPoint_[
f[fp]]);
483 label edgeI = fEdges[fp];
485 if (edgeToDualPoint_[edgeI] != -1)
487 verts.append(edgeToDualPoint_[edgeI]);
494 label dualCell0 = findDualCell(own,
f[nextFp]);
495 label dualCell1 = findDualCell(nei,
f[nextFp]);
497 if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
500 if (edgeToDualPoint_[edgeI] == -1)
503 <<
"face:" << facei <<
" verts:" <<
f
504 <<
" points:" << UIndirectList<point>(mesh_.points(),
f)()
505 <<
" no feature edge between " <<
f[fp]
506 <<
" and " <<
f[nextFp] <<
" although have different"
507 <<
" dual cells." <<
endl
508 <<
"point " <<
f[fp] <<
" has dual cells "
509 << currentDualCell0 <<
" and " << currentDualCell1
510 <<
" ; point "<<
f[nextFp] <<
" has dual cells "
511 << dualCell0 <<
" and " << dualCell1
539 void Foam::meshDualiser::createFacesAroundBoundaryPoint
542 const label patchPointi,
543 const label startFacei,
544 polyTopoChange& meshMod,
548 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
551 const labelList& own = mesh_.faceOwner();
553 label pointi = pp.meshPoints()[patchPointi];
555 if (pointToDualPoint_[pointi] == -1)
561 label facei = startFacei;
563 DynamicList<label> verts(4);
570 if (donePFaces[index])
574 donePFaces[index] =
true;
577 verts.append(faceToDualPoint_[facei]);
579 label dualCelli = findDualCell(own[facei], pointi);
582 const face&
f = mesh_.faces()[facei];
585 label edgeI = mesh_.faceEdges()[facei][prevFp];
587 if (edgeToDualPoint_[edgeI] != -1)
589 verts.append(edgeToDualPoint_[edgeI]);
593 edgeFaceCirculator circ
606 while (mesh_.isInternalFace(circ.faceLabel()));
609 facei = circ.faceLabel();
611 if (facei < pp.start() || facei >= pp.start()+pp.size())
614 <<
"Walked from face on patch:" <<
patchi
615 <<
" to face:" << facei
616 <<
" fc:" << mesh_.faceCentres()[facei]
617 <<
" on patch:" <<
patches.whichPatch(facei)
622 if (dualCelli != findDualCell(own[facei], pointi))
625 <<
"Different dual cells but no feature edge"
626 <<
" in between point:" << pointi
627 <<
" coord:" << mesh_.points()[pointi]
634 label dualCelli = findDualCell(own[facei], pointi);
654 label facei = startFacei;
657 DynamicList<label> verts(mesh_.faces()[facei].size());
660 verts.append(pointToDualPoint_[pointi]);
663 const labelList& fEdges = mesh_.faceEdges()[facei];
664 label nextEdgeI = fEdges[
findIndex(mesh_.faces()[facei], pointi)];
665 if (edgeToDualPoint_[nextEdgeI] != -1)
667 verts.append(edgeToDualPoint_[nextEdgeI]);
675 if (donePFaces[index])
679 donePFaces[index] =
true;
682 verts.append(faceToDualPoint_[facei]);
685 const labelList& fEdges = mesh_.faceEdges()[facei];
686 const face&
f = mesh_.faces()[facei];
688 label edgeI = fEdges[prevFp];
690 if (edgeToDualPoint_[edgeI] != -1)
694 verts.append(edgeToDualPoint_[edgeI]);
700 findDualCell(own[facei], pointi),
707 verts.append(pointToDualPoint_[pointi]);
708 verts.append(edgeToDualPoint_[edgeI]);
712 edgeFaceCirculator circ
725 while (mesh_.isInternalFace(circ.faceLabel()));
728 facei = circ.faceLabel();
733 && facei >= pp.start()
734 && facei < pp.start()+pp.size()
737 if (verts.size() > 2)
745 findDualCell(own[facei], pointi),
760 pointToDualCells_(mesh_.
nPoints()),
761 pointToDualPoint_(mesh_.
nPoints(), -1),
762 cellToDualPoint_(mesh_.nCells()),
763 faceToDualPoint_(mesh_.nFaces(), -1),
764 edgeToDualPoint_(mesh_.nEdges(), -1)
772 const bool splitFace,
775 const labelList& singleCellFeaturePoints,
777 polyTopoChange& meshMod
787 PackedBoolList isBoundaryEdge(mesh_.
nEdges());
794 isBoundaryEdge.set(fEdges[i], 1);
811 featureFaceSet[featureFaces[i]] =
true;
818 <<
"In split-face-mode (splitFace=true) but not all faces"
819 <<
" marked as feature faces." <<
endl
820 <<
"First conflicting face:" << facei
828 featureEdgeSet[featureEdges[i]] =
true;
834 const edge&
e = mesh_.
edges()[edgeI];
836 <<
"In split-face-mode (splitFace=true) but not all edges"
837 <<
" marked as feature edges." <<
endl
838 <<
"First conflicting edge:" << edgeI
851 featureFaceSet[featureFaces[i]] =
true;
860 if (!featureFaceSet[facei])
863 <<
"Not all boundary faces marked as feature faces."
865 <<
"First conflicting face:" << facei
896 autoPtr<OFstream> dualCcStr;
899 dualCcStr.reset(
new OFstream(
"dualCc.obj"));
900 Pout<<
"Dumping centres of dual cells to " << dualCcStr().
name()
915 forAll(singleCellFeaturePoints, i)
917 label pointi = singleCellFeaturePoints[i];
919 pointToDualPoint_[pointi] = meshMod.addPoint
927 pointToDualCells_[pointi].
setSize(1);
928 pointToDualCells_[pointi][0] = meshMod.addCell
932 if (dualCcStr.valid())
939 forAll(multiCellFeaturePoints, i)
941 label pointi = multiCellFeaturePoints[i];
943 if (pointToDualCells_[pointi].size() > 0)
946 <<
"Point " << pointi <<
" at:" << mesh_.
points()[pointi]
947 <<
" is both in singleCellFeaturePoints"
948 <<
" and multiCellFeaturePoints."
952 pointToDualPoint_[pointi] = meshMod.addPoint
963 pointToDualCells_[pointi].
setSize(pCells.size());
967 pointToDualCells_[pointi][pCelli] = meshMod.addCell
971 if (dualCcStr.valid())
976 0.5*(mesh_.
points()[pointi]+cellCentres[pCells[pCelli]])
984 if (pointToDualCells_[pointi].empty())
986 pointToDualCells_[pointi].
setSize(1);
987 pointToDualCells_[pointi][0] = meshMod.addCell
992 if (dualCcStr.valid())
1003 forAll(cellToDualPoint_, celli)
1005 cellToDualPoint_[celli] = meshMod.addPoint
1017 label facei = featureFaces[i];
1019 faceToDualPoint_[facei] = meshMod.addPoint
1022 mesh_.
faces()[facei][0],
1031 if (faceToDualPoint_[facei] == -1)
1033 const face&
f = mesh_.
faces()[facei];
1037 label ownDualCell = findDualCell(own[facei],
f[fp]);
1038 label neiDualCell = findDualCell(nei[facei],
f[fp]);
1040 if (ownDualCell != neiDualCell)
1042 faceToDualPoint_[facei] = meshMod.addPoint
1059 label edgeI = featureEdges[i];
1061 const edge&
e = mesh_.
edges()[edgeI];
1063 edgeToDualPoint_[edgeI] = meshMod.addPoint
1079 if (edgeToDualPoint_[edgeI] == -1)
1081 const edge&
e = mesh_.
edges()[edgeI];
1088 label dualE0 = findDualCell(eCells[0],
e[0]);
1089 label dualE1 = findDualCell(eCells[0],
e[1]);
1091 for (
label i = 1; i < eCells.size(); i++)
1093 label newDualE0 = findDualCell(eCells[i],
e[0]);
1095 if (dualE0 != newDualE0)
1097 edgeToDualPoint_[edgeI] = meshMod.addPoint
1107 label newDualE1 = findDualCell(eCells[i],
e[1]);
1109 if (dualE1 != newDualE1)
1111 edgeToDualPoint_[edgeI] = meshMod.addPoint
1126 forAll(singleCellFeaturePoints, i)
1128 generateDualBoundaryEdges
1131 singleCellFeaturePoints[i],
1135 forAll(multiCellFeaturePoints, i)
1137 generateDualBoundaryEdges
1140 multiCellFeaturePoints[i],
1149 dumpPolyTopoChange(meshMod,
"generatedPoints_");
1169 boolList doneEFaces(eFaces.size(),
false);
1179 label startFacei = eFaces[i];
1188 createFacesAroundEdge
1203 dumpPolyTopoChange(meshMod,
"generatedFacesFromEdges_");
1212 forAll(faceToDualPoint_, facei)
1214 if (faceToDualPoint_[facei] != -1 && mesh_.
isInternalFace(facei))
1216 const face&
f = mesh_.
faces()[facei];
1226 bool foundStart =
false;
1232 edgeToDualPoint_[fEdges[fp]] != -1
1233 && !sameDualCell(facei,
f.nextLabel(fp))
1249 createFaceFromInternalFace
1262 dumpPolyTopoChange(meshMod,
"generatedFacesFromFeatFaces_");
1276 forAll(pointFaces, patchPointi)
1296 createFacesAroundBoundaryPoint
1311 dumpPolyTopoChange(meshMod,
"generatedFacesFromBndFaces_");
#define forAll(list, i)
Loop across all elements in list.
void setSize(const label)
Reset size of List.
virtual const fileName & name() const
Return the name of the stream.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
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.
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.
meshDualiser(const polyMesh &)
Construct from mesh.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const labelList & faceNeighbour() const
Return face neighbour.
virtual const pointField & points() const
Return raw points.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
const labelListList & pointCells() const
const vectorField & cellCentres() const
label nInternalFaces() const
const labelListList & edgeFaces() const
const labelListList & faceEdges() const
const labelListList & edgeCells() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const cellList & cells() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const fvPatchList & patches
List< label > labelList
A List of labels.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Ostream & endl(Ostream &os)
Add newline and flush stream.
errorManip< error > abort(error &err)
List< bool > boolList
Bool container classes.
void reverse(UList< T > &, const label n)
List< labelList > labelListList
A List of labelList.
defineTypeNameAndDebug(combustionModel, 0)
Field< vector > vectorField
Specialisation of Field<T> for vector.
prefixOSstream Pout(cout, "Pout")
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("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]=cellShape(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, &mergedCyclicPolyPatch::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]