69 SortableList<label> nbr(cFaces.size());
73 label facei = cFaces[i];
80 if (nbrCelli == celli)
109 oldToNew[cFaces[nbr.indices()[i]]] = newFacei++;
121 for (label facei = newFacei; facei <
faceOwner.
size(); facei++)
123 oldToNew[facei] = facei;
130 if (oldToNew[facei] == -1)
133 <<
"Did not determine new position" 134 <<
" for face " << facei
148 void Foam::polyDualMesh::getPointEdges
158 const face&
f =
patch.localFaces()[facei];
165 label edgeI = fEdges[i];
167 const edge&
e =
patch.edges()[edgeI];
172 label index =
f.
find(pointi);
174 if (
f.nextLabel(index) ==
e[1])
183 if (e0 != -1 && e1 != -1)
188 else if (
e[1] == pointi)
191 label index =
f.
find(pointi);
193 if (
f.nextLabel(index) ==
e[0])
202 if (e0 != -1 && e1 != -1)
210 <<
" vertices:" <<
patch.localFaces()[facei]
211 <<
" that uses point:" << pointi
219 const polyPatch&
patch,
220 const label patchToDualOffset,
231 label meshPointi =
patch.meshPoints()[pointi];
234 DynamicList<label> dualFace;
236 if (pointToDualPoint[meshPointi] >= 0)
239 dualFace.setCapacity(
pFaces.size()+2+1);
241 dualFace.append(pointToDualPoint[meshPointi]);
245 dualFace.setCapacity(
pFaces.size()+2);
249 if (edgeToDualPoint[
patch.meshEdges()[edgeI]] < 0)
255 dualFace.append(edgeToDualPoint[
patch.meshEdges()[edgeI]]);
257 label facei =
patch.edgeFaces()[edgeI][0];
263 getPointEdges(
patch, facei, pointi, e0, e1);
277 dualFace.append(facei + patchToDualOffset);
281 getPointEdges(
patch, facei, pointi, e0, e1);
292 if (edgeToDualPoint[
patch.meshEdges()[edgeI]] >= 0)
295 dualFace.append(edgeToDualPoint[
patch.meshEdges()[edgeI]]);
300 if (eFaces.size() == 1)
307 if (eFaces[0] == facei)
331 void Foam::polyDualMesh::collectPatchInternalFace
333 const polyPatch&
patch,
334 const label patchToDualOffset,
338 const label startEdgeI,
349 DynamicList<label> dualFace(
pFaces.size());
351 DynamicList<label> featEdgeIndices(dualFace.size());
354 label edgeI = startEdgeI;
355 label facei =
patch.edgeFaces()[edgeI][0];
361 getPointEdges(
patch, facei, pointi, e0, e1);
375 dualFace.append(facei + patchToDualOffset);
379 getPointEdges(
patch, facei, pointi, e0, e1);
390 if (edgeToDualPoint[meshEdges[edgeI]] >= 0)
393 dualFace.append(edgeToDualPoint[meshEdges[edgeI]]);
394 featEdgeIndices.append(dualFace.size()-1);
397 if (edgeI == startEdgeI)
405 if (eFaces[0] == facei)
415 dualFace2.transfer(dualFace);
417 featEdgeIndices2.transfer(featEdgeIndices);
424 forAll(featEdgeIndices2, i)
426 featEdgeIndices2[i] = dualFace2.size() -1 - featEdgeIndices2[i];
434 void Foam::polyDualMesh::splitFace
436 const polyPatch&
patch,
443 DynamicList<face>& dualFaces,
444 DynamicList<label>& dualOwner,
445 DynamicList<label>& dualNeighbour,
446 DynamicList<label>& dualRegion
452 label meshPointi =
patch.meshPoints()[pointi];
454 if (pointToDualPoint[meshPointi] >= 0)
458 if (featEdgeIndices.size() < 2)
461 dualFaces.append(face(dualFace));
462 dualOwner.append(meshPointi);
463 dualNeighbour.append(-1);
464 dualRegion.append(
patch.index());
471 forAll(featEdgeIndices, i)
473 label startFp = featEdgeIndices[i];
475 label endFp = featEdgeIndices[(i+1) % featEdgeIndices.size()];
482 sz = endFp - startFp + 2;
486 sz = endFp + dualFace.size() - startFp + 2;
491 subFace[0] = pointToDualPoint[
patch.meshPoints()[pointi]];
494 for (label subFp = 1; subFp < subFace.size(); subFp++)
496 subFace[subFp] = dualFace[startFp];
498 startFp = (startFp+1) % dualFace.size();
501 dualFaces.append(face(subFace));
502 dualOwner.append(meshPointi);
503 dualNeighbour.append(-1);
504 dualRegion.append(
patch.index());
511 if (featEdgeIndices.size() < 2)
514 dualFaces.append(face(dualFace));
515 dualOwner.append(meshPointi);
516 dualNeighbour.append(-1);
517 dualRegion.append(
patch.index());
526 DynamicList<label> subFace(dualFace.size());
528 forAll(featEdgeIndices, featI)
530 label startFp = featEdgeIndices[featI];
531 label endFp = featEdgeIndices[featEdgeIndices.fcIndex(featI)];
537 label vertI = dualFace[fp];
539 subFace.append(vertI);
546 fp = dualFace.fcIndex(fp);
549 if (subFace.size() > 2)
552 dualFaces.append(face(subFace));
553 dualOwner.append(meshPointi);
554 dualNeighbour.append(-1);
555 dualRegion.append(
patch.index());
561 if (subFace.size() > 2)
564 dualFaces.append(face(subFace));
565 dualOwner.append(meshPointi);
566 dualNeighbour.append(-1);
567 dualRegion.append(
patch.index());
577 void Foam::polyDualMesh::dualPatch
579 const polyPatch&
patch,
580 const label patchToDualOffset,
586 DynamicList<face>& dualFaces,
587 DynamicList<label>& dualOwner,
588 DynamicList<label>& dualNeighbour,
589 DynamicList<label>& dualRegion
601 bitSet donePoint(
patch.nPoints(),
false);
607 forAll(doneEdgeSide, patchEdgeI)
611 if (eFaces.size() == 1)
613 const edge&
e =
patch.edges()[patchEdgeI];
617 label bitMask = 1<<eI;
619 if ((doneEdgeSide[patchEdgeI] & bitMask) == 0)
624 label pointi =
e[eI];
626 label edgeI = patchEdgeI;
642 if (
patch.edges()[edgeI][0] == pointi)
644 doneEdgeSide[edgeI] |= 1;
648 doneEdgeSide[edgeI] |= 2;
651 dualFaces.append(face(dualFace));
652 dualOwner.append(
patch.meshPoints()[pointi]);
653 dualNeighbour.append(-1);
654 dualRegion.append(
patch.index());
656 doneEdgeSide[patchEdgeI] |= bitMask;
657 donePoint.set(pointi);
670 if (!donePoint.test(pointi))
674 collectPatchInternalFace
680 patch.pointEdges()[pointi][0],
708 donePoint[pointi] =
true;
714 void Foam::polyDualMesh::calcDual
716 const polyMesh&
mesh,
740 allBoundary.meshEdges
748 pointSet nonManifoldPoints
752 meshEdges.
size() / 100
755 allBoundary.checkPointManifold(
true, &nonManifoldPoints);
757 if (nonManifoldPoints.size())
759 nonManifoldPoints.write();
762 <<
"There are " << nonManifoldPoints.size() <<
" points where" 763 <<
" the outside of the mesh is non-manifold." <<
nl 764 <<
"Such a mesh cannot be converted to a dual." <<
nl 765 <<
"Writing points at non-manifold sites to pointSet " 766 << nonManifoldPoints.name()
786 + featureEdges.size()
787 + featurePoints.size()
790 label dualPointi = 0;
796 cellPoint_.
setSize(cellCentres.size());
798 forAll(cellCentres, celli)
800 cellPoint_[celli] = dualPointi;
801 dualPoints[dualPointi++] = cellCentres[celli];
809 for (label facei = nIntFaces; facei <
mesh.
nFaces(); facei++)
811 boundaryFacePoint_[facei - nIntFaces] = dualPointi;
812 dualPoints[dualPointi++] = faceCentres[facei];
821 forAll(meshEdges, patchEdgeI)
823 label edgeI = meshEdges[patchEdgeI];
825 edgeToDualPoint[edgeI] = -1;
830 label edgeI = featureEdges[i];
834 edgeToDualPoint[edgeI] = dualPointi;
835 dualPoints[dualPointi++] =
e.centre(
mesh.
points());
853 pointToDualPoint[meshPoints[i]] = -2;
864 pointToDualPoint[meshPoints[loop[j]]] = -1;
871 label pointi = featurePoints[i];
873 pointToDualPoint[pointi] = dualPointi;
874 dualPoints[dualPointi++] =
mesh.
points()[pointi];
881 DynamicList<face> dynDualFaces(
mesh.
nEdges());
882 DynamicList<label> dynDualOwner(
mesh.
nEdges());
883 DynamicList<label> dynDualNeighbour(
mesh.
nEdges());
884 DynamicList<label> dynDualRegion(
mesh.
nEdges());
890 forAll(meshEdges, patchEdgeI)
892 label edgeI = meshEdges[patchEdgeI];
897 label neighbour = -1;
911 const labelList& patchFaces = allBoundary.edgeFaces()[patchEdgeI];
913 if (patchFaces.size() != 2)
916 <<
"Cannot handle edges with " << patchFaces.size()
917 <<
" connected boundary faces." 921 label face0 = patchFaces[0] + nIntFaces;
924 label face1 = patchFaces[1] + nIntFaces;
929 label startFacei = -1;
932 label index = f0.
find(neighbour);
934 if (f0.nextLabel(index) == owner)
947 DynamicList<label> dualFace;
949 if (edgeToDualPoint[edgeI] >= 0)
954 dualFace.append(edgeToDualPoint[edgeI]);
962 dualFace.append(
mesh.
nCells() + startFacei - nIntFaces);
965 label facei = startFacei;
986 if (facei == endFacei)
1004 dynDualFaces.append(face(dualFace.shrink()));
1005 dynDualOwner.append(owner);
1006 dynDualNeighbour.append(neighbour);
1007 dynDualRegion.append(-1);
1011 const face&
f = dynDualFaces.
last();
1012 const vector areaNorm =
f.areaNormal(dualPoints);
1014 if (((
mesh.
points()[owner] - dualPoints[
f[0]]) & areaNorm) > 0)
1017 <<
" on boundary edge:" << edgeI
1029 forAll(edgeToDualPoint, edgeI)
1031 if (edgeToDualPoint[edgeI] == -2)
1040 label neighbour = -1;
1056 label celli = eCells[0];
1066 label index = f0.
find(neighbour);
1068 bool f0OrderOk = (f0.nextLabel(index) == owner);
1070 label startFacei = -1;
1084 label facei = startFacei;
1089 dualFace.append(celli);
1105 if (facei == startFacei)
1120 dynDualFaces.
append(face(dualFace.shrink()));
1121 dynDualOwner.append(owner);
1122 dynDualNeighbour.append(neighbour);
1123 dynDualRegion.append(-1);
1127 const face&
f = dynDualFaces.
last();
1128 const vector areaNorm =
f.areaNormal(dualPoints);
1130 if (((
mesh.
points()[owner] - dualPoints[
f[0]]) & areaNorm) > 0)
1133 <<
" on internal edge:" << edgeI
1145 dynDualFaces.shrink();
1146 dynDualOwner.shrink();
1147 dynDualNeighbour.shrink();
1148 dynDualRegion.shrink();
1150 OFstream str(
"dualInternalFaces.obj");
1151 Pout<<
"polyDualMesh::calcDual : dumping internal faces to " 1154 forAll(dualPoints, dualPointi)
1158 forAll(dynDualFaces, dualFacei)
1160 const face&
f = dynDualFaces[dualFacei];
1165 str<<
' ' <<
f[fp]+1;
1171 const label nInternalFaces = dynDualFaces.size();
1178 const polyPatch& pp =
patches[patchi];
1199 faceList dualFaces(std::move(dynDualFaces));
1200 labelList dualOwner(std::move(dynDualOwner));
1201 labelList dualNeighbour(std::move(dynDualNeighbour));
1202 labelList dualRegion(std::move(dynDualRegion));
1208 OFstream str(
"dualFaces.obj");
1209 Pout<<
"polyDualMesh::calcDual : dumping all faces to " 1212 forAll(dualPoints, dualPointi)
1216 forAll(dualFaces, dualFacei)
1218 const face&
f = dualFaces[dualFacei];
1223 str<<
' ' <<
f[fp]+1;
1234 dualCells[celli].setSize(0);
1239 label celli = dualOwner[facei];
1243 label sz = cFaces.
size();
1244 cFaces.setSize(sz+1);
1247 forAll(dualNeighbour, facei)
1249 label celli = dualNeighbour[facei];
1255 label sz = cFaces.
size();
1256 cFaces.setSize(sz+1);
1291 forAll(dualRegion, facei)
1293 if (dualRegion[facei] >= 0)
1295 patchSizes[dualRegion[facei]]++;
1301 label facei = nInternalFaces;
1305 patchStarts[patchi] = facei;
1307 facei += patchSizes[patchi];
1311 Pout<<
"nFaces:" << dualFaces.size()
1312 <<
" patchSizes:" << patchSizes
1313 <<
" patchStarts:" << patchStarts
1322 const polyPatch& pp =
patches[patchi];
1336 addPatches(dualPatches);
1354 Foam::polyDualMesh::polyDualMesh(
const IOobject&
io)
1362 time().findInstance(meshDir(),
"cellPoint"),
1373 "boundaryFacePoint",
1374 time().findInstance(meshDir(),
"boundaryFacePoint"),
1385 Foam::polyDualMesh::polyDualMesh
1398 time().findInstance(meshDir(),
"faces"),
1410 "boundaryFacePoint",
1411 time().findInstance(meshDir(),
"faces"),
1420 calcDual(
mesh, featureEdges, featurePoints);
1425 Foam::polyDualMesh::polyDualMesh
1428 const scalar featureCos
1437 time().findInstance(meshDir(),
"faces"),
1449 "boundaryFacePoint",
1450 time().findInstance(meshDir(),
"faces"),
1462 calcDual(
mesh, featureEdges, featurePoints);
1469 const scalar featureCos,
1487 labelList allRegion(allBoundary.size());
1506 const vectorField& faceNormals = allBoundary.faceNormals();
1507 const labelList& meshPoints = allBoundary.meshPoints();
1509 bitSet isFeatureEdge(edgeFaces.size(),
false);
1513 const labelList& eFaces = edgeFaces[edgeI];
1515 if (eFaces.size() != 2)
1518 const edge&
e = allBoundary.edges()[edgeI];
1521 << meshPoints[
e[0]] <<
' ' << meshPoints[
e[1]]
1524 <<
" has more than 2 faces connected to it:" 1527 isFeatureEdge.set(edgeI);
1529 else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1531 isFeatureEdge.set(edgeI);
1535 (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1539 isFeatureEdge.set(edgeI);
1549 DynamicList<label> allFeaturePoints(pointEdges.size());
1551 forAll(pointEdges, pointi)
1553 const labelList& pEdges = pointEdges[pointi];
1555 label nFeatEdges = 0;
1559 if (isFeatureEdge.test(pEdges[i]))
1567 allFeaturePoints.append(allBoundary.meshPoints()[pointi]);
1570 featurePoints.
transfer(allFeaturePoints);
1574 OFstream str(
"featurePoints.obj");
1576 Pout<<
"polyDualMesh::calcFeatures : dumping feature points to " 1581 label pointi = featurePoints[i];
1590 allBoundary.meshEdges
1603 DynamicList<label> allFeatureEdges(isFeatureEdge.size());
1604 forAll(isFeatureEdge, edgeI)
1606 if (isFeatureEdge.test(edgeI))
1609 allFeatureEdges.append(meshEdges[edgeI]);
1612 featureEdges.
transfer(allFeatureEdges);
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
List< labelList > labelListList
A List of labelList.
const labelListList & cellEdges() const
void size(const label n)
Older name for setAddressableSize.
errorManipArg< error, int > exit(error &err, const int errNo=1)
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
label nPoints() const noexcept
Number of mesh points.
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.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
virtual const labelList & faceNeighbour() const
Return face neighbour.
const labelListList & pointEdges() const
constexpr char nl
The newline '\n' character (0x0a)
static void calcFeatures(const polyMesh &, const scalar featureCos, labelList &featureEdges, labelList &featurePoints)
Helper function to create feature edges and points based on.
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
List< face > faceList
A List of faces.
Ostream & endl(Ostream &os)
Add newline and flush stream.
label size() const noexcept
The number of elements in table.
const cellList & cells() const
label nFaces() const noexcept
Number of mesh faces.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
const labelListList & edgeCells() const
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
A list of faces which address into the list of points.
A List obtained as a section of another List.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
void setSize(const label n)
Alias for resize()
label size() const noexcept
The number of elements in the list.
virtual const labelList & faceOwner() const
Return face owner.
label nInternalFaces() const noexcept
Number of internal faces.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
virtual const faceList & faces() const
Return raw faces.
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
const T * set(const label i) const
Return const pointer to element (can be nullptr), or nullptr for out-of-range access (ie...
label nEdges() const
Number of mesh edges.
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field...
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]
T & last()
Access last element of the list, position [size()-1].
~polyDualMesh()
Destructor.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const vectorField & faceCentres() const
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
#define WarningInFunction
Report a warning using Foam::Warning.
label nCells() const noexcept
Number of mesh cells.
label start() const
Return start label of this patch in the polyMesh face list.
const polyBoundaryMesh & patches
const std::string patch
OpenFOAM patch number as a std::string.
Field< vector > vectorField
Specialisation of Field<T> for vector.
Mesh consisting of general polyhedral cells.
List< label > labelList
A List of labels.
static autoPtr< T > New(Args &&... args)
Construct autoPtr with forwarding arguments.
A patch is a list of labels that address the faces in the global face list.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Defines the attributes of an object for which implicit objectRegistry management is supported...
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< cell > cellList
A List of cells.
static constexpr const zero Zero
Global zero (0)