53 void Foam::triSurfaceTools::calcRefineStatus
55 const triSurface& surf,
57 List<refineType>& refine
60 if (refine[facei] == RED)
69 const labelList& myNeighbours = surf.faceFaces()[facei];
71 forAll(myNeighbours, myNeighbourI)
73 label neighbourFacei = myNeighbours[myNeighbourI];
75 if (refine[neighbourFacei] == GREEN)
78 calcRefineStatus(surf, neighbourFacei, refine);
80 else if (refine[neighbourFacei] == NONE)
82 refine[neighbourFacei] = GREEN;
90 void Foam::triSurfaceTools::greenRefine
92 const triSurface& surf,
95 const label newPointi,
96 DynamicList<labelledTri>& newFaces
99 const labelledTri& f = surf.localFaces()[facei];
100 const edge& e = surf.edges()[edgeI];
105 label fp1 = f.fcIndex(fp0);
106 label fp2 = f.fcIndex(fp1);
161 const triSurface& surf,
162 const List<refineType>& refineStatus
166 DynamicList<point> newPoints(surf.nPoints());
167 forAll(surf.localPoints(), pointi)
169 newPoints.append(surf.localPoints()[pointi]);
171 label newVertI = surf.nPoints();
174 DynamicList<labelledTri> newFaces(surf.size());
180 forAll(refineStatus, facei)
182 if (refineStatus[facei] == RED)
185 const labelList& fEdges = surf.faceEdges()[facei];
189 label edgeI = fEdges[i];
191 if (edgeMid[edgeI] == -1)
193 const edge& e = surf.edges()[edgeI];
200 surf.localPoints()[e.start()]
201 + surf.localPoints()[e.end()]
204 edgeMid[edgeI] = newVertI++;
211 const edgeList& edges = surf.edges();
219 edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
230 edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
241 edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
263 const label edgeI = fEdges[i];
267 if ((otherFacei != -1) && (refineStatus[otherFacei] == GREEN))
283 forAll(refineStatus, facei)
285 if (refineStatus[facei] == NONE)
287 newFaces.append(surf.localFaces()[facei]);
297 allPoints.transfer(newPoints);
300 return triSurface(newFaces, surf.patches(),
allPoints,
true);
306 Foam::scalar Foam::triSurfaceTools::faceCosAngle
314 const vector common(pEnd - pStart);
315 const vector base0(pLeft - pStart);
316 const vector base1(pRight - pStart);
318 vector n0(common ^ base0);
321 vector n1(base1 ^ common);
332 void Foam::triSurfaceTools::protectNeighbours
334 const triSurface& surf,
350 const labelList& myEdges = surf.pointEdges()[vertI];
353 const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
357 label facei = myFaces[myFacei];
359 if ((faceStatus[facei] ==
ANYEDGE) || (faceStatus[facei] >= 0))
361 faceStatus[facei] =
NOEDGE;
375 const triSurface& surf,
379 const edge& e = surf.edges()[edgeI];
380 label v1 = e.start();
384 const labelList& myFaces = surf.edgeFaces()[edgeI];
390 facesToBeCollapsed.insert(myFaces[myFacei]);
397 const labelList& v1Faces = surf.pointFaces()[v1];
401 label face1I = v1Faces[v1Facei];
415 facesToBeCollapsed.insert(face1I);
416 facesToBeCollapsed.insert(face2I);
421 return facesToBeCollapsed;
428 const triSurface& surf,
433 const labelList& myFaces = surf.pointFaces()[vertI];
437 label face1I = myFaces[myFacei];
439 if (faceUsed.found(face1I))
449 void Foam::triSurfaceTools::getMergedEdges
451 const triSurface& surf,
454 HashTable<
label,
label, Hash<label>>& edgeToEdge,
455 HashTable<
label,
label, Hash<label>>& edgeToFace
458 const edge& e = surf.edges()[edgeI];
459 label v1 = e.start();
462 const labelList& v1Faces = surf.pointFaces()[v1];
463 const labelList& v2Faces = surf.pointFaces()[v2];
470 if (!collapsedFaces.found(v2Faces[v2Facei]))
472 v2FacesHash.insert(v2Faces[v2Facei]);
479 label face1I = v1Faces[v1Facei];
481 if (collapsedFaces.found(face1I))
505 label commonVert = vert1I;
506 label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
510 face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
519 edgeToEdge.insert(edge1I, edge2I);
520 edgeToEdge.insert(edge2I, edge1I);
522 edgeToFace.insert(edge1I, face2I);
523 edgeToFace.insert(edge2I, face1I);
531 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
533 const triSurface& surf,
537 const HashTable<
label,
label, Hash<label>>& edgeToEdge,
538 const HashTable<
label,
label, Hash<label>>& edgeToFace,
543 const pointField& localPoints = surf.localPoints();
545 label A = surf.edges()[edgeI].start();
546 label B = surf.edges()[edgeI].end();
553 if (edgeToEdge.found(edgeI))
556 label edge2I = edgeToEdge[edgeI];
557 face2I = edgeToFace[edgeI];
566 if ((face2I != -1) && !collapsedFaces.found(face2I))
578 cosAngle = faceCosAngle
588 cosAngle = faceCosAngle
598 cosAngle = faceCosAngle
608 cosAngle = faceCosAngle
619 <<
"face " << facei <<
" does not use vertex " 627 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
629 const triSurface& surf,
633 const HashTable<
label,
label, Hash<label>>& edgeToEdge,
634 const HashTable<
label,
label, Hash<label>>& edgeToFace
637 const labelList& v1Faces = surf.pointFaces()[v1];
643 label facei = v1Faces[v1Facei];
645 if (collapsedFaces.found(facei))
650 const labelList& myEdges = surf.faceEdges()[facei];
654 label edgeI = myEdges[myEdgeI];
681 bool Foam::triSurfaceTools::collapseCreatesFold
683 const triSurface& surf,
687 const HashTable<
label,
label, Hash<label>>& edgeToEdge,
688 const HashTable<
label,
label, Hash<label>>& edgeToFace,
692 const labelList& v1Faces = surf.pointFaces()[v1];
696 label facei = v1Faces[v1Facei];
698 if (collapsedFaces.found(facei))
703 const labelList& myEdges = surf.faceEdges()[facei];
707 label edgeI = myEdges[myEdgeI];
844 const label excludeEdgeI,
845 const label excludePointi,
847 const point& triPoint,
848 const plane& cutPlane,
853 const labelledTri& f = s[triI];
854 const labelList& fEdges = s.faceEdges()[triI];
857 FixedList<scalar, 3> d;
862 d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal();
871 if (
mag(d[i]) < 1e-6)
880 if (excludePointi != -1)
892 label fp1 = f.fcIndex(fp0);
893 label fp2 = f.fcIndex(fp1);
899 cut.setPoint(points[f[fp1]]);
901 cut.setIndex(s.localFaces()[triI][fp1]);
903 else if (d[fp2] == 0.0)
906 cut.setPoint(points[f[fp2]]);
908 cut.setIndex(s.localFaces()[triI][fp2]);
912 (d[fp1] < 0 && d[fp2] < 0)
913 || (d[fp1] > 0 && d[fp2] > 0)
924 (d[fp2]*points[f[fp1]] - d[fp1]*points[f[fp2]])
928 cut.setIndex(fEdges[fp1]);
934 FixedList<surfaceLocation, 2> inters;
939 label fp1 = f.fcIndex(fp0);
946 <<
"problem : triangle has three intersections." <<
nl 947 <<
"triangle:" << f.tri(points)
950 inters[interI].setHit();
951 inters[interI].setPoint(points[f[fp0]]);
953 inters[interI].setIndex(s.localFaces()[triI][fp0]);
958 (d[fp0] < 0 && d[fp1] > 0)
959 || (d[fp0] > 0 && d[fp1] < 0)
965 <<
"problem : triangle has three intersections." <<
nl 966 <<
"triangle:" << f.tri(points)
969 inters[interI].setHit();
970 inters[interI].setPoint
972 (d[fp0]*points[f[fp1]] - d[fp1]*points[f[fp0]])
976 inters[interI].setIndex(fEdges[fp0]);
986 else if (interI == 1)
991 else if (interI == 2)
997 && inters[0].index() == excludeEdgeI
1005 && inters[1].index() == excludeEdgeI
1015 magSqr(inters[0].rawPoint() - toPoint)
1016 <
magSqr(inters[1].rawPoint() - toPoint)
1033 void Foam::triSurfaceTools::snapToEnd
1035 const triSurface& s,
1036 const surfaceLocation& end,
1037 surfaceLocation& current
1045 if (current.index() == end.index())
1063 const labelList& fEdges = s.faceEdges()[current.index()];
1065 if (
findIndex(fEdges, end.index()) != -1)
1079 if (current.index() == end.index())
1093 const edge& e = s.edges()[end.index()];
1095 if (current.index() == e[0] || current.index() == e[1])
1128 const edge& e = s.
edges()[current.index()];
1130 if (end.index() == e[0] || end.index() == e[1])
1144 if (current.index() == end.index())
1165 const triSurface& s,
1167 const surfaceLocation& start,
1168 const label excludeEdgeI,
1169 const label excludePointi,
1170 const surfaceLocation& end,
1171 const plane& cutPlane
1174 surfaceLocation nearest;
1180 label triI = eFaces[i];
1183 if (triI != start.triangle())
1190 nearest.triangle() = triI;
1197 surfaceLocation cutInfo = cutEdge
1209 if (excludeEdgeI != -1 && !cutInfo.hit())
1212 <<
"Triangle:" << triI
1213 <<
" excludeEdge:" << excludeEdgeI
1214 <<
" point:" << start.rawPoint()
1215 <<
" plane:" << cutPlane
1221 scalar distSqr =
magSqr(cutInfo.rawPoint()-end.rawPoint());
1223 if (distSqr < minDistSqr)
1225 minDistSqr = distSqr;
1227 nearest.triangle() = triI;
1235 if (nearest.triangle() == -1)
1260 const point& pt = pts[pointi];
1262 outFile<<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
1264 Pout<<
"Written " << pts.
size() <<
" vertices to file " << fName <<
endl;
1279 forAll(markedVerts, vertI)
1281 if (markedVerts[vertI])
1285 outFile<<
"v " << pt.
x() <<
' ' << pt.
y() <<
' ' << pt.
z() <<
endl;
1290 Pout<<
"Written " << nVerts <<
" vertices to file " << fName <<
endl;
1309 label face1I = myFaces[0];
1311 if (myFaces.
size() == 2)
1313 face2I = myFaces[1];
1324 forAll(startFaces, startFacei)
1326 edgeTris[nTris++] = startFaces[startFacei];
1329 forAll(endFaces, endFacei)
1331 label facei = endFaces[endFacei];
1333 if ((facei != face1I) && (facei != face2I))
1335 edgeTris[nTris++] = facei;
1360 const edge& e = edges[v1Edges[v1EdgeI]];
1368 const edge& e = edges[v2Edges[v2EdgeI]];
1372 vertexNeighbours.
insert(vertI);
1374 return vertexNeighbours.
toc();
1436 if (myFaces.
size() != 2)
1442 if (facei == myFaces[0])
1471 <<
"Edge " << surf.
edges()[edgeI] <<
" not in face " 1500 else if (vertI == f[1])
1505 else if (vertI == f[2])
1530 label edgeI = myEdges[myEdgeI];
1534 if ((e.
start() != vertI) && (e.
end() != vertI))
1541 <<
"Cannot find vertex " << vertI <<
" in edges of face " << facei
1561 label vertI = f[fp];
1563 if (vertI != e.
start() && vertI != e.
end())
1570 <<
"Cannot find vertex opposite edge " << edgeI <<
" vertices " << e
1589 label edgeI = v1Edges[v1EdgeI];
1592 if ((e.
start() == v2) || (e.
end() == v2))
1610 if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1613 <<
"Duplicate edge labels : e0:" << e0I <<
" e1:" << e1I
1622 label facei = eFaces[eFacei];
1629 || (myEdges[1] == e1I)
1630 || (myEdges[2] == e1I)
1636 || (myEdges[1] == e2I)
1637 || (myEdges[2] == e2I)
1697 return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1719 forAll(localPoints, pointi)
1721 pointMap[pointi] = pointi;
1727 forAll(collapseEdgeLabels, collapseEdgeI)
1729 const label edgeI = collapseEdgeLabels[collapseEdgeI];
1731 if ((edgeI < 0) || (edgeI >= surf.
nEdges()))
1734 <<
"Edge label outside valid range." <<
endl 1735 <<
"edge label:" << edgeI <<
endl 1736 <<
"total number of edges:" << surf.
nEdges() <<
endl 1740 const labelList& neighbours = edgeFaces[edgeI];
1742 if (neighbours.
size() == 2)
1744 const label stat0 = faceStatus[neighbours[0]];
1745 const label stat1 = faceStatus[neighbours[1]];
1750 ((stat0 ==
ANYEDGE) || (stat0 == edgeI))
1751 && ((stat1 ==
ANYEDGE) || (stat1 == edgeI))
1754 const edge& e = edges[edgeI];
1760 || (pointMap[e.
end()] != e.
end())
1764 <<
"points already mapped. Double collapse." <<
endl 1765 <<
"edgeI:" << edgeI
1766 <<
" start:" << e.
start()
1767 <<
" end:" << e.
end()
1768 <<
" pointMap[start]:" << pointMap[e.
start()]
1769 <<
" pointMap[end]:" << pointMap[e.
end()]
1774 pointMap[e.
start()] = minVert;
1775 pointMap[e.
end()] = minVert;
1778 newPoints[minVert] = edgeMids[edgeI];
1781 protectNeighbours(surf, e.
start(), faceStatus);
1782 protectNeighbours(surf, e.
end(), faceStatus);
1804 forAll(collapseFaces, collapseI)
1806 faceStatus[collapseFaces[collapseI]] =
COLLAPSED;
1821 forAll(localFaces, facei)
1825 const label a = pointMap[f[0]];
1826 const label b = pointMap[f[1]];
1827 const label c = pointMap[f[2]];
1831 (a != b) && (a != c) && (b != c)
1836 newTris[newTriI++] =
labelledTri(a, b, c, f.region());
1855 tempSurf.localFaces(),
1857 tempSurf.localPoints()
1873 forAll(refineFaces, refineFacei)
1875 calcRefineStatus(surf, refineFaces[refineFacei], refineStatus);
1879 return doRefine(surf, refineStatus);
1901 forAll(refineEdges, refineEdgeI)
1903 label edgeI = refineEdges[refineEdgeI];
1907 bool neighbourIsRefined=
false;
1911 if (refineStatus[myFaces[myFacei]] != NONE)
1913 neighbourIsRefined =
true;
1918 if (!neighbourIsRefined)
1946 refineStatus[myFaces[myFacei]] = GREEN;
1956 if (refineStatus[facei] == NONE)
1981 scalar minLength = GREAT;
1982 label minIndex = -1;
1985 const edge& e = surf.
edges()[edgeIndices[i]];
1994 if (length < minLength)
2000 return edgeIndices[minIndex];
2011 scalar maxLength = -GREAT;
2012 label maxIndex = -1;
2015 const edge& e = surf.
edges()[edgeIndices[i]];
2024 if (length > maxLength)
2030 return edgeIndices[maxIndex];
2038 const scalar mergeTol
2060 label newTriangleI = 0;
2066 label newA = pointMap[f[0]];
2067 label newB = pointMap[f[1]];
2068 label newC = pointMap[f[2]];
2070 if ((newA != newB) && (newA != newC) && (newB != newC))
2072 newTriangles[newTriangleI++] =
2076 newTriangles.
setSize(newTriangleI);
2097 const label nearestFacei,
2098 const point& nearestPt
2104 label nearType, nearLabel;
2127 return edgeNormal/(
mag(edgeNormal) + VSMALL);
2141 const point& sample,
2142 const point& nearestPoint,
2148 if (eFaces.
size() != 2)
2160 vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2162 if (((sample - nearestPoint) & n) > 0)
2178 const point& sample,
2179 const label nearestFacei
2186 label nearType, nearLabel;
2194 vector sampleNearestVec = (sample - nearestPoint);
2197 scalar
c = sampleNearestVec & surf.
faceNormals()[nearestFacei];
2269 return edgeSide(surf, sample, nearestPoint, edgeI);
2279 label nearPointi = localF[nearLabel];
2283 const point& base = localPoints[nearPointi];
2288 label minEdgeI = -1;
2292 label edgeI = pEdges[i];
2294 const edge& e = edges[edgeI];
2299 vector eVec(localPoints[otherPointi] - base);
2300 scalar magEVec =
mag(eVec);
2302 if (magEVec > VSMALL)
2307 const point perturbPoint = base + eVec;
2311 if (distSqr < minDistSqr)
2313 minDistSqr = distSqr;
2322 <<
"Problem: did not find edge closer than " << minDistSqr
2326 return edgeSide(surf, sample, nearestPoint, minEdgeI);
2347 label newPatchi = 0;
2355 label nTriTotal = 0;
2357 forAll(patch, patchFacei)
2359 const face& f = patch[patchFacei];
2367 forAll(triFaces, triFacei)
2369 const face& f = triFaces[triFacei];
2379 Pout<< patch.
name() <<
" : generated " << nTriTotal
2380 <<
" triangles from " << patch.size() <<
" faces with" 2381 <<
" new patchid " << newPatchi <<
endl;
2394 rawSurface.localFaces(),
2395 rawSurface.localPoints()
2408 surface.patches()[newPatchi].
name() = patch.
name();
2409 surface.patches()[newPatchi].geometricType() = patch.type();
2434 label newPointi = 0;
2438 newPoints[newPointi++] = points[pointi];
2440 forAll(faceCentres, facei)
2442 newPoints[newPointi++] = faceCentres[facei];
2452 label newPatchi = 0;
2459 label nTriTotal = 0;
2461 forAll(patch, patchFacei)
2464 const face& f = patch[patchFacei];
2473 triangles.append(
labelledTri(f[fp], f[fp1], fc, newPatchi));
2481 Pout<< patch.
name() <<
" : generated " << nTriTotal
2482 <<
" triangles from " << patch.size() <<
" faces with" 2483 <<
" new patchid " << newPatchi <<
endl;
2511 surface.patches()[newPatchi].
name() = patch.
name();
2512 surface.patches()[newPatchi].geometricType() = patch.type();
2529 geompackVertices[doubleI++] = pts[i][0];
2530 geompackVertices[doubleI++] = pts[i][1];
2543 geompackVertices.begin(),
2545 triangle_node.
begin(),
2546 triangle_neighbor.begin()
2552 <<
"Failed dtris2 with vertices:" << pts.
size()
2557 triangle_node.setSize(3*nTris);
2558 triangle_neighbor.setSize(3*nTris);
2567 triangle_node[3*i]-1,
2568 triangle_node[3*i+1]-1,
2569 triangle_node[3*i+2]-1,
2577 points[i][0] = pts[i][0];
2578 points[i][1] = pts[i][1];
2596 edge[0] = tri.
c()-tri.
b();
2597 edge[1] = tri.
a()-tri.
c();
2598 edge[2] = tri.
b()-tri.
a();
2600 vector triangleFaceNormal = edge[1] ^ edge[2];
2604 for (
label i=0; i<3; i++)
2606 normal[i] = triangleFaceNormal ^ edge[i];
2607 normal[i] /=
mag(normal[i]) + VSMALL;
2610 weights[0] = ((p-tri.
b()) & normal[0]) /
max(VSMALL, normal[0] & edge[1]);
2611 weights[1] = ((p-tri.
c()) & normal[1]) /
max(VSMALL, normal[1] & edge[2]);
2612 weights[2] = ((p-tri.
a()) & normal[2]) /
max(VSMALL, normal[2] & edge[0]);
2626 allVerts.setSize(samplePts.
size());
2627 allWeights.setSize(samplePts.
size());
2633 const point& samplePt = samplePts[i];
2639 scalar minDistance = GREAT;
2647 label nearType, nearLabel;
2673 else if (nearest.
distance() < minDistance)
2681 verts[0] = f[nearLabel];
2684 weights[1] = -GREAT;
2686 weights[2] = -GREAT;
2695 verts[0] = f[nearLabel];
2696 verts[1] = f[f.
fcIndex(nearLabel)];
2699 const point& p0 = points[verts[0]];
2700 const point& p1 = points[verts[1]];
2715 weights[2] = -GREAT;
2754 const point& trianglePoint
2760 label index, elemType;
2799 const plane& cutPlane
2807 snapToEnd(s, end, nearest);
2840 nearest = visitFaces
2855 nearest = visitFaces
2866 snapToEnd(s, end, nearest);
2876 const plane& cutPlane,
label end() const
Return end vertex label.
label nPoints() const
Return number of points supporting patch faces.
A HashTable with keys but without contents.
#define forAll(list, i)
Loop across all elements in list.
A triangle primitive used to calculate face normals and swept volumes.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
A class for handling file names.
A face is a list of labels corresponding to mesh vertices.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
const Point & c() const
Return third vertex.
bool hit() const
Is there a hit.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
void size(const label)
Override size to be inconsistent with allocated storage.
void setIndex(const label index)
const Point & rawPoint() const
Return point with no checking.
Ostream & endl(Ostream &os)
Add newline and flush stream.
const Field< PointType > & points() const
Return reference to global points.
const Field< PointType > & pointNormals() const
Return point normals for patch.
label otherVertex(const label a) const
Given one vertex, return the other.
Vector< scalar > vector
A scalar version of the templated Vector.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
const vectorField & faceCentres() const
bool insert(const Key &key)
Insert a new entry.
void setPoint(const Point &p)
virtual const pointField & points() const
Return raw points.
const labelListList & pointEdges() const
Return point-edge addressing.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
scalar distance() const
Return distance to hit.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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))
vectorField pointField
pointField is a vectorField.
triPointRef::proxType & elementType()
label start() const
Return start label of this patch in the polyMesh face list.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
pointHit nearestPointClassify(const point &p, const pointField &points, label &nearType, label &nearLabel) const
Return nearest point to face and classify it:
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
const geometricSurfacePatchList & patches() const
const labelListList & faceEdges() const
Return face-edge addressing.
void append(const T &)
Append an element at the end of the list.
edgeList edges() const
Return edges in face point ordering,.
const Point & a() const
Return first vertex.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
iterator begin()
Return an iterator to begin traversing the UList.
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,&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]
List< label > labelList
A List of labels.
Triangle with additional region number.
const word & name() const
Return name.
const Point & rawPoint() const
Return point with no checking.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
const polyMesh & mesh() const
Return the mesh reference.
List< Key > toc() const
Return the table of contents.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
label start() const
Return start vertex label.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
label nEdges() const
Return number of edges in patch.
A normal distribution model.
const Point & b() const
Return second vertex.
void setSize(const label)
Reset size of List.
int dtris2(int point_num, double point_xy[], int *tri_num, int tri_vert[], int tri_nabe[])
const Field< PointType > & faceNormals() const
Return face normals for patch.
vector point
Point is a vector.
triPointRef tri(const pointField &) const
Return the triangle.
label nTriangles() const
Number of triangles after splitting.
label index() const
Return index.
const dimensionedScalar c
Speed of light in a vacuum.
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
triangle< point, const point & > triPointRef
label triangles(const pointField &points, label &triI, faceList &triFaces) const
Split into triangles using existing points.
dimensioned< scalar > mag(const dimensioned< Type > &)
void setSize(const label)
Dummy setSize function.
const labelListList & edgeFaces() const
Return edge-face addressing.
bool hit() const
Is there a hit.
Mesh consisting of general polyhedral cells.
Contains information about location on a triSurface.
const labelListList & pointFaces() const
Return point-face addressing.
A patch is a list of labels that address the faces in the global face list.
Triangulated surface description with patch information.
label nInternalFaces() const
const Field< PointType > & localPoints() const
Return pointField of points in patch.
const word & name() const
Return name.
pointHit nearestPointClassify(const point &p, label &nearType, label &nearLabel) const
Find the nearest point to p on the triangle and classify it: