33 template<
class SourcePatch,
class TargetPatch>
35 minCutRatio_ = 10*small;
37 template<
class SourcePatch,
class TargetPatch>
44 template<
class SourcePatch,
class TargetPatch>
45 template<
unsigned Size>
48 const cutTriList<Size>& tris,
52 OFstream obj(typeName +
"_" + name +
".vtk");
54 obj <<
"# vtk DataFile Version 2.0" <<
endl 55 <<
"triSurface" <<
endl 57 <<
"DATASET POLYDATA" <<
endl 58 <<
"POINTS " << 3*tris.size() <<
" float" <<
endl;
60 for (
label i = 0; i < tris.size(); ++ i)
62 for (
label j = 0; j < 3; ++ j)
64 const vector& x = tris[i][j];
65 obj << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
69 obj <<
"POLYGONS " << tris.size() <<
' ' << 4*tris.size() <<
endl;
71 for (
label i = 0; i < tris.size(); ++ i)
73 obj << 3 <<
' ' << 3*i <<
' ' << 3*i + 1 <<
' ' << 3*i + 2 <<
endl;
78 template<
class SourcePatch,
class TargetPatch>
86 OFstream obj(typeName +
"_" + name +
".obj");
88 for (
label i = 0; i < f.size(); ++ i)
90 const vector& x = ps[f[i]];
91 obj <<
"v " << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
95 for (
label i = 0; i < f.size(); ++ i)
103 template<
class SourcePatch,
class TargetPatch>
107 const FixedList<point, 4>& srcTri,
108 const FixedList<point, 4>& srcNrm
112 for (
label i = 0; i < srcN - 1; ++ i)
114 l =
max(l,
mag(srcTri[i] - srcTri[i + 1]));
117 const label nu = 20, nv = 20;
118 const scalar u0 = 0, u1 = 1, v0 = -l, v1 = l;
120 OFstream obj(typeName +
"_projection.obj");
122 for (
label i = 0; i < srcN; ++ i)
124 const point& p0 = srcTri[i], p1 = srcTri[(i + 1) % srcN];
125 const vector& n0 = srcNrm[i], n1 = srcNrm[(i + 1) % srcN];
127 for (
label iu = 0; iu <=
nu; ++ iu)
129 const scalar u = u0 + (u1 -
u0)*scalar(iu)/
nu;
130 for (
label iv = 0; iv <= nv; ++ iv)
132 const scalar v = v0 + (v1 - v0)*scalar(iv)/nv;
133 const vector x = p0 + (p1 - p0)*u + (n0 + (n1 - n0)*u)*v;
134 obj <<
"v " << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
139 for (
label i = 0; i < srcN; ++ i)
141 for (
label iu = 0; iu <
nu; ++ iu)
143 for (
label iv = 0; iv < nv; ++ iv)
146 << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 1 <<
' ' 147 << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 2 <<
' ' 148 << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 2 <<
' ' 149 << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 1<<
endl;
156 template<
class SourcePatch,
class TargetPatch>
160 FixedList<point, 4>& srcTri,
161 FixedList<point, 4>& srcNrm,
162 const FixedList<point, 3>& tgtTri
166 const vector np = plane(tgtTri[0], tgtTri[1], tgtTri[2]).normal();
169 FixedList<scalar, 4> dots;
170 for (
label i = 0; i < 3; ++ i)
172 dots[i] = srcNrm[i] & np;
178 label nR = 0, iR = -1, iF = -1;
179 for (
label i = 0; i < 3; ++ i)
181 if (dots[i] > maxDot_)
194 if (nR == 0 || nR == 3)
216 for (
label j = 2; j >= iR; -- j)
218 srcTri[j + 1] = srcTri[j];
219 srcNrm[j + 1] = srcNrm[j];
220 dots[j + 1] = dots[j];
223 const label jR = iR + 1;
224 const label i0 = (iR + 3) % 4,
j1 = (jR + 1) % 4;
226 const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]);
227 const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[
j1]);
229 srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi;
230 srcTri[jR] = srcTri[jR] + (srcTri[
j1] - srcTri[jR])*wj;
232 srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi;
233 srcNrm[jR] = srcNrm[jR] + (srcNrm[
j1] - srcNrm[jR])*wj;
254 const label iR = (iF + 1) % 3, jR = (iF + 2) % 3;
255 const label i0 = (iR + 2) % 3,
j1 = (jR + 1) % 3;
257 const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]);
258 const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[
j1]);
260 srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi;
261 srcTri[jR] = srcTri[jR] + (srcTri[
j1] - srcTri[jR])*wj;
263 srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi;
264 srcNrm[jR] = srcNrm[jR] + (srcNrm[
j1] - srcNrm[jR])*wj;
274 template<
class SourcePatch,
class TargetPatch>
281 const FixedList<point, 3>& tgtTri
285 const plane tgtPln(tgtTri[0], tgtTri[1], tgtTri[2]);
286 const point& pp = tgtPln.refPoint();
287 const vector& np = tgtPln.normal();
292 const scalar v0 = ((pp - p0) & np)/(n0 & np);
293 const scalar v1 = ((pp - p1) & np)/(n1 & np);
294 Pair<point> cutPnts(p0 + v0*n0, p1 + v1*n1);
295 Pair<vector> cutNrms((p1 - p0 + n1*v0) ^ n0, (p1 - p0 - n0*v1) ^ n1);
298 for (
label i = 0; i < 3; ++ i)
301 const point& l0 = tgtTri[i], l1 = tgtTri[(i + 1)%3];
305 const vector a = l0 - l1, b = p1 - p0, c = n0, d = n1 - n0;
306 const vector ka = k ^ a, ba = b ^ a, ca = c ^ a, da = d ^ a;
307 const scalar A = d & ba;
308 const scalar B = (c & ba) - (d & ka);
309 const scalar C = - c & ka;
312 const Roots<2> us = quadraticEqn(A, B, C).roots();
315 Pair<bool> cuts(
false,
false);
316 for (
label j = 0; j < 2; ++ j)
320 const vector den = ca + da*us[j];
323 const vector vNum = ka - ba*us[j];
324 const vector wNum = (- k + b*us[j]) ^ (c + d*us[j]);
325 vs[j] = (vNum & den)/
magSqr(den);
326 ws[j] = (wNum & den)/
magSqr(den);
327 ns[j] = (b ^ (c + d*us[j])) + (n1 ^ n0)*vs[j];
328 const bool cutu = 0 < us[j] && us[j] < 1;
329 const bool cutw = 0 < ws[j] && ws[j] < 1;
330 cuts[j] = cutu && cutw;
337 if (cuts[0] != cuts[1])
339 const label j = cuts[0] ? 0 : 1;
340 const scalar na = ns[j] & a;
341 cutPnts[na < 0] = l0 + (l1 - l0)*ws[j];
342 cutNrms[na < 0] = ns[j];
350 const vector cutDelta = cutPnts[1] - cutPnts[0];
351 const bool coincident =
mag(cutDelta) < minCutRatio_*
mag(p1 - p0);
352 return plane(cutPnts[0], coincident ? cutNrms[0] : np ^ cutDelta);
358 template<
class SourcePatch,
class TargetPatch>
361 const label srcFacei,
365 const label debugTgtFace =
366 (- 1 - debug) % this->tgtPatch_.size() == tgtFacei
367 ? (- 1 - debug) / this->tgtPatch_.size() + 1 : 0;
368 const bool debugPrint = debug > 0 || debugTgtFace > 0;
369 const bool debugWrite = debug > 1 || debugTgtFace > 1;
373 Info<<
"Inter area between source face #" << srcFacei
374 <<
" and target face #" << tgtFacei << (debugWrite ?
"\n" :
": ");
378 const pointField& srcPoints = this->srcPatch_.localPoints();
379 const pointField& tgtPoints = this->tgtPatch_.localPoints();
380 const vectorField& srcPointNormals = this->srcPatch_.pointNormals();
383 const face& srcFace = this->srcPatch_.localFaces()[srcFacei];
384 const face& tgtFace = this->tgtPatch_.localFaces()[tgtFacei];
389 writeFaceOBJ(srcFace, srcPoints,
"source");
390 writeFaceOBJ(tgtFace, tgtPoints,
"target");
394 writeFaceOBJ(tgtFace, tgtPoints,
"target" +
name(tgtFacei));
400 faceAreaIntersect::triangulate(srcFace, srcPoints, triMode, srcFaceTris);
401 faceAreaIntersect::triangulate(tgtFace, tgtPoints, triMode, tgtFaceTris);
404 scalar areaMag =
Zero;
412 tgtPoints[(*tgtIter)[0]],
413 tgtPoints[(*tgtIter)[this->reverseTarget_ ? 2 : 1]],
414 tgtPoints[(*tgtIter)[this->reverseTarget_ ? 1 : 2]]
423 srcPoints[(*srcIter)[0]],
424 srcPoints[(*srcIter)[1]],
425 srcPoints[(*srcIter)[2]],
431 srcPointNormals[(*srcIter)[0]],
432 srcPointNormals[(*srcIter)[1]],
433 srcPointNormals[(*srcIter)[2]],
438 const label srcN = getSourceProjection(srcTri, srcNrm, tgtTri);
447 writeProjectionOBJ(srcN, srcTri, srcNrm);
451 cutTriList<8> cutTris;
452 cutTris.append(tgtTri);
457 writeCutTrisVTK(cutTris,
"tris0");
464 i < srcN - 1 && (debugWrite || cutTris.size());
469 const plane cutPlane =
479 cutTriList<8> cutTrisTmp;
481 for (
label j = 0; j < cutTris.size(); ++j)
491 Swap(cutTris, cutTrisTmp);
496 writeCutTrisVTK(cutTris,
"tris" +
name(i + 1));
501 const plane cutPlane =
510 cutTriList<8> cutTrisTmp;
511 for (
label i = 0; i < cutTris.size(); ++ i)
526 if (debugWrite || debugTgtFace)
537 Swap(cutTris, cutTrisTmp);
540 if (debugTgtFace && cutTris.size())
542 static label writeCutTrisVTKIndex = 0;
546 "target" +
name(tgtFacei) +
"_" 547 +
"tris" +
name(writeCutTrisVTKIndex ++)
552 writeCutTrisVTK(cutTris,
"tris" +
name(srcN));
554 Info <<
"view triangles then press ENTER to continue ...";
564 const scalar standardAreaMag =
571 Info<<
"standard=" << standardAreaMag <<
", swept=" << areaMag
579 template<
class SourcePatch,
class TargetPatch>
587 template<
class SourcePatch,
class TargetPatch>
597 template<
class SourcePatch,
class TargetPatch>
Swept face area weighted Arbitrary Mesh Interface (AMI) method.
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 face is a list of labels corresponding to mesh vertices.
A 1D vector of objects of type <T> with a fixed size <Size>.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Ostream & endl(Ostream &os)
Add newline and flush stream.
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 ...
cut::opAddResult< AboveOp, BelowOp >::type triCut(const FixedList< point, 3 > &tri, const FixedList< scalar, 3 > &level, const AboveOp &aboveOp, const BelowOp &belowOp)
Cut a triangle along the zero plane defined by the given levels. Templated.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
dimensionedScalar j1(const dimensionedScalar &ds)
Face area weighted Arbitrary Mesh Interface (AMI) method.
vectorField pointField
pointField is a vectorField.
dimensionedScalar cos(const dimensionedScalar &ds)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
virtual scalar interArea(const label srcFacei, const label tgtFacei) const
Area of intersection between source and target faces.
vector point
Point is a vector.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual ~sweptFaceAreaWeightAMI()
Destructor.