55 void Foam::triSurfaceTools::calcRefineStatus
57 const triSurface& surf,
59 List<refineType>& refine
62 if (refine[facei] == RED)
71 const labelList& myNeighbours = surf.faceFaces()[facei];
73 for (
const label neighbourFacei : myNeighbours)
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,
96 DynamicList<labelledTri>& newFaces
99 const labelledTri&
f = surf.localFaces()[facei];
100 const edge&
e = surf.edges()[edgeI];
161 const triSurface& surf,
162 const List<refineType>& refineStatus
166 label newVertI = surf.
nPoints();
168 DynamicList<point> newPoints(newVertI);
169 newPoints.append(surf.localPoints());
172 DynamicList<labelledTri> newFaces(surf.size());
178 forAll(refineStatus, facei)
180 if (refineStatus[facei] == RED)
183 const labelList& fEdges = surf.faceEdges()[facei];
185 for (
const label edgei : fEdges)
187 if (edgeMid[edgei] == -1)
189 const edge&
e = surf.edges()[edgei];
192 newPoints.append(
e.centre(surf.localPoints()));
193 edgeMid[edgei] = newVertI++;
200 const edgeList& edges = surf.edges();
208 edges[fEdges[0]].commonVertex(edges[fEdges[1]]),
219 edges[fEdges[1]].commonVertex(edges[fEdges[2]]),
230 edges[fEdges[2]].commonVertex(edges[fEdges[0]]),
250 for (
const label edgei : fEdges)
252 label otherFacei =
otherFace(surf, facei, edgei);
254 if ((otherFacei != -1) && (refineStatus[otherFacei] == GREEN))
270 forAll(refineStatus, facei)
272 if (refineStatus[facei] == NONE)
274 newFaces.append(surf.localFaces()[facei]);
287 return triSurface(newFaces, surf.patches(),
allPoints,
true);
293 Foam::scalar Foam::triSurfaceTools::faceCosAngle
301 const vector common(pEnd - pStart);
302 const vector base0(pLeft - pStart);
303 const vector base1(pRight - pStart);
316 void Foam::triSurfaceTools::protectNeighbours
318 const triSurface& surf,
331 for (
const label edgei : surf.pointEdges()[vertI])
333 for (
const label facei : surf.edgeFaces()[edgei])
335 if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
337 faceStatus[facei] = NOEDGE;
351 const triSurface& surf,
355 const edge&
e = surf.edges()[edgeI];
356 const label v1 =
e.start();
357 const label v2 =
e.end();
360 const labelList& myFaces = surf.edgeFaces()[edgeI];
363 facesToBeCollapsed.
insert(myFaces);
369 const labelList& v1Faces = surf.pointFaces()[v1];
371 for (
const label face1I : v1Faces)
373 label otherEdgeI = oppositeEdge(surf, face1I, v1);
376 label face2I =
otherFace(surf, face1I, otherEdgeI);
381 if (oppositeVertex(surf, face2I, otherEdgeI) == v2)
385 facesToBeCollapsed.insert(face1I);
386 facesToBeCollapsed.insert(face2I);
391 return facesToBeCollapsed;
396 Foam::label Foam::triSurfaceTools::vertexUsesFace
398 const triSurface& surf,
403 for (
const label face1I : surf.pointFaces()[vertI])
405 if (faceUsed.found(face1I))
415 void Foam::triSurfaceTools::getMergedEdges
417 const triSurface& surf,
420 Map<label>& edgeToEdge,
421 Map<label>& edgeToFace
424 const edge&
e = surf.edges()[edgeI];
425 const label v1 =
e.start();
426 const label v2 =
e.end();
428 const labelList& v1Faces = surf.pointFaces()[v1];
429 const labelList& v2Faces = surf.pointFaces()[v2];
434 for (
const label facei : v2Faces)
436 if (!collapsedFaces.found(facei))
438 v2FacesHash.
insert(facei);
443 for (
const label face1I: v1Faces)
445 if (collapsedFaces.found(face1I))
469 label commonVert = vert1I;
470 label face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
474 face2I = vertexUsesFace(surf, v2FacesHash, commonVert);
480 label edge1I =
getEdge(surf, v1, commonVert);
481 label edge2I =
getEdge(surf, v2, commonVert);
483 edgeToEdge.insert(edge1I, edge2I);
484 edgeToEdge.insert(edge2I, edge1I);
486 edgeToFace.insert(edge1I, face2I);
487 edgeToFace.insert(edge2I, face1I);
495 Foam::scalar Foam::triSurfaceTools::edgeCosAngle
497 const triSurface& surf,
501 const Map<label>& edgeToEdge,
502 const Map<label>& edgeToFace,
507 const pointField& localPoints = surf.localPoints();
509 label
A = surf.edges()[edgeI].start();
510 label
B = surf.edges()[edgeI].end();
511 label
C = oppositeVertex(surf, facei, edgeI);
517 if (edgeToEdge.found(edgeI))
520 label edge2I = edgeToEdge[edgeI];
521 face2I = edgeToFace[edgeI];
523 D = oppositeVertex(surf, face2I, edge2I);
530 if ((face2I != -1) && !collapsedFaces.found(face2I))
532 D = oppositeVertex(surf, face2I, edgeI);
542 cosAngle = faceCosAngle
552 cosAngle = faceCosAngle
562 cosAngle = faceCosAngle
572 cosAngle = faceCosAngle
583 <<
"face " << facei <<
" does not use vertex " 591 Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
593 const triSurface& surf,
597 const Map<label>& edgeToEdge,
598 const Map<label>& edgeToFace
601 const labelList& v1Faces = surf.pointFaces()[v1];
605 for (
const label facei : v1Faces)
607 if (collapsedFaces.found(facei))
612 for (
const label edgeI : surf.faceEdges()[facei])
639 bool Foam::triSurfaceTools::collapseCreatesFold
641 const triSurface& surf,
645 const Map<label>& edgeToEdge,
646 const Map<label>& edgeToFace,
650 const labelList& v1Faces = surf.pointFaces()[v1];
652 for (
const label facei : v1Faces)
654 if (collapsedFaces.found(facei))
659 const labelList& myEdges = surf.faceEdges()[facei];
661 for (
const label edgeI : myEdges)
796 const label excludeEdgeI,
797 const label excludePointi,
799 const point& triPoint,
800 const plane& cutPlane,
805 const labelledTri&
f =
s[triI];
806 const labelList& fEdges =
s.faceEdges()[triI];
809 FixedList<scalar, 3> d;
814 d[fp] = cutPlane.signedDistance(
points[
f[fp]]);
823 if (
mag(d[i]) < 1
e-6)
832 if (excludePointi != -1)
836 label fp0 =
s.localFaces()[triI].find(excludePointi);
852 cut.setIndex(
s.localFaces()[triI][fp1]);
854 else if (d[fp2] == 0.0)
858 cut.setIndex(
s.localFaces()[triI][fp2]);
862 (d[fp1] < 0 && d[fp2] < 0)
863 || (d[fp1] > 0 && d[fp2] > 0)
877 cut.setIndex(fEdges[fp1]);
883 FixedList<surfaceLocation, 2> inters;
895 <<
"problem : triangle has three intersections." <<
nl 896 <<
"triangle:" <<
f.tri(
points)
899 inters[interI].hitPoint(
points[
f[fp0]]);
901 inters[interI].setIndex(
s.localFaces()[triI][fp0]);
906 (d[fp0] < 0 && d[fp1] > 0)
907 || (d[fp0] > 0 && d[fp1] < 0)
913 <<
"problem : triangle has three intersections." <<
nl 914 <<
"triangle:" <<
f.tri(
points)
917 inters[interI].hitPoint
923 inters[interI].setIndex(fEdges[fp0]);
933 else if (interI == 1)
938 else if (interI == 2)
944 && inters[0].index() == excludeEdgeI
952 && inters[1].index() == excludeEdgeI
980 void Foam::triSurfaceTools::snapToEnd
983 const surfaceLocation&
end,
984 surfaceLocation& current
992 if (current.index() ==
end.index())
1010 const labelList& fEdges =
s.faceEdges()[current.index()];
1012 if (fEdges.found(
end.index()))
1026 if (current.index() ==
end.index())
1040 const edge&
e =
s.edges()[
end.index()];
1042 if (current.index() ==
e[0] || current.index() ==
e[1])
1075 const edge&
e =
s.edges()[current.index()];
1077 if (
end.index() ==
e[0] ||
end.index() ==
e[1])
1091 if (current.index() ==
end.index())
1112 const triSurface&
s,
1114 const surfaceLocation& start,
1115 const label excludeEdgeI,
1116 const label excludePointi,
1117 const surfaceLocation&
end,
1118 const plane& cutPlane
1121 surfaceLocation nearest;
1125 for (
const label triI : eFaces)
1128 if (triI != start.triangle())
1135 nearest.triangle() = triI;
1142 surfaceLocation cutInfo = cutEdge
1154 if (excludeEdgeI != -1 && !cutInfo.hit())
1157 <<
"Triangle:" << triI
1158 <<
" excludeEdge:" << excludeEdgeI
1159 <<
" point:" << start.point()
1160 <<
" plane:" << cutPlane
1166 scalar distSqr = cutInfo.point().distSqr(
end.point());
1168 if (distSqr < minDistSqr)
1170 minDistSqr = distSqr;
1172 nearest.triangle() = triI;
1180 if (nearest.triangle() == -1)
1205 outFile<<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
endl;
1207 Pout<<
"Written " <<
pts.
size() <<
" vertices to file " << fName <<
endl;
1214 const triSurface& surf,
1215 const fileName& fName,
1219 OFstream outFile(fName);
1222 forAll(markedVerts, vertI)
1224 if (markedVerts[vertI])
1226 const point& pt = surf.localPoints()[vertI];
1228 outFile<<
"v " << pt.x() <<
' ' << pt.y() <<
' ' << pt.z() <<
endl;
1233 Pout<<
"Written " << nVerts <<
" vertices to file " << fName <<
endl;
1252 label face1I = myFaces[0];
1254 if (myFaces.
size() == 2)
1256 face2I = myFaces[1];
1267 for (
const label facei : startFaces)
1269 edgeTris[nTris++] = facei;
1272 for (
const label facei : endFaces)
1274 if ((facei != face1I) && (facei != face2I))
1276 edgeTris[nTris++] = facei;
1285 const triSurface& surf,
1289 const edgeList& edges = surf.edges();
1290 const label v1 =
e.start();
1291 const label v2 =
e.end();
1296 const labelList& v1Edges = surf.pointEdges()[v1];
1298 for (
const label edgei : v1Edges)
1300 vertexNeighbours.insert(edges[edgei].otherVertex(v1));
1303 const labelList& v2Edges = surf.pointEdges()[v2];
1305 for (
const label edgei : v2Edges)
1307 vertexNeighbours.insert(edges[edgei].otherVertex(v2));
1309 return vertexNeighbours.toc();
1316 const triSurface& surf,
1321 const labelList& myFaces = surf.edgeFaces()[edgeI];
1323 if (myFaces.size() != 2)
1327 else if (facei == myFaces[0])
1341 const triSurface& surf,
1348 const labelList& eFaces = surf.faceEdges()[facei];
1350 label i0 = eFaces.
find(edgeI);
1355 <<
"Edge " << surf.edges()[edgeI] <<
" not in face " 1359 label i1 = eFaces.fcIndex(i0);
1360 label i2 = eFaces.fcIndex(i1);
1384 else if (vertI ==
f[1])
1389 else if (vertI ==
f[2])
1397 <<
"Vertex " << vertI <<
" not in face " <<
f <<
nl 1406 const triSurface& surf,
1411 const labelList& myEdges = surf.faceEdges()[facei];
1413 for (
const label edgei : myEdges)
1415 const edge&
e = surf.edges()[edgei];
1417 if (!
e.found(vertI))
1424 <<
"Cannot find vertex " << vertI <<
" in edges of face " << facei
1443 for (
const label pointi :
f)
1445 if (!
e.found(pointi))
1452 <<
"Cannot find vertex opposite edge " << edgeI <<
" vertices " <<
e 1469 for (
const label edgei : v1Edges)
1485 const triSurface& surf,
1491 if ((e0I == e1I) || (e0I == e2I) || (e1I == e2I))
1494 <<
"Duplicate edge labels : e0:" << e0I <<
" e1:" << e1I
1499 const labelList& eFaces = surf.edgeFaces()[e0I];
1501 for (
const label facei : eFaces)
1503 const labelList& myEdges = surf.faceEdges()[facei];
1508 || (myEdges[1] == e1I)
1509 || (myEdges[2] == e1I)
1515 || (myEdges[1] == e2I)
1516 || (myEdges[2] == e2I)
1530 const triSurface& surf,
1538 const edge&
e = surf.edges()[edgeI];
1539 edgeMids[edgeI] =
e.centre(surf.localPoints());
1543 labelList faceStatus(surf.size(), ANYEDGE);
1570 return collapseEdges(surf, collapsableEdges, edgeMids, faceStatus);
1595 for (
const label edgei : collapseEdgeLabels)
1597 if (edgei < 0 || edgei >= surf.
nEdges())
1600 <<
"Edge label outside valid range." <<
endl 1601 <<
"edge label:" << edgei <<
endl 1602 <<
"total number of edges:" << surf.
nEdges() <<
endl 1606 const labelList& neighbours = edgeFaces[edgei];
1608 if (neighbours.size() == 2)
1610 const label stat0 = faceStatus[neighbours[0]];
1611 const label stat1 = faceStatus[neighbours[1]];
1616 ((stat0 == ANYEDGE) || (stat0 == edgei))
1617 && ((stat1 == ANYEDGE) || (stat1 == edgei))
1620 const edge&
e = edges[edgei];
1625 (pointMap[
e.start()] !=
e.start())
1626 || (pointMap[
e.end()] !=
e.end())
1630 <<
"points already mapped. Double collapse." <<
endl 1631 <<
"edgei:" << edgei
1632 <<
" start:" <<
e.start()
1633 <<
" end:" <<
e.end()
1634 <<
" pointMap[start]:" << pointMap[
e.start()]
1635 <<
" pointMap[end]:" << pointMap[
e.end()]
1639 const label minVert =
min(
e.start(),
e.end());
1640 pointMap[
e.start()] = minVert;
1641 pointMap[
e.end()] = minVert;
1644 newPoints[minVert] = edgeMids[edgei];
1647 protectNeighbours(surf,
e.start(), faceStatus);
1648 protectNeighbours(surf,
e.end(), faceStatus);
1652 oppositeVertex(surf, neighbours[0], edgei),
1658 oppositeVertex(surf, neighbours[1], edgei),
1670 forAll(collapseFaces, collapseI)
1672 faceStatus[collapseFaces[collapseI]] = COLLAPSED;
1680 List<labelledTri> newTriangles(surf.
size());
1683 const List<labelledTri>& localFaces = surf.
localFaces();
1686 forAll(localFaces, facei)
1688 if (faceStatus[facei] != COLLAPSED)
1691 labelledTri
f(localFaces[facei]);
1694 f[0] = pointMap[
f[0]];
1695 f[1] = pointMap[
f[1]];
1696 f[2] = pointMap[
f[2]];
1700 newTriangles[nNewTris++] =
f;
1704 newTriangles.
resize(nNewTris);
1709 triSurface tempSurf(newTriangles, surf.
patches(), newPoints);
1714 tempSurf.localFaces(),
1716 tempSurf.localPoints()
1732 forAll(refineFaces, refineFacei)
1734 calcRefineStatus(surf, refineFaces[refineFacei], refineStatus);
1738 return doRefine(surf, refineStatus);
1760 for (
const label edgei : refineEdges)
1764 bool neighbourIsRefined=
false;
1766 for (
const label facei : myFaces)
1768 if (refineStatus[facei] != NONE)
1770 neighbourIsRefined =
true;
1775 if (!neighbourIsRefined)
1778 const edge&
e = surf.
edges()[edgei];
1783 for (
const label facei : myFaces)
1796 refineStatus[facei] = GREEN;
1806 if (refineStatus[facei] == NONE)
1815 return triSurface(newFaces, surf.
patches(), newPoints,
true);
1831 scalar minLen = GREAT;
1834 for (
const label edgei : edgeIndices)
1840 if (length < minLen)
1858 scalar maxLen = -GREAT;
1861 for (
const label edgei : edgeIndices)
1867 if (length > maxLen)
1882 const scalar mergeTol
1909 f[0] = pointMap[
f[0]];
1910 f[1] = pointMap[
f[1]];
1911 f[2] = pointMap[
f[2]];
1915 newTriangles[nNewTris++] =
f;
1918 newTriangles.
resize(nNewTris);
1939 const label nearestFacei,
1940 const point& nearestPt
1946 label nearType, nearLabel;
1948 f.nearestPointClassify(nearestPt,
points, nearType, nearLabel);
1958 label edgeI = surf.
faceEdges()[nearestFacei][nearLabel];
1965 for (
const label facei : eFaces)
1982 const triSurface& surf,
1983 const point& sample,
1984 const point& nearestPoint,
1988 const labelList& eFaces = surf.edgeFaces()[edgeI];
1990 if (eFaces.size() != 2)
1997 const vectorField& faceNormals = surf.faceNormals();
2002 vector n = 0.5*(faceNormals[eFaces[0]] + faceNormals[eFaces[1]]);
2004 if (((sample - nearestPoint) &
n) > 0)
2019 const triSurface& surf,
2020 const point& sample,
2021 const label nearestFacei
2028 label nearType, nearLabel;
2030 pointHit pHit =
f.nearestPointClassify(sample,
points, nearType, nearLabel);
2032 const point& nearestPoint = pHit.point();
2036 vector sampleNearestVec = (sample - nearestPoint);
2039 scalar
c = sampleNearestVec & surf.faceNormals()[nearestFacei];
2087 label edgeI = surf.faceEdges()[nearestFacei][nearLabel];
2111 return edgeSide(surf, sample, nearestPoint, edgeI);
2121 label nearPointi = localF[nearLabel];
2123 const edgeList& edges = surf.edges();
2124 const pointField& localPoints = surf.localPoints();
2125 const point& base = localPoints[nearPointi];
2127 const labelList& pEdges = surf.pointEdges()[nearPointi];
2130 label minEdgeI = -1;
2134 label edgeI = pEdges[i];
2136 const edge&
e = edges[edgeI];
2138 label otherPointi =
e.otherVertex(nearPointi);
2141 vector eVec(localPoints[otherPointi] - base);
2142 scalar magEVec =
mag(eVec);
2144 if (magEVec > VSMALL)
2149 const point perturbPoint = base + eVec;
2153 if (distSqr < minDistSqr)
2155 minDistSqr = distSqr;
2164 <<
"Problem: did not find edge closer than " << minDistSqr
2168 return edgeSide(surf, sample, nearestPoint, minEdgeI);
2175 const polyBoundaryMesh&
bMesh,
2189 for (
const label patchi : includePatches)
2193 for (
const face&
f :
patch)
2195 nTris +=
f.nTriangles(
points);
2201 label newPatchi = 0;
2204 for (
const label patchi : includePatches)
2209 label nTriTotal = 0;
2212 for (
const face&
f :
patch)
2218 f.triangles(
points, nTri, triFaces);
2220 for (
const face&
f : triFaces)
2223 triangles[nTris++] = labelledTri(
f[0],
f[1],
f[2], newPatchi);
2233 Pout<<
patch.name() <<
" : generated " << nTriTotal
2234 <<
" triangles from " <<
patch.size() <<
" faces with" 2235 <<
" new patchid " << newPatchi <<
endl;
2248 rawSurface.localFaces(),
2249 rawSurface.localPoints()
2257 for (
const label patchi : includePatches)
2262 surface.patches()[newPatchi].geometricType() =
patch.type();
2273 const polyBoundaryMesh&
bMesh,
2275 const boundBox& bBox,
2284 label newPatchi = 0;
2286 for (
const label patchi : includePatches)
2291 label nTriTotal = 0;
2295 const face&
f =
patch[patchFacei];
2297 if (bBox.containsAny(
points,
f))
2303 f.triangles(
points, nTri, triFaces);
2305 forAll(triFaces, triFacei)
2307 const face&
f = triFaces[triFacei];
2309 triangles.append(labelledTri(
f[0],
f[1],
f[2], newPatchi));
2318 Pout<<
patch.name() <<
" : generated " << nTriTotal
2319 <<
" triangles from " <<
patch.size() <<
" faces with" 2320 <<
" new patchid " << newPatchi <<
endl;
2333 rawSurface.localFaces(),
2334 rawSurface.localPoints()
2342 for (
const label patchi : includePatches)
2347 surface.patches()[newPatchi].geometricType() =
patch.type();
2378 forAll(faceCentres, facei)
2380 newPoints[
newPointi++] = faceCentres[facei];
2387 label newPatchi = 0;
2389 for (
const label patchi : includePatches)
2393 label nTriTotal = 0;
2398 const face&
f =
patch[patchFacei];
2407 triangles.append(labelledTri(
f[fp],
f[fp1], fc, newPatchi));
2415 Pout<<
patch.name() <<
" : generated " << nTriTotal
2416 <<
" triangles from " <<
patch.size() <<
" faces with" 2417 <<
" new patchid " << newPatchi <<
endl;
2426 triSurface rawSurface(
triangles, newPoints);
2431 rawSurface.localFaces(),
2432 rawSurface.localPoints()
2440 for (
const label patchi : includePatches)
2445 surface.patches()[newPatchi].geometricType() =
patch.type();
2458 List<double> geompackVertices(2*
pts.
size());
2462 geompackVertices[doubleI++] = pt[0];
2463 geompackVertices[doubleI++] = pt[1];
2468 List<int> triangle_node(m2*3*
pts.
size());
2469 List<int> triangle_neighbor(m2*3*
pts.
size());
2476 geompackVertices.begin(),
2478 triangle_node.begin(),
2479 triangle_neighbor.begin()
2485 <<
"Failed dtris2 with vertices:" <<
pts.
size()
2490 triangle_node.setSize(3*nTris);
2491 triangle_neighbor.setSize(3*nTris);
2494 List<labelledTri> faces(nTris);
2498 faces[i] = labelledTri
2500 triangle_node[3*i]-1,
2501 triangle_node[3*i+1]-1,
2502 triangle_node[3*i+2]-1,
2523 FixedList<scalar, 3>& weights
2528 FixedList<vector, 3> edge;
2529 edge[0] = tri.c()-tri.b();
2530 edge[1] = tri.a()-tri.c();
2531 edge[2] = tri.b()-tri.a();
2533 const vector triangleFaceNormal = edge[1] ^ edge[2];
2536 FixedList<vector, 3> normal;
2537 for (label i=0; i<3; i++)
2539 normal[i] =
normalised(triangleFaceNormal ^ edge[i]);
2542 weights[0] = ((
p-tri.b()) & normal[0]) /
max(VSMALL, normal[0] & edge[1]);
2543 weights[1] = ((
p-tri.c()) & normal[1]) /
max(VSMALL, normal[1] & edge[2]);
2544 weights[2] = ((
p-tri.a()) & normal[2]) /
max(VSMALL, normal[2] &
edge[0]);
2558 allVerts.setSize(samplePts.
size());
2559 allWeights.setSize(samplePts.
size());
2565 const point& samplePt = samplePts[i];
2570 scalar minDistance = GREAT;
2576 label nearType, nearLabel;
2578 pointHit nearest = tri.nearestPointClassify
2592 calcInterpolationWeights(tri, nearest.
point(), weights);
2602 else if (nearest.
distance() < minDistance)
2610 verts[0] =
f[nearLabel];
2613 weights[1] = -GREAT;
2615 weights[2] = -GREAT;
2624 verts[0] =
f[nearLabel];
2640 weights[2] = -GREAT;
2654 calcInterpolationWeights(tri, nearest.
point(), weights);
2681 const FaceType&
f = surf[facei];
2684 for (
const label pointi :
f)
2686 if (pointi < 0 || pointi >= surf.
points().size())
2691 <<
"triangle " << facei <<
" vertices " <<
f 2692 <<
" uses point indices outside point range 0.." 2699 if (
f[0] ==
f[1] ||
f[0] ==
f[2] ||
f[1] ==
f[2])
2704 <<
"triangle " << facei
2705 <<
" uses non-unique vertices " <<
f 2717 for (
const label nbrFacei : fFaces)
2719 if (nbrFacei <= facei)
2725 const FaceType& nbrF = surf[nbrFacei];
2730 (
f[0] == nbrF[0] ||
f[0] == nbrF[1] ||
f[0] == nbrF[2])
2731 && (
f[1] == nbrF[0] ||
f[1] == nbrF[1] ||
f[1] == nbrF[2])
2732 && (
f[2] == nbrF[0] ||
f[2] == nbrF[1] ||
f[2] == nbrF[2])
2738 <<
"triangle " << facei <<
" vertices " <<
f 2739 <<
" has the same vertices as triangle " << nbrFacei
2740 <<
" vertices " << nbrF
2754 const MeshedSurface<face>& surf,
2759 typedef face FaceType;
2760 const FaceType&
f = surf[facei];
2768 <<
" is not a triangle, it has " <<
f.
size()
2769 <<
" indices" <<
endl;
2775 for (
const label pointi :
f)
2777 if (pointi < 0 || pointi >= surf.points().size())
2782 <<
"triangle " << facei <<
" vertices " <<
f 2783 <<
" uses point indices outside point range 0.." 2790 if (
f[0] ==
f[1] ||
f[0] ==
f[2] ||
f[1] ==
f[2])
2795 <<
"triangle " << facei
2796 <<
" uses non-unique vertices " <<
f 2797 <<
" coords:" <<
f.points(surf.points()) <<
endl;
2804 const labelList& fFaces = surf.faceFaces()[facei];
2808 for (
const label nbrFacei : fFaces)
2810 if (nbrFacei <= facei)
2816 const FaceType& nbrF = surf[nbrFacei];
2821 (
f[0] == nbrF[0] ||
f[0] == nbrF[1] ||
f[0] == nbrF[2])
2822 && (
f[1] == nbrF[0] ||
f[1] == nbrF[1] ||
f[1] == nbrF[2])
2823 && (
f[2] == nbrF[0] ||
f[2] == nbrF[1] ||
f[2] == nbrF[2])
2829 <<
"triangle " << facei <<
" vertices " <<
f 2830 <<
" has the same vertices as triangle " << nbrFacei
2831 <<
" vertices " << nbrF
2832 <<
" coords:" <<
f.points(surf.points()) <<
endl;
2850 const point& trianglePoint
2856 label index, elemType;
2876 nearest.
setIndex(
s.faceEdges()[triI][index]);
2882 nearest.
setIndex(
s.localFaces()[triI][index]);
2892 const triSurface&
s,
2893 const surfaceLocation& start,
2894 const surfaceLocation&
end,
2895 const plane& cutPlane
2899 surfaceLocation nearest = start;
2903 snapToEnd(
s,
end, nearest);
2928 nearest.triangle() = start.index();
2934 const labelList& eFaces =
s.edgeFaces()[start.index()];
2936 nearest = visitFaces
2951 nearest = visitFaces
2970 const triSurface&
s,
2971 const surfaceLocation& endInfo,
2972 const plane& cutPlane,
2973 surfaceLocation& hitInfo
2988 hitInfo = trackToEdge
3003 if (hitInfo.hit() || hitInfo.triangle() == -1)
const labelListList & pointEdges() const
Return point-edge addressing.
void setMiss() noexcept
Set the hit status off.
label nPoints() const
Number of points supporting patch faces.
triPointRef::proxType & elementType() noexcept
const List< face_type > & localFaces() const
Return patch faces addressing into local point list.
void size(const label n)
Older name for setAddressableSize.
A triangle primitive used to calculate face normals and swept volumes. Uses referred points...
A class for handling file names.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const Field< point_type > & localPoints() const
Return pointField of points in patch.
void resize(const label len)
Adjust allocated size of list.
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.
List< edge > edgeList
List of edge.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Output to file stream as an OSstream, normally using std::ofstream for the actual output...
constexpr char nl
The newline '\n' character (0x0a)
Vector2D< scalar > vector2D
A 2D vector of scalars obtained from the generic Vector2D.
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool found(const T &val, label pos=0) const
Same as contains()
void setHit() noexcept
Set the hit status on.
std::vector< Triangle > triangles
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
label mergePoints(const PointList &points, labelList &pointToUnique, labelList &uniquePoints, const scalar mergeTol=SMALL, const bool verbose=false)
Calculate merge mapping, preserving the original point order. All points closer/equal mergeTol are to...
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
bool hit() const noexcept
Is there a hit.
const Field< point_type > & faceNormals() const
Return face unit normals for patch.
List< face > faceList
List of faces.
static label getEdge(List< DynamicList< label >> &pe, DynamicList< edge > &es, const label pointi, const label nextPointi)
A list of faces which address into the list of points.
const labelListList & faceFaces() const
Return face-face addressing.
const geometricSurfacePatchList & patches() const noexcept
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
void setPoint(const point_type &p)
Set the point.
labelledTri face_type
The face type (same as the underlying PrimitivePatch)
const labelListList & edgeFaces() const
Return edge-face addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
PointFrompoint toPoint(const Foam::point &p)
A triFace with additional (region) index.
Describes the interaction of a object and a (templated) point. It carries the info of a successful hi...
const Field< point_type > & points() const noexcept
Return reference to global points.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
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...
const wordList surface
Standard surface field types (scalar, vector, tensor, etc)
const labelListList & pointFaces() const
Return point-face addressing.
void setMiss(const bool eligible) noexcept
Set the hit status off and set the eligible miss status.
label nEdges() const
Number of edges in patch.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
const Field< point_type > & pointNormals() const
Return point normals for patch.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Geometric merging of points. See below.
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]
void setIndex(const label index) noexcept
Set the index.
const vectorField & faceCentres() const
int dtris2(int point_num, double point_xy[], int *tri_num, int tri_vert[], int tri_nabe[])
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar c
Speed of light in a vacuum.
scalar distance() const noexcept
Return distance to hit.
const std::string patch
OpenFOAM patch number as a std::string.
triangle< point, const point & > triPointRef
A triangle using referred points.
const point_type & point() const noexcept
Return the point, no checks.
const labelListList & faceEdges() const
Return face-edge addressing.
Field< vector > vectorField
Specialisation of Field<T> for vector.
Mesh consisting of general polyhedral cells.
Contains information about location on a triSurface.
const dimensionedScalar & D
List< label > labelList
A List of labels.
Triangulated surface description with patch information.
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))
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
PrimitivePatch< List< face >, const pointField > bMesh
Holder of faceList and points. (v.s. e.g. primitivePatch which references points) ...
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
const volScalarField & p0
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
PointHit< point > pointHit
A PointHit with a 3D point.
static constexpr const zero Zero
Global zero (0)