41 const Foam::scalar Foam::sweptFaceAreaWeightAMI::minCutRatio_ = 10*small;
43 const Foam::scalar Foam::sweptFaceAreaWeightAMI::maxDot_ =
49 template<
unsigned Size>
50 void Foam::sweptFaceAreaWeightAMI::writeCutTrisVTK
52 const cutTriList<Size>& tris,
56 OFstream obj(typeName +
"_" + name +
".vtk");
58 obj <<
"# vtk DataFile Version 2.0" <<
endl 59 <<
"triSurface" <<
endl 61 <<
"DATASET POLYDATA" <<
endl 62 <<
"POINTS " << 3*tris.size() <<
" float" <<
endl;
64 for (
label i = 0; i < tris.size(); ++ i)
66 for (
label j = 0; j < 3; ++ j)
68 const vector& x = tris[i][j];
69 obj << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
73 obj <<
"POLYGONS " << tris.size() <<
' ' << 4*tris.size() <<
endl;
75 for (
label i = 0; i < tris.size(); ++ i)
77 obj << 3 <<
' ' << 3*i <<
' ' << 3*i + 1 <<
' ' << 3*i + 2 <<
endl;
82 void Foam::sweptFaceAreaWeightAMI::writeFaceOBJ
89 OFstream obj(typeName +
"_" + name +
".obj");
91 for (
label i = 0; i < f.size(); ++ i)
93 const vector& x = ps[f[i]];
94 obj <<
"v " << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
98 for (
label i = 0; i < f.size(); ++ i)
106 void Foam::sweptFaceAreaWeightAMI::writeProjectionOBJ
109 const FixedList<point, 4>& srcTri,
110 const FixedList<point, 4>& srcNrm
114 for (
label i = 0; i < srcN - 1; ++ i)
116 l =
max(l,
mag(srcTri[i] - srcTri[i + 1]));
119 const label nu = 20, nv = 20;
120 const scalar u0 = 0, u1 = 1, v0 = -l, v1 = l;
122 OFstream obj(typeName +
"_projection.obj");
124 for (
label i = 0; i < srcN; ++ i)
126 const point& p0 = srcTri[i], p1 = srcTri[(i + 1) % srcN];
127 const vector& n0 = srcNrm[i], n1 = srcNrm[(i + 1) % srcN];
129 for (
label iu = 0; iu <= nu; ++ iu)
131 const scalar u = u0 + (u1 -
u0)*scalar(iu)/nu;
132 for (
label iv = 0; iv <= nv; ++ iv)
134 const scalar v = v0 + (v1 - v0)*scalar(iv)/nv;
135 const vector x = p0 + (p1 - p0)*u + (n0 + (n1 - n0)*u)*v;
136 obj <<
"v " << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
141 for (
label i = 0; i < srcN; ++ i)
143 for (
label iu = 0; iu < nu; ++ iu)
145 for (
label iv = 0; iv < nv; ++ iv)
148 << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 1 <<
' ' 149 << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 2 <<
' ' 150 << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 2 <<
' ' 151 << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 1<<
endl;
158 Foam::label Foam::sweptFaceAreaWeightAMI::getSourceProjection
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 Foam::plane Foam::sweptFaceAreaWeightAMI::getCutPlane
280 const FixedList<point, 3>& tgtTri
284 const plane tgtPln(tgtTri[0], tgtTri[1], tgtTri[2]);
285 const point& pp = tgtPln.refPoint();
286 const vector& np = tgtPln.normal();
291 const scalar v0 = ((pp - p0) & np)/(n0 & np);
292 const scalar v1 = ((pp - p1) & np)/(n1 & np);
293 Pair<point> cutPnts(p0 + v0*n0, p1 + v1*n1);
294 Pair<vector> cutNrms((p1 - p0 + n1*v0) ^ n0, (p1 - p0 - n0*v1) ^ n1);
297 for (
label i = 0; i < 3; ++ i)
300 const point& l0 = tgtTri[i], l1 = tgtTri[(i + 1)%3];
304 const vector a = l0 - l1, b = p1 - p0, c = n0, d = n1 - n0;
305 const vector ka = k ^ a, ba = b ^ a, ca = c ^ a, da = d ^ a;
306 const scalar A = d & ba;
307 const scalar B = (c & ba) - (d & ka);
308 const scalar C = - c & ka;
311 const Roots<2> us = quadraticEqn(A, B, C).roots();
314 Pair<bool> cuts(
false,
false);
315 for (
label j = 0; j < 2; ++ j)
319 const vector den = ca + da*us[j];
322 const vector vNum = ka - ba*us[j];
323 const vector wNum = (- k + b*us[j]) ^ (c + d*us[j]);
324 vs[j] = (vNum & den)/
magSqr(den);
325 ws[j] = (wNum & den)/
magSqr(den);
326 ns[j] = (b ^ (c + d*us[j])) + (n1 ^ n0)*vs[j];
327 const bool cutu = 0 < us[j] && us[j] < 1;
328 const bool cutw = 0 < ws[j] && ws[j] < 1;
329 cuts[j] = cutu && cutw;
336 if (cuts[0] != cuts[1])
338 const label j = cuts[0] ? 0 : 1;
339 const scalar na = ns[j] & a;
340 cutPnts[na < 0] = l0 + (l1 - l0)*ws[j];
341 cutNrms[na < 0] = ns[j];
349 const vector cutDelta = cutPnts[1] - cutPnts[0];
350 const bool coincident =
mag(cutDelta) < minCutRatio_*
mag(p1 - p0);
351 return plane(cutPnts[0], coincident ? cutNrms[0] : np ^ cutDelta);
359 const label srcFacei,
363 const label debugTgtFace =
364 (- 1 - debug) % this->tgtPatch_.size() == tgtFacei
365 ? (- 1 - debug) / this->tgtPatch_.size() + 1 : 0;
366 const bool debugPrint = debug > 0 || debugTgtFace > 0;
367 const bool debugWrite = debug > 1 || debugTgtFace > 1;
371 Info<<
"Inter area between source face #" << srcFacei
372 <<
" and target face #" << tgtFacei << (debugWrite ?
"\n" :
": ");
376 const pointField& srcPoints = this->srcPatch_.localPoints();
377 const pointField& tgtPoints = this->tgtPatch_.localPoints();
378 const vectorField& srcPointNormals = this->srcPatch_.pointNormals();
381 const face& srcFace = this->srcPatch_.localFaces()[srcFacei];
382 const face& tgtFace = this->tgtPatch_.localFaces()[tgtFacei];
387 writeFaceOBJ(srcFace, srcPoints,
"source");
388 writeFaceOBJ(tgtFace, tgtPoints,
"target");
392 writeFaceOBJ(tgtFace, tgtPoints,
"target" +
name(tgtFacei));
402 scalar areaMag =
Zero;
410 tgtPoints[(*tgtIter)[0]],
411 tgtPoints[(*tgtIter)[this->reverseTarget_ ? 2 : 1]],
412 tgtPoints[(*tgtIter)[this->reverseTarget_ ? 1 : 2]]
421 srcPoints[(*srcIter)[0]],
422 srcPoints[(*srcIter)[1]],
423 srcPoints[(*srcIter)[2]],
429 srcPointNormals[(*srcIter)[0]],
430 srcPointNormals[(*srcIter)[1]],
431 srcPointNormals[(*srcIter)[2]],
436 const label srcN = getSourceProjection(srcTri, srcNrm, tgtTri);
445 writeProjectionOBJ(srcN, srcTri, srcNrm);
449 cutTriList<8> cutTris;
450 cutTris.append(tgtTri);
455 writeCutTrisVTK(cutTris,
"tris0");
462 i < srcN - 1 && (debugWrite || cutTris.size());
467 const plane cutPlane =
477 cutTriList<8> cutTrisTmp;
479 for (
label j = 0; j < cutTris.size(); ++j)
489 Swap(cutTris, cutTrisTmp);
494 writeCutTrisVTK(cutTris,
"tris" +
name(i + 1));
499 const plane cutPlane =
508 cutTriList<8> cutTrisTmp;
509 for (
label i = 0; i < cutTris.size(); ++ i)
524 if (debugWrite || debugTgtFace)
535 Swap(cutTris, cutTrisTmp);
538 if (debugTgtFace && cutTris.size())
540 static label writeCutTrisVTKIndex = 0;
544 "target" +
name(tgtFacei) +
"_" 545 +
"tris" +
name(writeCutTrisVTKIndex ++)
550 writeCutTrisVTK(cutTris,
"tris" +
name(srcN));
552 Info <<
"view triangles then press ENTER to continue ...";
562 const scalar standardAreaMag =
569 Info<<
"standard=" << standardAreaMag <<
", swept=" << areaMag
577 Foam::scalar Foam::sweptFaceAreaWeightAMI::minWeight()
const 583 Foam::scalar Foam::sweptFaceAreaWeightAMI::maxWalkAngle()
const virtual scalar interArea(const label srcFacei, const label tgtFacei) const
Area of intersection between source and target faces.
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 > &)
Unit conversion functions.
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 ...
Macros for easy insertion into run-time selection tables.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
dimensionedScalar j1(const dimensionedScalar &ds)
vectorField pointField
pointField is a vectorField.
dimensionedScalar cos(const dimensionedScalar &ds)
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.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
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 > &)
static void triangulate(const face &f, const pointField &points, const triangulationMode &triMode, faceList &faceTris)
Triangulate a face using the given triangulation mode.
static const Vector< scalar > zero
virtual ~sweptFaceAreaWeightAMI()
Destructor.