54 const labelList& faceOwner,
55 const labelList& faceNeighbour,
56 const cellList&
cells,
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
148void 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)
331void 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];
434void 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());
577void 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;
714void 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];
1324 dualPatches[patchi] = pp.
clone
1332 addPatches(dualPatches);
1358 time().findInstance(meshDir(),
"cellPoint"),
1369 "boundaryFacePoint",
1370 time().findInstance(meshDir(),
"boundaryFacePoint"),
1394 time().findInstance(meshDir(),
"faces"),
1406 "boundaryFacePoint",
1407 time().findInstance(meshDir(),
"faces"),
1416 calcDual(
mesh, featureEdges, featurePoints);
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());
1505 bitSet isFeatureEdge(edgeFaces.
size(),
false);
1509 const labelList& eFaces = edgeFaces[edgeI];
1511 if (eFaces.
size() != 2)
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);
1535 isFeatureEdge.
set(edgeI);
1547 forAll(pointEdges, pointi)
1549 const labelList& pEdges = pointEdges[pointi];
1551 label nFeatEdges = 0;
1555 if (isFeatureEdge.
test(pEdges[i]))
1566 featurePoints.
transfer(allFeaturePoints);
1572 Pout<<
"polyDualMesh::calcFeatures : dumping feature points to "
1577 label pointi = featurePoints[i];
1600 forAll(isFeatureEdge, edgeI)
1602 if (isFeatureEdge.
test(edgeI))
1605 allFeatureEdges.
append(meshEdges[edgeI]);
1608 featureEdges.
transfer(allFeatureEdges);
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
label size() const noexcept
The number of elements in table.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
void transfer(List< T > &list)
void setSize(const label n)
Alias for resize()
void append(const T &val)
Append an element at the end of the list.
Output to file stream, using an OSstream.
virtual const fileName & name() const
Read/write access to the name of the stream.
virtual const fileName & name() const
Get the name of the stream.
label size() const noexcept
Number of entries.
A list of faces which address into the list of points.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const labelListList & edgeFaces() const
Return edge-face addressing.
PtrList< T > clone(Args &&... args) const
Make a copy by cloning each of the list elements.
A List obtained as a section of another List.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
label size() const noexcept
The number of elements in the list.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
void set(const bitSet &bitset)
Set specified bits from another bitset.
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Creates dual of polyMesh.
~polyDualMesh()
Destructor.
static void calcFeatures(const polyMesh &, const scalar featureCos, labelList &featureEdges, labelList &featurePoints)
Helper function to create feature edges and points based on.
Mesh consisting of general polyhedral cells.
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.
A patch is a list of labels that address the faces in the global face list.
label start() const
Return start label of this patch in the polyMesh face list.
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
const labelListList & pointEdges() const
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & cellEdges() const
const vectorField & faceCentres() const
label nInternalFaces() const noexcept
Number of internal faces.
const vectorField & cellCentres() const
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
label nFaces() const noexcept
Number of mesh faces.
const labelListList & edgeCells() const
const cellList & cells() const
label nEdges() const
Number of mesh edges.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
#define WarningInFunction
Report a warning using Foam::Warning.
const std::string patch
OpenFOAM patch number as a std::string.
void inplaceRenumber(const labelUList &oldToNew, IntListType &input)
Inplace renumber the values (not the indices) of a list.
List< label > labelList
A List of labels.
List< cell > cellList
A List of cells.
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
vectorField pointField
pointField is a vectorField.
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
List< labelList > labelListList
A List of labelList.
errorManip< error > abort(error &err)
static constexpr const zero Zero
Global zero (0)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
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.