94#include "CGAL/version.h"
95#if defined(CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1050211000)
96#define BOOST_BIND_GLOBAL_PLACEHOLDERS
99#include <CGAL/AABB_tree.h>
100#include <CGAL/AABB_traits.h>
101#include <CGAL/AABB_face_graph_triangle_primitive.h>
104typedef CGAL::AABB_face_graph_triangle_primitive
108typedef CGAL::AABB_traits<K, Primitive> Traits;
109typedef CGAL::AABB_tree<Traits> Tree;
111typedef boost::optional
113 Tree::Intersection_and_primitive_id<Segment>::Type
114> Segment_intersection;
123bool intersectSurfaces
129 bool hasMoved =
false;
131 for (label iter = 0; iter < 10; iter++)
133 Info<<
"Determining intersections of surface edges with itself" <<
endl;
144 max(1
e-3*edgeIntersections::minEdgeLength(surf), SMALL)
182 Info<<
"Surface has been moved. Writing to " << newFile <<
endl;
192bool intersectSurfaces
200 bool hasMoved1 =
false;
201 bool hasMoved2 =
false;
203 for (label iter = 0; iter < 10; iter++)
205 Info<<
"Determining intersections of surf1 edges with surf2"
217 max(1
e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
251 Info<<
"Determining intersections of surf2 edges with surf1"
260 max(1
e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
294 if (nIters1 == 0 && nIters2 == 0)
296 Info<<
"** Resolved all intersections to be proper edge-face pierce"
305 Info<<
"Surface 1 has been moved. Writing to " << newFile
307 surf1.
write(newFile);
313 Info<<
"Surface 2 has been moved. Writing to " << newFile
315 surf2.
write(newFile);
318 return hasMoved1 || hasMoved2;
322label calcNormalDirection
325 const vector& otherNormal,
335 label nDir = ((
cross & fC0tofE0) > 0.0 ? 1 : -1);
337 nDir *= ((otherNormal & fC0tofE0) > 0.0 ? -1 : 1);
366 Info<<
"Determining intersections of surf1 edges with surf2 faces"
373 max(1
e-3*edgeIntersections::minEdgeLength(surf1), SMALL)
378 Info<<
"Determining intersections of surf2 edges with surf1 faces"
385 max(1
e-3*edgeIntersections::minEdgeLength(surf2), SMALL)
400 const label startEdgeI,
401 const label startFaceI,
405 const labelList& eFaces =
s.edgeFaces()[startEdgeI];
407 if (eFaces.
size() == 2)
410 if (eFaces[0] == startFaceI)
412 nextFaceI = eFaces[1];
414 else if (eFaces[1] == startFaceI)
416 nextFaceI = eFaces[0];
421 <<
"problem" <<
exit(FatalError);
428 label index =
s.pointFaces()[pointI].find(nextFaceI);
430 if (pFacesZone[index] == -1)
433 pFacesZone[index] = zoneI;
436 const labelList& fEdges =
s.faceEdges()[nextFaceI];
438 label nextEdgeI = -1;
442 label edgeI = fEdges[i];
443 const edge&
e =
s.edges()[edgeI];
445 if (edgeI != startEdgeI && (
e[0] == pointI ||
e[1] == pointI))
456 <<
"Problem: cannot find edge out of " << fEdges
457 <<
"on face " << nextFaceI <<
" that uses point " << pointI
458 <<
" and is not edge " << startEdgeI <<
abort(FatalError);
487 label nNonManifold = 0;
498 for (; index < pFacesZone.
size(); index++)
500 if (pFacesZone[index] == -1)
503 pFacesZone[index] = zoneI;
505 label faceI =
pFaces[index];
511 label edgeI = fEdges[fEdgeI];
512 const edge&
e = edges[edgeI];
514 if (
e[0] == pointI ||
e[1] == pointI)
534 for (label zoneI = 1; zoneI <
nZones; zoneI++)
536 label newPointI = newPoints.size();
537 newPointMap.append(
s.meshPoints()[pointI]);
538 newPoints.append(
s.points()[
s.meshPoints()[pointI]]);
542 if (pFacesZone[index] == zoneI)
544 label faceI =
pFaces[index];
548 if (localF[fp] == pointI)
550 newFaces[faceI][fp] = newPointI;
561 Info<<
"Duplicating " << nNonManifold <<
" points out of " <<
s.nPoints()
563 if (nNonManifold > 0)
565 triSurface dupSurf(newFaces,
s.patches(), newPoints,
true);
570 const labelList& dupMp = dupSurf.meshPoints();
575 label dupMeshPointI = dupMp[pointI];
576 label meshPointI = newPointMap[dupMeshPointI];
577 dupPointMap[pointI] = mpm[meshPointI];
581 forAll(dupPointMap, pointI)
583 const point& dupPt = dupSurf.points()[dupMp[pointI]];
584 label sLocalPointI = dupPointMap[pointI];
585 label sMeshPointI =
s.meshPoints()[sLocalPointI];
586 const point& sPt =
s.points()[sMeshPointI];
588 if (
mag(dupPt-sPt) > SMALL)
627 std::vector<Segment_intersection> segments;
630 const edge&
e = edges[edgeI];
635 K::Segment_3 segment_query
642 tree.all_intersections(segment_query, std::back_inserter(segments));
645 for (
const Segment_intersection& intersect : segments)
650 const Point* ptPtr = boost::get<Point>(&(intersect->first))
655 CGAL::to_double(ptPtr->x()),
656 CGAL::to_double(ptPtr->y()),
657 CGAL::to_double(ptPtr->z())
660 #if defined (CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1041400000)
661 Polyhedron::Face_handle
f = (intersect->second);
664 Polyhedron::Face_handle
f = (intersect->second).first;
667 intersections[edgeI].append
677 classifications[edgeI].append(-1);
682 const Segment* sPtr = boost::get<Segment>(&(intersect->first))
685 #if defined (CGAL_VERSION_NR) && (CGAL_VERSION_NR < 1041400000)
686 Polyhedron::Face_handle
f = (intersect->second);
689 Polyhedron::Face_handle
f = (intersect->second).first;
701 CGAL::to_double(sPtr->source().x()),
702 CGAL::to_double(sPtr->source().y()),
703 CGAL::to_double(sPtr->source().z())
708 CGAL::to_double(sPtr->target().x()),
709 CGAL::to_double(sPtr->target().y()),
710 CGAL::to_double(sPtr->target().z())
714 intersections[edgeI].append
724 classifications[edgeI].append(2);
753 const Tree tree(
p.facets_begin(),
p.facets_end(),
p);
761 for (label iter = 0; iter < 5; iter++)
764 Info<<
"Determining intersections of " << surf1.
nEdges()
765 <<
" edges with surface of " << label(tree.size()) <<
" triangles"
767 cuts = edgeIntersectionsCGAL
774 Info<<
"Determined intersections:" <<
nl
775 <<
" edges : " << surf1.
nEdges() <<
nl
776 <<
" non-degenerate cuts : " << cuts.first() <<
nl
777 <<
" degenerate cuts : " << cuts.second() <<
nl
780 if (cuts.second() == 0)
785 Info<<
"Shuffling conflicting points" <<
endl;
790 const point p05(0.5, 0.5, 0.5);
799 const edge&
e = edges[edgeI];
803 surf1Points[
mp[
e[eI]]] += surf1PointTol[
e[eI]]*d;
826 <<
"problem" <<
exit(FatalError);
837 forAll(subEdges, subEdgeI)
839 const edge& subE = subEdges[subEdgeI];
842 label start = pointMap[subE[0]];
843 label
end = pointMap[subE[1]];
844 const labelList& pEdges = pointEdges[start];
847 label edgeI = pEdges[pEdgeI];
848 const edge&
e = edges[edgeI];
850 if (
e.otherVertex(start) == end)
852 if (edgeMap[subEdgeI] == -1)
854 edgeMap[subEdgeI] = edgeI;
856 else if (edgeMap[subEdgeI] != edgeI)
859 << subE <<
" points:"
861 <<
" matches to " << edgeI
863 <<
" but also " << edgeMap[subEdgeI]
865 << edges[edgeMap[subEdgeI]].line(surf.
localPoints())
872 if (edgeMap[subEdgeI] == -1)
875 << subE <<
" at:" << subSurf.
localPoints()[subE[0]]
898 Info<<
"Intersect surface 1 edges with surface 2:" <<
nl;
899 Info<<
" constructing CGAL surface ..." <<
endl;
903 Info<<
" constructing CGAL tree ..." <<
endl;
904 const Tree tree(
p.facets_begin(),
p.facets_end(),
p);
906 edgeIntersectionsCGAL
916 Info<<
"Intersect surface 2 edges with surface 1:" <<
nl;
917 Info<<
" constructing CGAL surface ..." <<
endl;
921 Info<<
" constructing CGAL tree ..." <<
endl;
922 const Tree tree(
p.facets_begin(),
p.facets_end(),
p);
924 edgeIntersectionsCGAL
938 max(1
e-8*edgeIntersections::minEdgeLength(surf1), SMALL)
942 max(1
e-8*edgeIntersections::minEdgeLength(surf2), SMALL)
951 for (label iter = 0; iter < 10; iter++)
954 cuts1 = edgeIntersectionsAndShuffleCGAL
964 cuts2 = edgeIntersectionsAndShuffleCGAL
982void calcEdgeCutsBitsCGAL
995 PatchTools::checkOrientation(surf1,
false, &orientationEdge);
996 nZones1 = PatchTools::markZones(surf1, orientationEdge, zone1);
998 Info<<
"Surface triangles " << surf1.
size()
999 <<
" number of zones (orientation or non-manifold):"
1007 PatchTools::checkOrientation(surf2,
false, &orientationEdge);
1008 nZones2 = PatchTools::markZones(surf2, orientationEdge, zone2);
1010 Info<<
"Surface triangles " << surf2.
size()
1011 <<
" number of zones (orientation or non-manifold):"
1016 if (nZones1 == 1 && nZones2 == 1)
1041 for (label zone1I = 0; zone1I < nZones1; zone1I++)
1048 if (zone1[faceI] == zone1I)
1050 includeMap1[faceI] =
true;
1069 dupNonManifoldPoints(subSurf1, pointMap1);
1074 subSurf1.meshPoints(),
1089 for (label zone2I = 0; zone2I < nZones2; zone2I++)
1096 if (zone2[faceI] == zone2I)
1098 includeMap2[faceI] =
true;
1118 subSurf2.meshPoints(),
1123 if (!subBb2.overlaps(subBb1))
1130 dupNonManifoldPoints(subSurf2, pointMap2);
1160 label subMeshPointI = subSurf2.meshPoints()[i];
1161 label meshPointI = surf2.
meshPoints()[pointMap2[i]];
1162 points2[meshPointI] = subSurf2.points()[subMeshPointI];
1168 edgeCuts1.
merge(subEdgeCuts1, edgeMap1, faceMap2);
1169 edgeCuts2.
merge(subEdgeCuts2, edgeMap2, faceMap1);
1178 label subMeshPointI = subSurf1.meshPoints()[i];
1179 label meshPointI = surf1.
meshPoints()[pointMap1[i]];
1180 points1[meshPointI] = subSurf1.points()[subMeshPointI];
1245 const bool surf1Baffle,
1246 const bool surf2Baffle,
1247 const bool invertedSpace,
1263 label nFeatEds = inter.cutEdges().size();
1281 const label cutEdgeI = iter.val();
1283 const edge& fE = inter.cutEdges()[cutEdgeI];
1291 edgeDirections[cutEdgeI] = fE.
vec(inter.cutPoints());
1293 normals.append(norm1);
1294 eNormals.
append(normals.size() - 1);
1298 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1304 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1311 edgeDirections[cutEdgeI],
1313 inter.cutPoints()[fE.
start()]
1318 normals.append(norm2);
1319 eNormals.
append(normals.size() - 1);
1323 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1329 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1337 edgeDirections[cutEdgeI],
1339 inter.cutPoints()[fE.
start()]
1347 normals.append(norm2);
1351 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1357 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1365 edgeDirections[cutEdgeI],
1367 inter.cutPoints()[fE.
start()]
1372 eNormals.
append(normals.size() - 1);
1377 normals.append(norm1);
1381 normalVolumeTypes.append(extendedFeatureEdgeMesh::BOTH);
1387 normalVolumeTypes.append(extendedFeatureEdgeMesh::INSIDE);
1395 edgeDirections[cutEdgeI],
1397 inter.cutPoints()[fE.
start()]
1402 eNormals.
append(normals.size() - 1);
1407 label internalStart = -1;
1408 label nIntOrExt = 0;
1411 label nMultiple = 0;
1415 label nEdNorms = edgeNormals[eI].size();
1421 else if (nEdNorms == 2)
1423 const vector& n0(normals[edgeNormals[eI][0]]);
1424 const vector& n1(normals[edgeNormals[eI][1]]);
1426 if ((n0 & n1) > extendedFeatureEdgeMesh::cosNormalAngleTol_)
1435 else if (nEdNorms > 2)
1441 if (action == booleanSurface::UNION)
1451 internalStart = nIntOrExt;
1454 else if (action == booleanSurface::INTERSECTION)
1459 internalStart = nIntOrExt;
1467 else if (action == booleanSurface::DIFFERENCE)
1470 internalStart = nIntOrExt;
1475 <<
"Unsupported booleanSurface:booleanOpType and space "
1476 << action <<
" " << invertedSpace
1477 <<
abort(FatalError);
1490 forAll(edgeNormalsTmp, i)
1492 edgeNormalsTmp[i] = edgeNormals[i];
1495 forAll(normalDirectionsTmp, i)
1497 normalDirectionsTmp[i] = normalDirections[i];
1515 nIntOrExt + nFlat + nOpen,
1518 normalVolumeTypesTmp,
1520 normalDirectionsTmp,
1529int main(
int argc,
char *argv[])
1533 "Generates a extendedFeatureEdgeMesh for the interface created by"
1534 " a boolean operation on two surfaces."
1536 " [Compiled with CGAL]"
1538 " [Compiled without CGAL]"
1542 argList::noParallel();
1543 argList::addArgument
1546 "One of (intersection | union | difference)"
1548 argList::addArgument(
"surface1",
"The input surface file 1");
1549 argList::addArgument(
"surface2",
"The input surface file 2");
1555 "Geometry scaling factor (both surfaces)"
1557 argList::addBoolOption
1560 "Mark surface 1 as a baffle"
1563 argList::addBoolOption
1566 "Mark surface 2 as a baffle"
1569 argList::addBoolOption
1572 "Perturb surface points to escape degenerate intersections"
1575 argList::addBoolOption
1579 "Do not use CGAL algorithms"
1581 "Ignored, compiled without CGAL"
1585 argList::addBoolOption
1588 "Do the surfaces have inverted space orientation, "
1589 "i.e. a point at infinity is considered inside. "
1590 "This is only sensible for union and intersection."
1596 "((surface1 volumeType) .. (surfaceN volumeType))",
1597 "Trim resulting intersection with additional surfaces;"
1598 " volumeType is 'inside' (keep (parts of) edges that are inside)"
1599 ", 'outside' (keep (parts of) edges that are outside) or"
1600 " 'mixed' (keep all)"
1610 { booleanSurface::INTERSECTION,
"intersection" },
1611 { booleanSurface::UNION,
"union" },
1612 { booleanSurface::DIFFERENCE,
"difference" }
1615 if (!validActions.found(action))
1618 <<
"Unsupported action " << action <<
endl
1619 <<
"Supported actions:" << validActions <<
nl
1620 <<
abort(FatalError);
1627 Info<<
"Trimming edges with " << surfaceAndSide <<
endl;
1635 Info<<
"Reading surface " << surf1Name <<
endl;
1642 triSurfaceMesh::meshSubDir,
1646 if (scaleFactor > 0)
1648 Info<<
"Scaling : " << scaleFactor <<
nl;
1652 Info<< surf1Name <<
" statistics:" <<
endl;
1657 Info<<
"Reading surface " << surf2Name <<
endl;
1664 triSurfaceMesh::meshSubDir,
1668 if (scaleFactor > 0)
1670 Info<<
"Scaling : " << scaleFactor <<
nl;
1674 Info<< surf2Name <<
" statistics:" <<
endl;
1678 const bool surf1Baffle =
args.
found(
"surf1Baffle");
1679 const bool surf2Baffle =
args.
found(
"surf2Baffle");
1684 const bool invertedSpace =
args.
found(
"invertedSpace");
1686 if (invertedSpace && validActions[action] == booleanSurface::DIFFERENCE)
1689 <<
"Inverted space only makes sense for union or intersection."
1690 <<
exit(FatalError);
1698 calcEdgeCutsBitsCGAL
1724 +
fileName(surf2Name).nameLessExt()
1735 sFeatFileName +
".extendedFeatureEdgeMesh",
1737 "extendedFeatureEdgeMesh",
1742 booleanSurface::booleanOpTypeNames[action],
1755 forAll(surfaceAndSide, trimI)
1757 const word& trimName = surfaceAndSide[trimI].
first();
1760 volumeType::names[surfaceAndSide[trimI].second()]
1763 Info<<
"Reading trim surface " << trimName <<
endl;
1770 triSurfaceMesh::meshSubDir,
1772 IOobject::MUST_READ,
1777 Info<< trimName <<
" statistics:" <<
endl;
1778 trimSurf.writeStats(Info);
1796 Info<<
nl <<
"Writing extendedFeatureEdgeMesh to "
1808 sFeatFileName +
".eMesh",
1810 triSurfaceMesh::meshSubDir,
1820 Info<<
nl <<
"Writing featureEdgeMesh to "
1821 << bfeMesh.objectRelPath() <<
endl;
1823 bfeMesh.regIOobject::write();
CGAL data structures used for triSurface handling.
CGAL::Segment_3< K > Segment
CGAL::Polyhedron_3< K, My_items > Polyhedron
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.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
T & first() noexcept
The first element of the list, position [0].
Defines the attributes of an object for which implicit objectRegistry management is supported,...
fileName objectRelPath() const
The object path relative to the root.
fileName path() const
The complete path.
A HashTable to objects of type <T> with a label key.
const T & second() const noexcept
Return second element, which is also the last element.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
label nEdges() const
Number of edges in patch.
label nPoints() const
Number of points supporting patch faces.
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & localPoints() const
Return pointField of points in patch.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
A List with indirect addressing. Like IndirectList but does not store addressing.
T & first()
Return the first element of the list.
void size(const label n)
Older name for setAddressableSize.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
bool found(const word &optName) const
Return true if the named option is found.
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
booleanOpType
Enumeration listing the possible volume operator types.
A bounding box defined in terms of min/max extrema points.
Holder of intersections of edges of a surface with another surface. Optionally shuffles around points...
const labelListList & classification() const
For every intersection the classification status.
void merge(const edgeIntersections &, const labelList &edgeMap, const labelList &faceMap, const bool merge=true)
Merge (or override) edge intersection for a subset.
label removeDegenerates(const label nIters, const triSurface &surf1, const triSurfaceSearch &query2, const scalarField &surf1PointTol, pointField &points1)
Resolve ties. Shuffles points so all edge - face intersections.
const pointField & points() const noexcept
Return points.
const edgeList & edges() const noexcept
Return edges.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
vector vec(const UList< point > &pts) const
Return the vector (end - start)
linePointRef line(const UList< point > &pts) const
Return edge line.
label start() const
Return start (first) vertex label.
void writeObj(const fileName &prefix) const
Write all components of the extendedEdgeMesh as obj files.
virtual void writeStats(Ostream &os) const
Dump some information.
virtual bool write(const bool valid=true) const
Give precedence to the regIOobject write.
A class for handling file names.
A triFace with additional (region) index.
Basic surface-surface intersection description. Constructed from two surfaces it creates a descriptio...
IOoject and searching on triSurface.
Helper class to search on triSurface.
Triangulated surface description with patch information.
triSurface subsetMesh(const UList< bool > &include, labelList &pointMap, labelList &faceMap) const
Return a new surface subsetted on the selected faces.
void write(Ostream &os) const
Write to Ostream in simple OpenFOAM format.
void writeStats(Ostream &os) const
Write some statistics.
virtual void movePoints(const pointField &pts)
Move points.
virtual void scalePoints(const scalar scaleFactor)
Scale points. A non-positive factor is ignored.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
An enumeration wrapper for classification of a location as being inside/outside of a volume.
A class for handling words, derived from Foam::string.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
const dimensionedScalar mp
Proton mass.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
errorManip< error > abort(error &err)
void cross(FieldField< Field1, typename crossProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
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]
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.