42 const Foam::scalar Foam::sweptFaceAreaWeightAMI::minCutRatio_ = 10*small;
44 const Foam::scalar Foam::sweptFaceAreaWeightAMI::maxDot_ =
50 template<
unsigned Size>
51 void Foam::sweptFaceAreaWeightAMI::writeCutTrisVTK
53 const cutTriList<Size>& tris,
57 OFstream obj(typeName +
"_" + name +
".vtk");
59 obj <<
"# vtk DataFile Version 2.0" <<
endl 60 << triSurface::typeName <<
endl 62 <<
"DATASET POLYDATA" <<
endl 63 <<
"POINTS " << 3*tris.size() <<
" float" <<
endl;
65 for (
label i = 0; i < tris.size(); ++ i)
67 for (
label j = 0; j < 3; ++ j)
69 const vector& x = tris[i][j];
70 obj << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
74 obj <<
"POLYGONS " << tris.size() <<
' ' << 4*tris.size() <<
endl;
76 for (
label i = 0; i < tris.size(); ++ i)
78 obj << 3 <<
' ' << 3*i <<
' ' << 3*i + 1 <<
' ' << 3*i + 2 <<
endl;
83 void Foam::sweptFaceAreaWeightAMI::writeFaceOBJ
90 OFstream obj(typeName +
"_" + name +
".obj");
92 for (
label i = 0; i < f.size(); ++ i)
94 const vector& x = ps[f[i]];
95 obj <<
"v " << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
99 for (
label i = 0; i < f.size(); ++ i)
107 void Foam::sweptFaceAreaWeightAMI::writeProjectionOBJ
110 const FixedList<point, 4>& srcTri,
111 const FixedList<point, 4>& srcNrm
115 for (
label i = 0; i < srcN - 1; ++ i)
117 l =
max(l,
mag(srcTri[i] - srcTri[i + 1]));
120 const label nu = 20, nv = 20;
121 const scalar u0 = 0, u1 = 1, v0 = -l, v1 = l;
123 OFstream obj(typeName +
"_projection.obj");
125 for (
label i = 0; i < srcN; ++ i)
127 const point& p0 = srcTri[i], p1 = srcTri[(i + 1) % srcN];
128 const vector& n0 = srcNrm[i], n1 = srcNrm[(i + 1) % srcN];
130 for (
label iu = 0; iu <= nu; ++ iu)
132 const scalar u = u0 + (u1 -
u0)*scalar(iu)/nu;
133 for (
label iv = 0; iv <= nv; ++ iv)
135 const scalar v = v0 + (v1 - v0)*scalar(iv)/nv;
136 const vector x = p0 + (p1 - p0)*u + (n0 + (n1 - n0)*u)*v;
137 obj <<
"v " << x.
x() <<
' ' << x.y() <<
' ' << x.z() <<
endl;
142 for (
label i = 0; i < srcN; ++ i)
144 for (
label iu = 0; iu < nu; ++ iu)
146 for (
label iv = 0; iv < nv; ++ iv)
149 << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 1 <<
' ' 150 << i*(nu + 1)*(nv + 1) + (nv + 1)*iu + iv + 2 <<
' ' 151 << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 2 <<
' ' 152 << i*(nu + 1)*(nv + 1) + (nv + 1)*(iu + 1) + iv + 1<<
endl;
159 Foam::label Foam::sweptFaceAreaWeightAMI::getSourceProjection
161 FixedList<point, 4>& srcTri,
162 FixedList<point, 4>& srcNrm,
163 const FixedList<point, 3>& tgtTri
167 const vector np = plane(tgtTri[0], tgtTri[1], tgtTri[2]).normal();
170 FixedList<scalar, 4> dots;
171 for (
label i = 0; i < 3; ++ i)
173 dots[i] = srcNrm[i] & np;
179 label nR = 0, iR = -1, iF = -1;
180 for (
label i = 0; i < 3; ++ i)
182 if (dots[i] > maxDot_)
195 if (nR == 0 || nR == 3)
217 for (
label j = 2; j >= iR; -- j)
219 srcTri[j + 1] = srcTri[j];
220 srcNrm[j + 1] = srcNrm[j];
221 dots[j + 1] = dots[j];
224 const label jR = iR + 1;
225 const label i0 = (iR + 3) % 4,
j1 = (jR + 1) % 4;
227 const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]);
228 const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[
j1]);
230 srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi;
231 srcTri[jR] = srcTri[jR] + (srcTri[
j1] - srcTri[jR])*wj;
233 srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi;
234 srcNrm[jR] = srcNrm[jR] + (srcNrm[
j1] - srcNrm[jR])*wj;
255 const label iR = (iF + 1) % 3, jR = (iF + 2) % 3;
256 const label i0 = (iR + 2) % 3,
j1 = (jR + 1) % 3;
258 const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]);
259 const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[
j1]);
261 srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi;
262 srcTri[jR] = srcTri[jR] + (srcTri[
j1] - srcTri[jR])*wj;
264 srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi;
265 srcNrm[jR] = srcNrm[jR] + (srcNrm[
j1] - srcNrm[jR])*wj;
275 Foam::plane Foam::sweptFaceAreaWeightAMI::getCutPlane
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);
360 const label srcFacei,
364 const label debugTgtFace =
365 (- 1 - debug) % this->tgtPatch_.size() == tgtFacei
366 ? (- 1 - debug) / this->tgtPatch_.size() + 1 : 0;
367 const bool debugPrint = debug > 0 || debugTgtFace > 0;
368 const bool debugWrite = debug > 1 || debugTgtFace > 1;
372 Info<<
"Inter area between source face #" << srcFacei
373 <<
" and target face #" << tgtFacei << (debugWrite ?
"\n" :
": ");
377 const pointField& srcPoints = this->srcPatch_.localPoints();
378 const pointField& tgtPoints = this->tgtPatch_.localPoints();
379 const vectorField& srcPointNormals = this->srcPatch_.pointNormals();
382 const face& srcFace = this->srcPatch_.localFaces()[srcFacei];
383 const face& tgtFace = this->tgtPatch_.localFaces()[tgtFacei];
388 writeFaceOBJ(srcFace, srcPoints,
"source");
389 writeFaceOBJ(tgtFace, tgtPoints,
"target");
393 writeFaceOBJ(tgtFace, tgtPoints,
"target" +
name(tgtFacei));
403 scalar areaMag =
Zero;
411 tgtPoints[(*tgtIter)[0]],
412 tgtPoints[(*tgtIter)[this->reverseTarget_ ? 2 : 1]],
413 tgtPoints[(*tgtIter)[this->reverseTarget_ ? 1 : 2]]
422 srcPoints[(*srcIter)[0]],
423 srcPoints[(*srcIter)[1]],
424 srcPoints[(*srcIter)[2]],
430 srcPointNormals[(*srcIter)[0]],
431 srcPointNormals[(*srcIter)[1]],
432 srcPointNormals[(*srcIter)[2]],
437 const label srcN = getSourceProjection(srcTri, srcNrm, tgtTri);
446 writeProjectionOBJ(srcN, srcTri, srcNrm);
450 cutTriList<8> cutTris;
451 cutTris.append(tgtTri);
456 writeCutTrisVTK(cutTris,
"tris0");
463 i < srcN - 1 && (debugWrite || cutTris.size());
468 const plane cutPlane =
478 cutTriList<8> cutTrisTmp;
480 for (
label j = 0; j < cutTris.size(); ++j)
490 Swap(cutTris, cutTrisTmp);
495 writeCutTrisVTK(cutTris,
"tris" +
name(i + 1));
500 const plane cutPlane =
509 cutTriList<8> cutTrisTmp;
510 for (
label i = 0; i < cutTris.size(); ++ i)
525 if (debugWrite || debugTgtFace)
536 Swap(cutTris, cutTrisTmp);
539 if (debugTgtFace && cutTris.size())
541 static label writeCutTrisVTKIndex = 0;
545 "target" +
name(tgtFacei) +
"_" 546 +
"tris" +
name(writeCutTrisVTKIndex ++)
551 writeCutTrisVTK(cutTris,
"tris" +
name(srcN));
553 Info <<
"view triangles then press ENTER to continue ...";
563 const scalar standardAreaMag =
570 Info<<
"standard=" << standardAreaMag <<
", swept=" << areaMag
578 Foam::scalar Foam::sweptFaceAreaWeightAMI::minWeight()
const 584 Foam::scalar Foam::sweptFaceAreaWeightAMI::maxWalkAngle()
const virtual scalar interArea(const label srcFacei, const label tgtFacei) const
Area of intersection between source and target faces.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > 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>.
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of 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 ...
static void triangulate(const face &f, const pointField &points, const triangulationMode &triMode, triFaceList &faceTris)
Triangulate a face using the given triangulation mode.
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.
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 const Vector< scalar > zero
virtual ~sweptFaceAreaWeightAMI()
Destructor.