All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sweptFaceAreaWeightAMI.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "sweptFaceAreaWeightAMI.H"
27 #include "cut.H"
28 #include "linearEqn.H"
29 #include "quadraticEqn.H"
30 #include "unitConversion.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(sweptFaceAreaWeightAMI, 0);
38  addToRunTimeSelectionTable(AMIMethod, sweptFaceAreaWeightAMI, components);
39 }
40 
41 const Foam::scalar Foam::sweptFaceAreaWeightAMI::minCutRatio_ = 10*small;
42 
43 const Foam::scalar Foam::sweptFaceAreaWeightAMI::maxDot_ =
44  - cos(degToRad(89.9));
45 
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 template<unsigned Size>
50 void Foam::sweptFaceAreaWeightAMI::writeCutTrisVTK
51 (
52  const cutTriList<Size>& tris,
53  const word& name
54 ) const
55 {
56  OFstream obj(typeName + "_" + name + ".vtk");
57 
58  obj << "# vtk DataFile Version 2.0" << endl
59  << "triSurface" << endl
60  << "ASCII" << endl
61  << "DATASET POLYDATA" << endl
62  << "POINTS " << 3*tris.size() << " float" << endl;
63 
64  for (label i = 0; i < tris.size(); ++ i)
65  {
66  for (label j = 0; j < 3; ++ j)
67  {
68  const vector& x = tris[i][j];
69  obj << x.x() << ' ' << x.y() << ' ' << x.z() << endl;
70  }
71  }
72 
73  obj << "POLYGONS " << tris.size() << ' ' << 4*tris.size() << endl;
74 
75  for (label i = 0; i < tris.size(); ++ i)
76  {
77  obj << 3 << ' ' << 3*i << ' ' << 3*i + 1 << ' ' << 3*i + 2 << endl;
78  }
79 }
80 
81 
82 void Foam::sweptFaceAreaWeightAMI::writeFaceOBJ
83 (
84  const face& f,
85  const pointField& ps,
86  const string& name
87 ) const
88 {
89  OFstream obj(typeName + "_" + name + ".obj");
90 
91  for (label i = 0; i < f.size(); ++ i)
92  {
93  const vector& x = ps[f[i]];
94  obj << "v " << x.x() << ' ' << x.y() << ' ' << x.z() << endl;
95  }
96 
97  obj << 'f';
98  for (label i = 0; i < f.size(); ++ i)
99  {
100  obj << ' ' << i + 1;
101  }
102  obj << endl;
103 }
104 
105 
106 void Foam::sweptFaceAreaWeightAMI::writeProjectionOBJ
107 (
108  const label srcN,
109  const FixedList<point, 4>& srcTri,
110  const FixedList<point, 4>& srcNrm
111 ) const
112 {
113  scalar l = 0;
114  for (label i = 0; i < srcN - 1; ++ i)
115  {
116  l = max(l, mag(srcTri[i] - srcTri[i + 1]));
117  }
118 
119  const label nu = 20, nv = 20;
120  const scalar u0 = 0, u1 = 1, v0 = -l, v1 = l;
121 
122  OFstream obj(typeName + "_projection.obj");
123 
124  for (label i = 0; i < srcN; ++ i)
125  {
126  const point& p0 = srcTri[i], p1 = srcTri[(i + 1) % srcN];
127  const vector& n0 = srcNrm[i], n1 = srcNrm[(i + 1) % srcN];
128 
129  for (label iu = 0; iu <= nu; ++ iu)
130  {
131  const scalar u = u0 + (u1 - u0)*scalar(iu)/nu;
132  for (label iv = 0; iv <= nv; ++ iv)
133  {
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;
137  }
138  }
139  }
140 
141  for (label i = 0; i < srcN; ++ i)
142  {
143  for (label iu = 0; iu < nu; ++ iu)
144  {
145  for (label iv = 0; iv < nv; ++ iv)
146  {
147  obj << "f "
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;
152  }
153  }
154  }
155 }
156 
157 
158 Foam::label Foam::sweptFaceAreaWeightAMI::getSourceProjection
159 (
160  FixedList<point, 4>& srcTri,
161  FixedList<point, 4>& srcNrm,
162  const FixedList<point, 3>& tgtTri
163 ) const
164 {
165  // The target normal
166  const vector np = plane(tgtTri[0], tgtTri[1], tgtTri[2]).normal();
167 
168  // Dot products between the projection normals and the target plane
169  FixedList<scalar, 4> dots;
170  for (label i = 0; i < 3; ++ i)
171  {
172  dots[i] = srcNrm[i] & np;
173  }
174  dots[3] = 0;
175 
176  // Search though for the number of reversed directions and an index for
177  // both a reverse and a forward direction
178  label nR = 0, iR = -1, iF = -1;
179  for (label i = 0; i < 3; ++ i)
180  {
181  if (dots[i] > maxDot_)
182  {
183  ++ nR;
184  iR = i;
185  }
186  else
187  {
188  iF = i;
189  }
190  }
191 
192  // If all the normals hit in the forward or reverse direction then return
193  // the existing projection or no projection respectively
194  if (nR == 0 || nR == 3)
195  {
196  return 3 - nR;
197  }
198 
199  // One normal hits in the reverse direction
200  if (nR == 1)
201  {
202  /*
203  o j1
204  / \
205  / \
206  / \
207  / o jR
208  / / .
209  / / .
210  / / .
211  o - > - o . . . . (old iR)
212  i0 iR
213  */
214 
215  // Make a gap by duplicating the reverse the point
216  for (label j = 2; j >= iR; -- j)
217  {
218  srcTri[j + 1] = srcTri[j];
219  srcNrm[j + 1] = srcNrm[j];
220  dots[j + 1] = dots[j];
221  }
222 
223  const label jR = iR + 1;
224  const label i0 = (iR + 3) % 4, j1 = (jR + 1) % 4;
225 
226  const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]);
227  const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[j1]);
228 
229  srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi;
230  srcTri[jR] = srcTri[jR] + (srcTri[j1] - srcTri[jR])*wj;
231 
232  srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi;
233  srcNrm[jR] = srcNrm[jR] + (srcNrm[j1] - srcNrm[jR])*wj;
234 
235  return 4;
236  }
237 
238  // Two normals hit in the reverse direction
239  if (nR == 2)
240  {
241  /*
242  . (old jR)
243  . .
244  jR o .
245  / \ .
246  / \ .
247  / \ .
248  / \ .
249  / \ .
250  o - - > - - o . . (old iR)
251  iF iR
252  */
253 
254  const label iR = (iF + 1) % 3, jR = (iF + 2) % 3;
255  const label i0 = (iR + 2) % 3, j1 = (jR + 1) % 3;
256 
257  const scalar wi = (dots[iR] - maxDot_)/(dots[iR] - dots[i0]);
258  const scalar wj = (dots[jR] - maxDot_)/(dots[jR] - dots[j1]);
259 
260  srcTri[iR] = srcTri[iR] + (srcTri[i0] - srcTri[iR])*wi;
261  srcTri[jR] = srcTri[jR] + (srcTri[j1] - srcTri[jR])*wj;
262 
263  srcNrm[iR] = srcNrm[iR] + (srcNrm[i0] - srcNrm[iR])*wi;
264  srcNrm[jR] = srcNrm[jR] + (srcNrm[j1] - srcNrm[jR])*wj;
265 
266  return 3;
267  }
268 
269  // This cannot happen
270  return -1;
271 }
272 
273 
274 Foam::plane Foam::sweptFaceAreaWeightAMI::getCutPlane
275 (
276  const point& p0,
277  const point& p1,
278  const vector& n0,
279  const vector& n1,
280  const FixedList<point, 3>& tgtTri
281 ) const
282 {
283  // The target plane
284  const plane tgtPln(tgtTri[0], tgtTri[1], tgtTri[2]);
285  const point& pp = tgtPln.refPoint();
286  const vector& np = tgtPln.normal();
287 
288  // Calculate the bounding intersection points. These are the locations at
289  // which the bounding lines of the projected surface intersect with the
290  // target plane.
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);
295 
296  // Calculate edge intersections with the surface
297  for (label i = 0; i < 3; ++ i)
298  {
299  // The current target edge
300  const point& l0 = tgtTri[i], l1 = tgtTri[(i + 1)%3];
301 
302  // Form the quadratic in the surface parameter u
303  const vector k = l0 - p0;
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;
309 
310  // Solve and extract the other parameters
311  const Roots<2> us = quadraticEqn(A, B, C).roots();
312  Pair<scalar> vs, ws;
313  Pair<vector> ns;
314  Pair<bool> cuts(false, false);
315  for (label j = 0; j < 2; ++ j)
316  {
317  if (us.type(j) == rootType::real)
318  {
319  const vector den = ca + da*us[j];
320  if (magSqr(den) > vSmall)
321  {
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;
330  }
331  }
332  }
333 
334  // If we have just one intersection in bounds then store it in the
335  // result list based on its direction
336  if (cuts[0] != cuts[1])
337  {
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];
342  }
343  }
344 
345  // Generate and return the cut plane. If the cut points are not coincident
346  // then form a plane normal by crossing the displacement between the points
347  // by the target plane normal. If the points are coincident then use the
348  // projected surface normal evaluated at the first cut point.
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);
352 };
353 
354 
355 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
356 
358 (
359  const label srcFacei,
360  const label tgtFacei
361 ) const
362 {
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;
368 
369  if (debugPrint)
370  {
371  Info<< "Inter area between source face #" << srcFacei
372  << " and target face #" << tgtFacei << (debugWrite ? "\n" : ": ");
373  }
374 
375  // Patch data
376  const pointField& srcPoints = this->srcPatch_.localPoints();
377  const pointField& tgtPoints = this->tgtPatch_.localPoints();
378  const vectorField& srcPointNormals = this->srcPatch_.pointNormals();
379 
380  // Faces
381  const face& srcFace = this->srcPatch_.localFaces()[srcFacei];
382  const face& tgtFace = this->tgtPatch_.localFaces()[tgtFacei];
383 
384  // Write out the faces
385  if (debugWrite)
386  {
387  writeFaceOBJ(srcFace, srcPoints, "source");
388  writeFaceOBJ(tgtFace, tgtPoints, "target");
389  }
390  if (debugTgtFace)
391  {
392  writeFaceOBJ(tgtFace, tgtPoints, "target" + name(tgtFacei));
393  }
394 
395  // Triangulate the faces
396  const faceAreaIntersect::triangulationMode triMode = this->triMode_;
397  faceList srcFaceTris, tgtFaceTris;
398  faceAreaIntersect::triangulate(srcFace, srcPoints, triMode, srcFaceTris);
399  faceAreaIntersect::triangulate(tgtFace, tgtPoints, triMode, tgtFaceTris);
400 
401  // Area sum
402  scalar areaMag = Zero;
403 
404  // Loop the target triangles
405  forAllConstIter(faceList, tgtFaceTris, tgtIter)
406  {
407  const FixedList<point, 3>
408  tgtTri =
409  {
410  tgtPoints[(*tgtIter)[0]],
411  tgtPoints[(*tgtIter)[this->reverseTarget_ ? 2 : 1]],
412  tgtPoints[(*tgtIter)[this->reverseTarget_ ? 1 : 2]]
413  };
414 
415  // Loop the source triangles
416  forAllConstIter(faceList, srcFaceTris, srcIter)
417  {
419  srcTri =
420  {
421  srcPoints[(*srcIter)[0]],
422  srcPoints[(*srcIter)[1]],
423  srcPoints[(*srcIter)[2]],
425  };
427  srcNrm =
428  {
429  srcPointNormals[(*srcIter)[0]],
430  srcPointNormals[(*srcIter)[1]],
431  srcPointNormals[(*srcIter)[2]],
432  vector::zero
433  };
434 
435  // Get the source projection modified for any reverse intersections
436  const label srcN = getSourceProjection(srcTri, srcNrm, tgtTri);
437  if (srcN <= 0)
438  {
439  continue;
440  }
441 
442  // Write the source projection
443  if (debugWrite)
444  {
445  writeProjectionOBJ(srcN, srcTri, srcNrm);
446  }
447 
448  // Create the initial cut triangle list
449  cutTriList<8> cutTris;
450  cutTris.append(tgtTri);
451 
452  // Write the initial triangle
453  if (debugWrite)
454  {
455  writeCutTrisVTK(cutTris, "tris0");
456  }
457 
458  // Do all but one of the cuts
459  for
460  (
461  label i = 0;
462  i < srcN - 1 && (debugWrite || cutTris.size());
463  ++ i
464  )
465  {
466  // Do the cut
467  const plane cutPlane =
468  getCutPlane
469  (
470  srcTri[i],
471  srcTri[i+1],
472  srcNrm[i],
473  srcNrm[i+1],
474  tgtTri
475  );
476 
477  cutTriList<8> cutTrisTmp;
478 
479  for (label j = 0; j < cutTris.size(); ++j)
480  {
481  triCut
482  (
483  cutTris[j],
484  cutPlane,
485  cut::noOp(),
486  cut::appendOp<cutTriList<8>>(cutTrisTmp)
487  );
488  }
489  Swap(cutTris, cutTrisTmp);
490 
491  // Write the triangles resulting from the cut
492  if (debugWrite)
493  {
494  writeCutTrisVTK(cutTris, "tris" + name(i + 1));
495  }
496  }
497 
498  // Do the last cut
499  const plane cutPlane =
500  getCutPlane
501  (
502  srcTri[srcN - 1],
503  srcTri[0],
504  srcNrm[srcN - 1],
505  srcNrm[0],
506  tgtTri
507  );
508  cutTriList<8> cutTrisTmp;
509  for (label i = 0; i < cutTris.size(); ++ i)
510  {
511  // Sum the area of the cut triangle
512  areaMag +=
513  mag
514  (
515  triCut
516  (
517  cutTris[i],
518  cutPlane,
519  cut::noOp(),
520  cut::areaOp()
521  )
522  );
523  // Store the cut triangle if it is needed for output
524  if (debugWrite || debugTgtFace)
525  {
526  triCut
527  (
528  cutTris[i],
529  cutPlane,
530  cut::noOp(),
531  cut::appendOp<cutTriList<8>>(cutTrisTmp)
532  );
533  }
534  }
535  Swap(cutTris, cutTrisTmp);
536 
537  // Write the triangles resulting from the cuts
538  if (debugTgtFace && cutTris.size())
539  {
540  static label writeCutTrisVTKIndex = 0;
541  writeCutTrisVTK
542  (
543  cutTris,
544  "target" + name(tgtFacei) + "_"
545  + "tris" + name(writeCutTrisVTKIndex ++)
546  );
547  }
548  if (debugWrite)
549  {
550  writeCutTrisVTK(cutTris, "tris" + name(srcN));
551 
552  Info << "view triangles then press ENTER to continue ...";
553  getchar();
554  }
555  }
556  }
557 
558  // Print the difference between this inter-area and that obtained by the
559  // standard algorithm built around faceAreaIntersect
560  if (debugPrint)
561  {
562  const scalar standardAreaMag =
564  (
565  srcFacei,
566  tgtFacei
567  );
568 
569  Info<<"standard=" << standardAreaMag << ", swept=" << areaMag
570  << endl;
571  }
572 
573  return areaMag;
574 }
575 
576 
577 Foam::scalar Foam::sweptFaceAreaWeightAMI::minWeight() const
578 {
579  return small;
580 }
581 
582 
583 Foam::scalar Foam::sweptFaceAreaWeightAMI::maxWalkAngle() const
584 {
585  return degToRad(180);
586 }
587 
588 
589 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
590 
592 {}
593 
594 
595 // ************************************************************************* //
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.
Definition: label.H:59
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Unit conversion functions.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
Macros for easy insertion into run-time selection tables.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
dimensionedScalar j1(const dimensionedScalar &ds)
scalar u0
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dimensionedScalar cos(const dimensionedScalar &ds)
void Swap(T &a, T &b)
Definition: Swap.H:43
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.
static const zero Zero
Definition: zero.H:97
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
const Cmpt & x() const
Definition: VectorI.H:75
dimensioned< scalar > magSqr(const dimensioned< Type > &)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
virtual scalar interArea(const label srcFacei, const label tgtFacei) const
Area of intersection between source and target faces.
vector point
Point is a vector.
Definition: point.H:41
messageStream Info
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.
Namespace for OpenFOAM.
virtual ~sweptFaceAreaWeightAMI()
Destructor.