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();
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;
1240 label celli = dualOwner[facei];
1244 dualCells[celli].push_back(facei);
1247 forAll(dualNeighbour, facei)
1249 label celli = dualNeighbour[facei];
1253 dualCells[celli].push_back(facei);
1287 forAll(dualRegion, facei)
1289 if (dualRegion[facei] >= 0)
1291 patchSizes[dualRegion[facei]]++;
1297 label facei = nInternalFaces;
1301 patchStarts[patchi] = facei;
1303 facei += patchSizes[patchi];
1307 Pout<<
"nFaces:" << dualFaces.size()
1308 <<
" patchSizes:" << patchSizes
1309 <<
" patchStarts:" << patchStarts
1332 addPatches(dualPatches);
1350 Foam::polyDualMesh::polyDualMesh(
const IOobject&
io)
1358 time().findInstance(meshDir(),
"cellPoint"),
1369 "boundaryFacePoint",
1370 time().findInstance(meshDir(),
"boundaryFacePoint"),
1381 Foam::polyDualMesh::polyDualMesh
1394 time().findInstance(meshDir(),
"faces"),
1406 "boundaryFacePoint",
1407 time().findInstance(meshDir(),
"faces"),
1416 calcDual(
mesh, featureEdges, featurePoints);
1421 Foam::polyDualMesh::polyDualMesh
1424 const scalar featureCos
1433 time().findInstance(meshDir(),
"faces"),
1445 "boundaryFacePoint",
1446 time().findInstance(meshDir(),
"faces"),
1458 calcDual(
mesh, featureEdges, featurePoints);
1465 const scalar featureCos,
1483 labelList allRegion(allBoundary.size());
1502 const vectorField& faceNormals = allBoundary.faceNormals();
1503 const labelList& meshPoints = allBoundary.meshPoints();
1505 bitSet isFeatureEdge(edgeFaces.size(),
false);
1509 const labelList& eFaces = edgeFaces[edgeI];
1511 if (eFaces.size() != 2)
1514 const edge&
e = allBoundary.edges()[edgeI];
1517 << meshPoints[
e[0]] <<
' ' << meshPoints[
e[1]]
1520 <<
" has more than 2 faces connected to it:" 1523 isFeatureEdge.set(edgeI);
1525 else if (allRegion[eFaces[0]] != allRegion[eFaces[1]])
1527 isFeatureEdge.set(edgeI);
1531 (faceNormals[eFaces[0]] & faceNormals[eFaces[1]])
1535 isFeatureEdge.set(edgeI);
1545 DynamicList<label> allFeaturePoints(pointEdges.size());
1547 forAll(pointEdges, pointi)
1549 const labelList& pEdges = pointEdges[pointi];
1551 label nFeatEdges = 0;
1555 if (isFeatureEdge.test(pEdges[i]))
1563 allFeaturePoints.append(allBoundary.meshPoints()[pointi]);
1566 featurePoints.
transfer(allFeaturePoints);
1570 OFstream str(
"featurePoints.obj");
1572 Pout<<
"polyDualMesh::calcFeatures : dumping feature points to " 1577 label pointi = featurePoints[i];
1586 allBoundary.meshEdges
1599 DynamicList<label> allFeatureEdges(isFeatureEdge.size());
1600 forAll(isFeatureEdge, edgeI)
1602 if (isFeatureEdge.test(edgeI))
1605 allFeatureEdges.append(meshEdges[edgeI]);
1608 featureEdges.
transfer(allFeatureEdges);
const labelListList & cellEdges() const
void size(const label n)
Older name for setAddressableSize.
List< cell > cellList
List of cell.
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.
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.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const cellList & cells() const
label nFaces() const noexcept
Number of mesh faces.
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.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
List< face > faceList
List of faces.
const labelListList & edgeCells() const
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
label size() const noexcept
The number of elements in table.
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.
void setSize(const label n)
Alias for resize()
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
label size() const noexcept
The number of entries 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.
label find(const T &val) const
Find index of the first occurrence of the value.
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 within a list.
#define WarningInFunction
Report a warning using Foam::Warning.
label nCells() const noexcept
Number of mesh cells.
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.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, IOobject::NO_REGISTER)
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.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
static constexpr const zero Zero
Global zero (0)