39 { triangulationMode::tmFan,
"fan" },
40 { triangulationMode::tmMesh,
"mesh" },
43 Foam::scalar Foam::faceAreaIntersect::tol = 1
e-6;
47 void Foam::faceAreaIntersect::triSliceWithPlane
51 FixedList<triPoints, 10>& tris,
57 FixedList<scalar, 3> d;
67 d[i] = pln.signedDistance(tri[i]);
69 if (
mag(d[i]) < tol*len)
95 || ((nPos == 2) && (nCoPlanar == 1))
96 || ((nPos == 1) && (nCoPlanar == 2))
111 else if ((nPos == 2) && (nCoPlanar == 0))
130 label i1 = d.fcIndex(i0);
131 label i2 = d.fcIndex(i1);
134 point p01 = planeIntersection(d, tri, i0, i1);
135 point p02 = planeIntersection(d, tri, i0, i2);
139 setTriPoints(tri[i1], tri[i2], p02, nTris, tris);
140 setTriPoints(tri[i1], p02, p01, nTris, tris);
163 label i1 = d.fcIndex(i0);
164 label i2 = d.fcIndex(i1);
167 point p01 = planeIntersection(d, tri, i1, i0);
168 point p02 = planeIntersection(d, tri, i2, i0);
171 setTriPoints(tri[i0], p01, p02, nTris, tris);
194 point p01 = planeIntersection(d, tri, i1, i0);
197 if (d.fcIndex(i0) == i1)
199 setTriPoints(tri[i0], p01, tri[i2], nTris, tris);
203 setTriPoints(tri[i0], tri[i2], p01, nTris, tris);
222 void Foam::faceAreaIntersect::triangleIntersect
224 const triPoints& src,
234 FixedList<triPoints, 10> workTris1;
235 label nWorkTris1 = 0;
237 FixedList<triPoints, 10> workTris2;
238 label nWorkTris2 = 0;
243 const scalar srcArea(src.mag());
244 if (srcArea < ROOTVSMALL)
250 const scalar t =
sqrt(srcArea);
257 scalar
s =
mag(tgt1 - tgt0);
266 const vector n0((tgt0 - tgt1)^(-
s*
n));
267 const scalar magSqrN0(
magSqr(n0));
268 if (magSqrN0 < ROOTVSMALL)
277 plane pl0(tgt0, n0/
Foam::sqrt(magSqrN0),
false);
278 triSliceWithPlane(src, pl0, workTris1, nWorkTris1, t);
292 scalar
s =
mag(tgt2 - tgt1);
297 const vector n1((tgt1 - tgt2)^(-
s*
n));
298 const scalar magSqrN1(
magSqr(n1));
300 if (magSqrN1 < ROOTVSMALL)
309 plane pl1(tgt1, n1/
Foam::sqrt(magSqrN1),
false);
313 for (label i = 0; i < nWorkTris1; ++i)
315 triSliceWithPlane(workTris1[i], pl1, workTris2, nWorkTris2, t);
330 scalar
s =
mag(tgt2 - tgt0);
335 const vector n2((tgt2 - tgt0)^(-
s*
n));
336 const scalar magSqrN2(
magSqr(n2));
338 if (magSqrN2 < ROOTVSMALL)
347 plane pl2(tgt2, n2/
Foam::sqrt(magSqrN2),
false);
351 for (label i = 0; i < nWorkTris2; ++i)
353 triSliceWithPlane(workTris2[i], pl2, workTris1, nWorkTris1, t);
363 for (label i = 0; i < nWorkTris1; ++i)
366 const scalar currArea = workTris1[i].mag();
370 centroid += currArea*workTris1[i].centre();
372 if (cacheTriangulation_)
374 triangles_.
append(workTris1[i]);
389 const DynamicList<face>& trisA,
390 const DynamicList<face>& trisB,
392 const bool cacheTriangulation
400 cacheTriangulation_(cacheTriangulation),
401 triangles_(cacheTriangulation ? 10 : 0)
411 const triangulationMode& triMode,
415 faceTris.
resize(
f.nTriangles());
419 case triangulationMode::tmFan:
421 for (label i = 0; i <
f.nTriangles(); ++i)
423 faceTris[i] =
face(3);
424 faceTris[i][0] =
f[0];
425 faceTris[i][1] =
f[i + 1];
426 faceTris[i][2] =
f[i + 2];
431 case triangulationMode::tmMesh:
433 const label nFaceTris =
f.nTriangles();
435 label nFaceTris1 = 0;
436 const label nFaceTris2 =
f.triangles(
points, nFaceTris1, faceTris);
438 if (nFaceTris != nFaceTris1 || nFaceTris != nFaceTris2)
441 <<
"The numbers of reported triangles in the face do not " 442 <<
"match that generated by the triangulation" 459 if (cacheTriangulation_)
465 centroid = vector::zero;
468 for (
const face& triA : trisA_)
470 triPoints tpA = getTriPoints(pointsA_, triA,
false);
472 for (
const face& triB : trisB_)
516 const scalar threshold
523 for (
const face& triA : trisA_)
525 const triPoints tpA = getTriPoints(pointsA_, triA,
false);
527 for (
const face& triB : trisB_)
556 if (
area > threshold)
errorManipArg< error, int > exit(error &err, const int errNo=1)
A face is a list of labels corresponding to mesh vertices.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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.
bool overlaps(const face &faceA, const face &faceB, const vector &n, const scalar threshold) const
Return area of intersection of faceA with faceB.
void append(const T &val)
Append an element at the end of the list.
dimensionedScalar sqrt(const dimensionedScalar &ds)
void calc(const face &faceA, const face &faceB, const vector &n, scalar &area, vector ¢roid) const
Return area of intersection of faceA with faceB and effective face centre.
#define forAll(list, i)
Loop across all elements in list.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
const wordList area
Standard area field types (scalar, vector, tensor, etc)
vector point
Point is a vector.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
static const Enum< triangulationMode > triangulationModeNames_
static void triangulate(const face &f, const pointField &points, const triangulationMode &triMode, faceList &faceTris)
Triangulate a face using the given triangulation mode.
faceAreaIntersect(const pointField &pointsA, const pointField &pointsB, const DynamicList< face > &trisA, const DynamicList< face > &trisB, const bool reverseB=false, const bool cacheTriangulation=false)
Construct from components.
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))
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)