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()),
140 bool Foam::meshDualiser::sameDualCell
146 if (!mesh_.isInternalFace(facei))
149 <<
"face:" << facei <<
" is not internal face."
153 label own = mesh_.faceOwner()[facei];
154 label nei = mesh_.faceNeighbour()[facei];
156 return findDualCell(own, pointi) == findDualCell(nei, pointi);
162 const label masterPointi,
163 const label masterEdgeI,
164 const label masterFacei,
166 const bool edgeOrder,
167 const label dualCell0,
168 const label dualCell1,
169 const DynamicList<label>& verts,
170 polyTopoChange& meshMod
175 if (edgeOrder != (dualCell0 < dualCell1))
181 bool zoneFlip =
false;
182 if (masterFacei != -1)
184 zoneID = mesh_.faceZones().whichZone(masterFacei);
188 const faceZone& fZone = mesh_.faceZones()[zoneID];
190 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
196 if (dualCell0 < dualCell1)
198 dualFacei = meshMod.addFace
225 dualFacei = meshMod.addFace
256 const label masterPointi,
257 const label masterEdgeI,
258 const label masterFacei,
260 const label dualCelli,
262 const DynamicList<label>& verts,
263 polyTopoChange& meshMod
269 bool zoneFlip =
false;
270 if (masterFacei != -1)
272 zoneID = mesh_.faceZones().whichZone(masterFacei);
276 const faceZone& fZone = mesh_.faceZones()[zoneID];
278 zoneFlip = fZone.flipMap()[fZone.whichFace(masterFacei)];
282 label dualFacei = meshMod.addFace
313 void Foam::meshDualiser::createFacesAroundEdge
315 const bool splitFace,
316 const PackedBoolList& isBoundaryEdge,
318 const label startFacei,
319 polyTopoChange& meshMod,
323 const edge&
e = mesh_.edges()[edgeI];
324 const labelList& eFaces = mesh_.edgeFaces()[edgeI];
328 mesh_.faces()[startFacei],
333 edgeFaceCirculator ie
339 isBoundaryEdge.get(edgeI) == 1
343 bool edgeOrder = ie.sameOrder(
e[0],
e[1]);
344 label startFaceLabel = ie.faceLabel();
355 DynamicList<label> verts(100);
357 if (edgeToDualPoint_[edgeI] != -1)
359 verts.append(edgeToDualPoint_[edgeI]);
361 if (faceToDualPoint_[ie.faceLabel()] != -1)
363 doneEFaces[
findIndex(eFaces, ie.faceLabel())] =
true;
364 verts.append(faceToDualPoint_[ie.faceLabel()]);
366 if (cellToDualPoint_[ie.cellLabel()] != -1)
368 verts.append(cellToDualPoint_[ie.cellLabel()]);
371 label currentDualCell0 = findDualCell(ie.cellLabel(),
e[0]);
372 label currentDualCell1 = findDualCell(ie.cellLabel(),
e[1]);
378 label facei = ie.faceLabel();
381 doneEFaces[
findIndex(eFaces, facei)] =
true;
383 if (faceToDualPoint_[facei] != -1)
385 verts.append(faceToDualPoint_[facei]);
388 label celli = ie.cellLabel();
398 label dualCell0 = findDualCell(celli,
e[0]);
399 label dualCell1 = findDualCell(celli,
e[1]);
407 dualCell0 != currentDualCell0
408 || dualCell1 != currentDualCell1
426 currentDualCell0 = dualCell0;
427 currentDualCell1 = dualCell1;
430 if (edgeToDualPoint_[edgeI] != -1)
432 verts.append(edgeToDualPoint_[edgeI]);
434 if (faceToDualPoint_[facei] != -1)
436 verts.append(faceToDualPoint_[facei]);
440 if (cellToDualPoint_[celli] != -1)
442 verts.append(cellToDualPoint_[celli]);
451 if (isBoundaryEdge.get(edgeI) == 0)
453 label startDual = faceToDualPoint_[startFaceLabel];
455 if (startDual != -1 &&
findIndex(verts, startDual) == -1)
457 verts.append(startDual);
482 void Foam::meshDualiser::createFaceFromInternalFace
486 polyTopoChange& meshMod
489 const face&
f = mesh_.faces()[facei];
490 const labelList& fEdges = mesh_.faceEdges()[facei];
491 label own = mesh_.faceOwner()[facei];
492 label nei = mesh_.faceNeighbour()[facei];
503 DynamicList<label> verts(100);
505 verts.append(faceToDualPoint_[facei]);
506 verts.append(edgeToDualPoint_[fEdges[fp]]);
512 label currentDualCell0 = findDualCell(own,
f[fp]);
513 label currentDualCell1 = findDualCell(nei,
f[fp]);
518 if (pointToDualPoint_[
f[fp]] != -1)
520 verts.append(pointToDualPoint_[
f[fp]]);
524 label edgeI = fEdges[fp];
526 if (edgeToDualPoint_[edgeI] != -1)
528 verts.append(edgeToDualPoint_[edgeI]);
535 label dualCell0 = findDualCell(own,
f[nextFp]);
536 label dualCell1 = findDualCell(nei,
f[nextFp]);
538 if (dualCell0 != currentDualCell0 || dualCell1 != currentDualCell1)
541 if (edgeToDualPoint_[edgeI] == -1)
544 <<
"face:" << facei <<
" verts:" <<
f
545 <<
" points:" << UIndirectList<point>(mesh_.points(),
f)()
546 <<
" no feature edge between " <<
f[fp]
547 <<
" and " <<
f[nextFp] <<
" although have different"
548 <<
" dual cells." <<
endl
549 <<
"point " <<
f[fp] <<
" has dual cells "
550 << currentDualCell0 <<
" and " << currentDualCell1
551 <<
" ; point "<<
f[nextFp] <<
" has dual cells "
552 << dualCell0 <<
" and " << dualCell1
580 void Foam::meshDualiser::createFacesAroundBoundaryPoint
583 const label patchPointi,
584 const label startFacei,
585 polyTopoChange& meshMod,
589 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
592 const labelList& own = mesh_.faceOwner();
594 label pointi = pp.meshPoints()[patchPointi];
596 if (pointToDualPoint_[pointi] == -1)
602 label facei = startFacei;
604 DynamicList<label> verts(4);
611 if (donePFaces[index])
615 donePFaces[index] =
true;
618 verts.append(faceToDualPoint_[facei]);
620 label dualCelli = findDualCell(own[facei], pointi);
623 const face&
f = mesh_.faces()[facei];
626 label edgeI = mesh_.faceEdges()[facei][prevFp];
628 if (edgeToDualPoint_[edgeI] != -1)
630 verts.append(edgeToDualPoint_[edgeI]);
634 edgeFaceCirculator circ
647 while (mesh_.isInternalFace(circ.faceLabel()));
650 facei = circ.faceLabel();
652 if (facei < pp.start() || facei >= pp.start()+pp.size())
655 <<
"Walked from face on patch:" <<
patchi
656 <<
" to face:" << facei
657 <<
" fc:" << mesh_.faceCentres()[facei]
658 <<
" on patch:" <<
patches.whichPatch(facei)
663 if (dualCelli != findDualCell(own[facei], pointi))
666 <<
"Different dual cells but no feature edge"
667 <<
" in between point:" << pointi
668 <<
" coord:" << mesh_.points()[pointi]
675 label dualCelli = findDualCell(own[facei], pointi);
695 label facei = startFacei;
698 DynamicList<label> verts(mesh_.faces()[facei].size());
701 verts.append(pointToDualPoint_[pointi]);
704 const labelList& fEdges = mesh_.faceEdges()[facei];
705 label nextEdgeI = fEdges[
findIndex(mesh_.faces()[facei], pointi)];
706 if (edgeToDualPoint_[nextEdgeI] != -1)
708 verts.append(edgeToDualPoint_[nextEdgeI]);
716 if (donePFaces[index])
720 donePFaces[index] =
true;
723 verts.append(faceToDualPoint_[facei]);
726 const labelList& fEdges = mesh_.faceEdges()[facei];
727 const face&
f = mesh_.faces()[facei];
729 label edgeI = fEdges[prevFp];
731 if (edgeToDualPoint_[edgeI] != -1)
735 verts.append(edgeToDualPoint_[edgeI]);
741 findDualCell(own[facei], pointi),
748 verts.append(pointToDualPoint_[pointi]);
749 verts.append(edgeToDualPoint_[edgeI]);
753 edgeFaceCirculator circ
766 while (mesh_.isInternalFace(circ.faceLabel()));
769 facei = circ.faceLabel();
774 && facei >= pp.start()
775 && facei < pp.start()+pp.size()
778 if (verts.size() > 2)
786 findDualCell(own[facei], pointi),
801 pointToDualCells_(mesh_.
nPoints()),
802 pointToDualPoint_(mesh_.
nPoints(), -1),
803 cellToDualPoint_(mesh_.nCells()),
804 faceToDualPoint_(mesh_.nFaces(), -1),
805 edgeToDualPoint_(mesh_.nEdges(), -1)
813 const bool splitFace,
816 const labelList& singleCellFeaturePoints,
818 polyTopoChange& meshMod
828 PackedBoolList isBoundaryEdge(mesh_.
nEdges());
835 isBoundaryEdge.set(fEdges[i], 1);
852 featureFaceSet[featureFaces[i]] =
true;
859 <<
"In split-face-mode (splitFace=true) but not all faces"
860 <<
" marked as feature faces." <<
endl
861 <<
"First conflicting face:" << facei
869 featureEdgeSet[featureEdges[i]] =
true;
875 const edge&
e = mesh_.
edges()[edgeI];
877 <<
"In split-face-mode (splitFace=true) but not all edges"
878 <<
" marked as feature edges." <<
endl
879 <<
"First conflicting edge:" << edgeI
892 featureFaceSet[featureFaces[i]] =
true;
901 if (!featureFaceSet[facei])
904 <<
"Not all boundary faces marked as feature faces."
906 <<
"First conflicting face:" << facei
937 autoPtr<OFstream> dualCcStr;
940 dualCcStr.reset(
new OFstream(
"dualCc.obj"));
941 Pout<<
"Dumping centres of dual cells to " << dualCcStr().
name()
956 forAll(singleCellFeaturePoints, i)
958 label pointi = singleCellFeaturePoints[i];
960 pointToDualPoint_[pointi] = meshMod.addPoint
969 pointToDualCells_[pointi].
setSize(1);
970 pointToDualCells_[pointi][0] = meshMod.addCell
978 if (dualCcStr.valid())
985 forAll(multiCellFeaturePoints, i)
987 label pointi = multiCellFeaturePoints[i];
989 if (pointToDualCells_[pointi].size() > 0)
992 <<
"Point " << pointi <<
" at:" << mesh_.
points()[pointi]
993 <<
" is both in singleCellFeaturePoints"
994 <<
" and multiCellFeaturePoints."
998 pointToDualPoint_[pointi] = meshMod.addPoint
1010 pointToDualCells_[pointi].
setSize(pCells.size());
1014 pointToDualCells_[pointi][pCelli] = meshMod.addCell
1022 if (dualCcStr.valid())
1027 0.5*(mesh_.
points()[pointi]+cellCentres[pCells[pCelli]])
1035 if (pointToDualCells_[pointi].empty())
1037 pointToDualCells_[pointi].
setSize(1);
1038 pointToDualCells_[pointi][0] = meshMod.addCell
1047 if (dualCcStr.valid())
1058 forAll(cellToDualPoint_, celli)
1060 cellToDualPoint_[celli] = meshMod.addPoint
1073 label facei = featureFaces[i];
1075 faceToDualPoint_[facei] = meshMod.addPoint
1078 mesh_.
faces()[facei][0],
1088 if (faceToDualPoint_[facei] == -1)
1090 const face&
f = mesh_.
faces()[facei];
1094 label ownDualCell = findDualCell(own[facei],
f[fp]);
1095 label neiDualCell = findDualCell(nei[facei],
f[fp]);
1097 if (ownDualCell != neiDualCell)
1099 faceToDualPoint_[facei] = meshMod.addPoint
1117 label edgeI = featureEdges[i];
1119 const edge&
e = mesh_.
edges()[edgeI];
1121 edgeToDualPoint_[edgeI] = meshMod.addPoint
1138 if (edgeToDualPoint_[edgeI] == -1)
1140 const edge&
e = mesh_.
edges()[edgeI];
1147 label dualE0 = findDualCell(eCells[0],
e[0]);
1148 label dualE1 = findDualCell(eCells[0],
e[1]);
1150 for (
label i = 1; i < eCells.size(); i++)
1152 label newDualE0 = findDualCell(eCells[i],
e[0]);
1154 if (dualE0 != newDualE0)
1156 edgeToDualPoint_[edgeI] = meshMod.addPoint
1167 label newDualE1 = findDualCell(eCells[i],
e[1]);
1169 if (dualE1 != newDualE1)
1171 edgeToDualPoint_[edgeI] = meshMod.addPoint
1187 forAll(singleCellFeaturePoints, i)
1189 generateDualBoundaryEdges
1192 singleCellFeaturePoints[i],
1196 forAll(multiCellFeaturePoints, i)
1198 generateDualBoundaryEdges
1201 multiCellFeaturePoints[i],
1210 dumpPolyTopoChange(meshMod,
"generatedPoints_");
1230 boolList doneEFaces(eFaces.size(),
false);
1240 label startFacei = eFaces[i];
1249 createFacesAroundEdge
1264 dumpPolyTopoChange(meshMod,
"generatedFacesFromEdges_");
1273 forAll(faceToDualPoint_, facei)
1275 if (faceToDualPoint_[facei] != -1 && mesh_.
isInternalFace(facei))
1277 const face&
f = mesh_.
faces()[facei];
1287 bool foundStart =
false;
1293 edgeToDualPoint_[fEdges[fp]] != -1
1294 && !sameDualCell(facei,
f.nextLabel(fp))
1310 createFaceFromInternalFace
1323 dumpPolyTopoChange(meshMod,
"generatedFacesFromFeatFaces_");
1337 forAll(pointFaces, patchPointi)
1357 createFacesAroundBoundaryPoint
1372 dumpPolyTopoChange(meshMod,
"generatedFacesFromBndFaces_");
#define forAll(list, i)
Loop across all elements in list.
void setSize(const label)
Reset size of List.
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
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 meshPointZones & pointZones() const
Return point zones.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const labelList & faceNeighbour() const
Return face neighbour.
virtual const pointField & points() const
Return raw points.
const meshCellZones & cellZones() const
Return cell zones.
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]