MPLICcell.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) 2020-2022 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 "MPLICcell.H"
27 #include "tetCell.H"
28 #include "cubicEqn.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 Foam::label Foam::MPLICcell::calcMatchAlphaCutCell
33 (
34  const MPLICcellStorage& cellInfo,
35  const bool tetDecom
36 )
37 {
38  // Clear all the temporary lists/fields
39  clear();
40 
41  const scalar cellAlpha = cellInfo.cellAlpha();
42 
43  const UIndirectList<scalar>& cellMagSfs = cellInfo.magSf();
44 
45  // Initialise fluxes from velocity point values
46  if (!unweighted_)
47  {
48  phiU
49  (
50  cellInfo.points(),
51  cellInfo.faces(),
52  cellInfo.cellFaces(),
53  cellInfo.pointsU()
54  );
55  }
56 
57  // Difference between point alpha values isn't large enough
58  if (mag(cellInfo.cellAlphaMax() - cellInfo.cellAlphaMin()) < vSmall)
59  {
60  return -1;
61  }
62 
63  // Finding between which two points cell cut lies
64  findPointAlphaBounds(cellInfo, tetDecom);
65 
66  // Direction of the cut isn't clear => use default scheme
67  if (mag(pCubicAlphas_.a() - pCubicAlphas_.d()) < vSmall)
68  {
69  return -1;
70  }
71 
72  // Calculates 2 intermediate volume fractions necessary for cubic polyfit
73  calcPointAlphaInterior(cellInfo, tetDecom);
74 
75  // Calculate coefficients of cubic polynomial (computed on normalised input)
76  const FixedList<scalar, 4> coeffs(solveVanderMatrix());
77 
78  // Direct solution of cubic equation and roots selection
79  findRoots(cellInfo, coeffs, tetDecom);
80 
81  // Use default scheme if the calculated alpha is more than 10% off from
82  // original cell alpha.
83  if (mag(cutAlpha_) > rootSmall && (1 - mag(cellAlpha/cutAlpha_)) > 0.1)
84  {
85  return 0;
86  }
87  else
88  {
89  unweighted_ ? calcAlphaf(cellMagSfs) : calcAlphaUf();
90  return 1;
91  }
92 }
93 
94 
95 void Foam::MPLICcell::findPointAlphaBounds
96 (
97  const MPLICcellStorage& cellInfo,
98  const bool tetDecom
99 )
100 {
101  const scalar cellAlpha = cellInfo.cellAlpha();
102 
103  cellPointsAlpha_ =
104  UIndirectList<scalar>(cellInfo.pointsAlpha(), cellInfo.cellPoints());
105 
106  if (tetDecom)
107  {
108  cellPointsAlpha_.append(cellAlpha);
109  }
110 
111  sort(cellPointsAlpha_);
112 
113  // Avoid useless cuts and keep only unique point values
114 
115  pointsAlpha_.clear();
116  pointsAlpha_.append(cellPointsAlpha_[0]);
117  for (label i=1; i<cellPointsAlpha_.size(); i++)
118  {
119  if (mag(cellPointsAlpha_[i-1] - cellPointsAlpha_[i]) > vSmall)
120  {
121  pointsAlpha_.append(cellPointsAlpha_[i]);
122  }
123  }
124 
125  const label nAlphas = pointsAlpha_.size();
126 
127  // If there is only one value (all the point are the same)
128  // the cut is not determined
129  if (nAlphas < 2)
130  {
131  pCubicAlphas_.a() = -1;
132  pCubicAlphas_.d() = -1;
133  cCubicAlphas_.a() = -1;
134  cCubicAlphas_.d() = -1;
135  }
136 
137  // If in the cell are at least two points with different values we can
138  // attempt to cut between them
139  else if (nAlphas == 2)
140  {
141  pCubicAlphas_.a() = pointsAlpha_.first();
142  pCubicAlphas_.d() = pointsAlpha_.last();
143  cCubicAlphas_.a() = 1;
144  cCubicAlphas_.d() = 0;
145  }
146 
147  // We know extreme of the volume fraction so for the points that have those
148  // it will be only 0 and 1, therefore we make only necessary cuts starting
149  // from the mid index point value
150  else
151  {
152  // Pick the mid point value
153  label index = label(round((nAlphas)/2.0)) - 1;
154 
155  // Calculate initial cell and point value
156  scalar target = pointsAlpha_[index];
157  scalar cutAlpha = calcAlpha(cellInfo, target, tetDecom);
158  scalar targetOld = target;
159  scalar cutAlphaOld = cutAlpha;
160 
161  for (label i = 1; i < nAlphas-1; ++i)
162  {
163  (cutAlpha >= cellAlpha) ? index++ : index--;
164 
165  // Special case
166  // Maximum point value and minimum volume fraction = 0
167  if (index == nAlphas - 1)
168  {
169  cCubicAlphas_.a() = cutAlpha;
170  cCubicAlphas_.d() = 0;
171  pCubicAlphas_.a() = target;
172  pCubicAlphas_.d() = pointsAlpha_[index];
173  break;
174  }
175 
176  // Special case
177  // Minimum point value and maximum volume fraction = 1
178  else if (index == 0)
179  {
180  cCubicAlphas_.a() = 1;
181  cCubicAlphas_.d() = cutAlpha;
182  pCubicAlphas_.a() = pointsAlpha_[index];
183  pCubicAlphas_.d() = target;
184  break;
185  }
186 
187  // Calculate new values
188  target = pointsAlpha_[index];
189  cutAlpha = calcAlpha(cellInfo, target, tetDecom);
190 
191  if (cutAlphaOld > cellAlpha && cutAlpha < cellAlpha)
192  {
193  cCubicAlphas_.a() = cutAlphaOld;
194  cCubicAlphas_.d() = cutAlpha;
195  pCubicAlphas_.a() = targetOld;
196  pCubicAlphas_.d() = target;
197  break;
198  }
199 
200  else if (cutAlphaOld < cellAlpha && cutAlpha > cellAlpha)
201  {
202  cCubicAlphas_.a() = cutAlpha;
203  cCubicAlphas_.d() = cutAlphaOld;
204  pCubicAlphas_.a() = target;
205  pCubicAlphas_.d() = targetOld;
206  break;
207  }
208 
209  // Store previous iteration values
210  targetOld = target;
211  cutAlphaOld = cutAlpha;
212  }
213  }
214 }
215 
216 
217 void Foam::MPLICcell::calcPointAlphaInterior
218 (
219  const MPLICcellStorage& cellInfo,
220  const bool tetDecom
221 )
222 {
223  for (label i=1; i<=2; i++)
224  {
225  pCubicAlphas_[i] =
226  pCubicAlphas_.a() + (pCubicAlphas_.d() - pCubicAlphas_.a())*(i/3.0);
227 
228  cCubicAlphas_[i] = calcAlpha(cellInfo, pCubicAlphas_[i], tetDecom);
229  }
230 }
231 
232 
233 Foam::FixedList<Foam::scalar, 4> Foam::MPLICcell::solveVanderMatrix() const
234 {
235  // The cubic polynomial of the volume is fit to a normalised coordinate,
236  // which is defined as follows (see the use of sFactor in findRoots below):
237  //
238  // x = (alpha - pCubicAlpha[0])/(pCubicAlpha[3] - pCubicAlpha[0])
239  //
240  // pCubicAlpha[0] and pCubicAlpha[3] are the bounds of the fit, and are
241  // actual iso-values on points of the cell. pCubicAlpha_1 and pCubicAlpha_2
242  // are interior values needed to complete the cubic fit. The 4 values
243  // are equidistant (see calcPointAlphaInterior above) so the corresponding
244  // values of the normalised coordinate x are always:
245  //
246  // x = [0 1/3 2/3 1]
247  //
248  // This means the Vandermonde matrix that is solved for the cubic
249  // polynomial's coefficients is always the same:
250  //
251  // V = [0 0 0 1]
252  // [1/27 1/9 1/3 1]
253  // [8/27 4/9 2/3 1]
254  // [1 1 1 1]
255  //
256  // This means it's inverse can be precomputed. This pre-computation is hard
257  // coded below.
258 
259  const vector4& b = cCubicAlphas_;
260 
261  return FixedList<scalar, 4>
262  {
263  scalar(- 4.5*b[0] + 13.5*b[1] - 13.5*b[2] + 4.5*b[3]),
264  scalar(9.0*b[0] - 22.5*b[1] + 18.0*b[2] - 4.5*b[3]),
265  scalar(- 5.5*b[0] + 9.0*b[1] - 4.5*b[2] + 1.0*b[3]),
266  scalar(1.0*b[0])
267  };
268 }
269 
270 
271 void Foam::MPLICcell::findRoots
272 (
273  const MPLICcellStorage& cellInfo,
274  const FixedList<scalar, 4>& coeff,
275  const bool tetDecom
276 )
277 {
278  const scalar cellAlpha = cellInfo.cellAlpha();
279 
280  // Solve cubic polynomial exactly
281  const Roots<3> roots =
282  cubicEqn(coeff[0], coeff[1], coeff[2], coeff[3] - cellAlpha).roots();
283 
284  // Find which root corresponds to the desired value
285  scalar rootOld = SMALL;
286  scalar target = 0;
287 
288  label nRoots = 0;
289  Roots<3> selectedRoots;
290 
291  const scalar pMax = cmptMax(pCubicAlphas_);
292  const scalar pMin = cmptMin(pCubicAlphas_);
293  const scalar sFactor = pCubicAlphas_.d() - pCubicAlphas_.a();
294 
295  forAll(roots, rooti)
296  {
297  // Scale the roots back into original scale
298  const scalar root = roots[rooti]*sFactor + pCubicAlphas_.a();
299 
300  // Pick up correct root
301  if (root < pMax && root > pMin && rootOld != root)
302  {
303  target = (target == 0) ? root : target;
304  selectedRoots[nRoots++] = root;
305  }
306 
307  rootOld = root;
308  }
309 
310  // Recompute alpha last for analytical approach
311  cutAlpha_ = calcAlpha(cellInfo, target, tetDecom);
312 
313  // In case the selection of the root failed compute all three volume
314  // fractions and choose the one with minimum error
315  scalar error = mag(cutAlpha_ - cellAlpha);
316 
317  // If error > 1e-3 check for better root
318  if (nRoots > 0 && error > 1e-3)
319  {
320  scalar minError = error;
321  label minIndex = 0;
322 
323  for (label rooti=1; rooti<nRoots; rooti++)
324  {
325  const scalar targeti =
326  calcAlpha(cellInfo, selectedRoots[rooti], tetDecom);
327 
328  error = mag(cellAlpha - targeti);
329 
330  if (error < minError)
331  {
332  minError = error;
333  minIndex = rooti;
334  }
335  }
336 
337  cutAlpha_ = calcAlpha(cellInfo, selectedRoots[minIndex], tetDecom);
338  }
339 }
340 
341 
342 Foam::scalar Foam::MPLICcell::calcAlpha
343 (
344  const MPLICcellStorage& cellInfo,
345  const scalar target,
346  const bool tetDecom
347 )
348 {
349  if (!tetDecom)
350  {
351  return calcCutCellVolumeAlpha(cellInfo, target);
352  }
353  else
354  {
355  return calcTetCutCellVolumeAlpha(cellInfo, target);
356  }
357 }
358 
359 
360 void Foam::MPLICcell::calcSubCellVolume()
361 {
362  vector cEst = subFaceCentres_[0];
363  for(label i = 1; i < subFaceCentres_.size(); i++)
364  {
365  cEst += subFaceCentres_[i];
366  }
367  cEst /= subFaceCentres_.size();
368 
369  subCellVolume_ = 0;
370  forAll(subFaceAreas_, i)
371  {
372  subCellVolume_ += subFaceAreas_[i] & (subFaceCentres_[i] - cEst);
373  }
374  subCellVolume_ /= 3.0;
375 }
376 
377 
378 Foam::scalar Foam::MPLICcell::calcCutCellVolumeAlpha
379 (
380  const MPLICcellStorage& cellInfo,
381  const scalar target
382 )
383 {
384  const scalar V = cellInfo.V();
385 
386  // Case when the cell has to be cut
387  if (cellInfo.cellAlphaMax() > target && cellInfo.cellAlphaMin() < target)
388  {
389  // Cut cell single cut if multicut detected use multicut
390  const bool status = singleCutCell(cellInfo, target);
391  if (!status && multiCut_)
392  {
393  multiCutCell(cellInfo, target);
394  }
395 
396  // Compute normal
397  cutNormal_ = normalised(cutSf_);
398 
399  // Calculate volume
400  if (subFaceCentres_.size() != 0)
401  {
402  calcSubCellVolume();
403  }
404 
405  // Snap negative volume cell to zero
406  if (subCellVolume_ <= 0)
407  {
408  resetFaceFields(cellInfo.size());
409  subCellVolume_ = 0;
410  return 0;
411  }
412 
413  return min(subCellVolume_, V)/V;
414  }
415  else if (target <= cellInfo.cellAlphaMin())
416  {
417  if (unweighted_)
418  {
419  subFaceMagSf_ = cellInfo.magSf();
420  }
421  else
422  {
423  alphaPhiU_ = phiU_;
424  }
425  subCellVolume_ = V;
426 
427  return 1;
428  }
429  else
430  {
431  resetFaceFields(cellInfo.size());
432  subCellVolume_ = 0;
433 
434  return 0;
435  }
436 }
437 
438 
439 Foam::scalar Foam::MPLICcell::calcTetCutCellVolumeAlpha
440 (
441  const MPLICcellStorage& cellInfo,
442  const scalar target
443 )
444 {
445  clear();
446  resetFaceFields(cellInfo.size());
447 
448  // Append cell centre value
449  pointsAlpha_ = cellInfo.pointsAlpha();
450  pointsAlpha_.append(target);
451 
452  // Overall volume
453  scalar cellVolume = 0;
454 
455  if (min(pointsAlpha_) < target && max(pointsAlpha_) > target)
456  {
457  // Cell centre is the first point of the tet for all tets in the cell
458  const vector& a = cellInfo.C();
459 
460  // Cell centre value is the first value of the tet
461  // for all tets in the cell
462  tetPointsAlpha_[0] = cellInfo.cellAlpha();
463 
464  if (!unweighted_)
465  {
466  tetPointsU_[0] = cellInfo.cellU();
467  }
468 
469  // Looping through all the faces
470  forAll(cellInfo, facei)
471  {
472  // Create copy of the face indexing in order to flip if necessary
473  face f = cellInfo.faces()[cellInfo.cellFaces()[facei]];
474 
475  // Work directly with all the faces pointing out of the cell
476  if (!cellInfo.isOwner()[facei])
477  {
478  f.flip();
479  }
480 
481  const label& bL = f[0];
482  const point& b = cellInfo.points()[bL];
483  tetPointsAlpha_[1] = cellInfo.pointsAlpha()[bL];
484  if (!unweighted_)
485  {
486  tetPointsU_[1] = cellInfo.pointsU()[bL];
487  }
488 
489  // Decomposing faces
490  for (label i = 1; i < f.size()-1; ++i)
491  {
492  // Labels for point c and d
493  const label cL = f[i];
494  const label dL = f[i + 1];
495 
496  // c, d points of tetrahedron
497  const point& c = cellInfo.points()[cL];
498  const point& d = cellInfo.points()[dL];
499 
500  // Tet point values
501  tetPointsAlpha_[2] = cellInfo.pointsAlpha()[cL];
502  tetPointsAlpha_[3] = cellInfo.pointsAlpha()[dL];
503 
504  if (!unweighted_)
505  {
506  tetPointsU_[2] = cellInfo.pointsU()[cL];
507  tetPointsU_[3] = cellInfo.pointsU()[dL];
508  }
509 
510  // Tet maximum and minimum point values
511  const scalar tetMax = max(tetPointsAlpha_);
512  const scalar tetMin = min(tetPointsAlpha_);
513 
514  // Contains all the geometric information
515  tetPointRef cellTet(a, b, c, d);
516 
517  // Integrate overall volume for consistency
518  cellVolume += cellTet.mag();
519 
520  // Tet cuts
521  if (tetMin < target && tetMax > target)
522  {
523  // Tet point
524  tetPoints_ = {a, b, c, d};
525  tetSf_[0] = cellTet.Sa();
526  tetSf_[1] = cellTet.Sb();
527  tetSf_[2] = cellTet.Sc();
528  tetSf_[3] = cellTet.Sd();
529  tetCf_[0] = triPointRef(b, c, d).centre();
530  tetCf_[1] = triPointRef(a, d, c).centre();
531  tetCf_[2] = triPointRef(a, b, d).centre();
532  tetCf_[3] = triPointRef(a, c, b).centre();
533 
534  // Geometric cut of tetrahedra cell
535  const bool ow = cellInfo.isOwner()[facei];
536 
537  // Tetrahedron cut
538  cutTetCell(target, facei,ow);
539 
540  if (subFaceCentres_.size() > 0)
541  {
542  calcSubCellVolume();
543  }
544  }
545 
546  // Fully submerged tet
547  else if (tetMin >= target)
548  {
549  subCellVolume_ += cellTet.mag();
550  if (unweighted_)
551  {
552  subFaceMagSf_[facei] += mag(cellTet.Sa());
553  }
554  else
555  {
556  const scalar phiU =
557  (
558  (1.0/3.0)*
559  (
560  tetPointsU_[1]
561  + tetPointsU_[2]
562  + tetPointsU_[3]
563  )
564  ) & cellTet.Sa();
565 
566  if (cellInfo.isOwner()[facei])
567  {
568  alphaPhiU_[facei] += phiU;
569  }
570  else
571  {
572  alphaPhiU_[facei] -= phiU;
573  }
574  }
575  }
576  }
577  }
578 
579  // Compute normal
580  cutNormal_ = normalised(cutSf_);
581 
582  // Snap negative volume cell to zero
583  if (subCellVolume_ <= 0)
584  {
585  resetFaceFields(cellInfo.size());
586  subCellVolume_ = 0;
587  return 0;
588  }
589  return min(subCellVolume_, cellVolume)/cellVolume;
590  }
591  else if (target <= min(pointsAlpha_))
592  {
593  if (unweighted_)
594  {
595  subFaceMagSf_ = cellInfo.magSf();
596  }
597  else
598  {
599  alphaPhiU_ = phiU_;
600  }
601  subCellVolume_ = cellVolume;
602  return 1;
603  }
604  else
605  {
606  resetFaceFields(cellInfo.size());
607  subCellVolume_ = 0;
608  return 0;
609  }
610 }
611 
612 
613 bool Foam::MPLICcell::singleCutCell
614 (
615  const MPLICcellStorage& cellInfo,
616  const scalar target
617 )
618 {
619  clear();
620  resetFaceFields(cellInfo.size());
621 
622  // Cut type
623  label cutType;
624 
625  // Any face has more then one cut?
626  bool moreCutsPerFace = 0;
627 
628  // Single cell cut
629  forAll(cellInfo, facei)
630  {
631  // Collect fully submerged faces
632  if (cellInfo.facesAlphaMin()[facei] >= target)
633  {
634  appendSfCf
635  (
636  cellInfo.Sf()[facei],
637  cellInfo.Cf()[facei],
638  cellInfo.magSf()[facei],
639  cellInfo.isOwner()[facei]
640  );
641 
642  if (unweighted_)
643  {
644  subFaceMagSf_[facei] = cellInfo.magSf()[facei];
645  }
646  else
647  {
648  alphaPhiU_[facei] = phiU_[facei];
649  }
650  continue;
651  }
652  else if (cellInfo.facesAlphaMax()[facei] < target)
653  {
654  continue;
655  }
656 
657  // Cut the face return label of next face and edge
658  cutType = faceCutter_.cutFace
659  (
660  cellInfo.faces()[cellInfo.cellFaces()[facei]],
661  cellInfo.points(),
662  cellInfo.pointsAlpha(),
663  cellInfo.pointsU(),
664  target,
665  cellInfo.isOwner()[facei]
666  );
667 
668  // Potentially multiple cuts through the cell
669  if (cutType == -1)
670  {
671  moreCutsPerFace = 1;
672  }
673 
674  else if (cutType == 1)
675  {
676  // Append to the cut list of points
677  cutPoints_.append(faceCutter_.cutPoints());
678 
679  // Append area vectors and face centers
680  if (faceCutter_.subPoints().size() > 2)
681  {
682  const vector Sf = faceCutter_.Sf();
683  const vector Cf = faceCutter_.Cf(Sf);
684  const scalar magSf = mag(Sf);
685  appendSfCf(Sf, Cf, magSf);
686 
687  if (unweighted_)
688  {
689  subFaceMagSf_[facei] += magSf;
690  }
691  else
692  {
693  alphaPhiU_[facei] += faceCutter_.alphaPhiU();
694  }
695  }
696  }
697  }
698 
699  // Assume it is multicut if triangle have opposite sign in any direction
700  bool cutOrientationDiffers = 0;
701  if (cutPoints_.size() > 2)
702  {
703  cutOrientationDiffers = cutStatusCalcSf();
704  const vector Cf = calcCutCf(cutSf_);
705  appendSfCf(cutSf_, Cf, mag(cutSf_));
706  }
707 
708  // Potentially multiple cuts through cell
709  if (cutOrientationDiffers || moreCutsPerFace)
710  {
711  return 0;
712  }
713 
714  // Only one cut through the cell
715  else
716  {
717  return 1;
718  }
719 }
720 
721 
722 bool Foam::MPLICcell::multiCutCell
723 (
724  const MPLICcellStorage& cellInfo,
725  const scalar target
726 )
727 {
728  clear();
729  resetFaceFields(cellInfo.size());
730 
731  // Prepare local addressing
732  if (!addressingCalculated_)
733  {
734  calcAddressing(cellInfo);
735  }
736 
737  // Keep track of cut edges
738  boolList isEdgeCutOld(cellInfo.cellEdges().size(), false);
739  boolList isEdgeCut(cellInfo.cellEdges().size(), false);
740 
741  // Keep track of the fully submerged subfaces
742  boolList submerged(cellInfo.size(), false);
743 
744  // Initialise the list of necessary labels
745  label facei, nextFace, faceEdgei, status;
746 
747  // Loop through all the cuts
748  // Number of cuts limited to number of faces
749  forAll(cellInfo, cutI)
750  {
751  faceEdgei = -1;
752  facei = 0;
753  nextFace = 0;
754  status = 0;
755 
756  // One cell cut
757  label j = 0;
758  while (j < cellInfo.size())
759  {
760  facei = (status == 0) ? j : nextFace;
761 
762  // Collect fully submerged faces
763  if (cellInfo.facesAlphaMin()[facei] >= target && !submerged[facei])
764  {
765  submerged[facei] = true;
766  appendSfCf
767  (
768  cellInfo.Sf()[facei],
769  cellInfo.Cf()[facei],
770  cellInfo.magSf()[facei],
771  cellInfo.isOwner()[facei]
772  );
773 
774  // Precompute face fields
775  if (unweighted_)
776  {
777  subFaceMagSf_[facei] = cellInfo.magSf()[facei];
778  }
779  else
780  {
781  alphaPhiU_[facei] = phiU_[facei];
782  }
783  }
784 
785  // Cut the face
786  status = faceCutter_.cutFace
787  (
788  cellInfo.faces()[cellInfo.cellFaces()[facei]],
789  localFaceEdges_[facei],
790  cellInfo.points(),
791  isEdgeCutOld,
792  isEdgeCut,
793  faceEdgei,
794  cellInfo.pointsAlpha(),
795  cellInfo.pointsU(),
796  facei,
797  target,
798  cellInfo.isOwner()[facei]
799  );
800 
801  // Get the next face and edge
802  if (status)
803  {
804  const label edgei = localFaceEdges_[facei][faceEdgei];
805  const labelList& edgeFaces = localEdgeFaces_[edgei];
806  nextFace = edgeFaces[edgeFaces[0] == facei];
807  faceEdgei = findIndex(localFaceEdges_[nextFace], edgei);
808  }
809 
810  // Append to the cut list of points
811  cutPoints_.append(faceCutter_.cutPoints());
812  cutEdges_.append(faceCutter_.cutEdges());
813 
814  // Append area vectors and face centers
815  if (faceCutter_.subPoints().size() > 2 && !submerged[facei])
816  {
817  const vector Sf = faceCutter_.Sf();
818  const vector Cf = faceCutter_.Cf(Sf);
819  const scalar magSf = mag(Sf);
820 
821  // The sub-faces are always pointing outwards
822  appendSfCf(Sf, Cf, magSf);
823 
824  if (unweighted_)
825  {
826  subFaceMagSf_[facei] += magSf;
827  }
828  else
829  {
830  alphaPhiU_[facei] += faceCutter_.alphaPhiU();
831  }
832  }
833 
834  // End on reaching first edge
835  if (cutEdges_.size() > 0 && cutEdges_.first() == cutEdges_.last())
836  {
837  break;
838  }
839  if (status == 0)
840  {
841  ++j;
842  }
843  }
844 
845  isEdgeCutOld = isEdgeCut;
846 
847  if (cutPoints_.size() == 0)
848  {
849  break;
850  }
851  else
852  {
853  // Append information from cut face
854  const vector Sf = calcCutSf();
855  const vector Cf = calcCutCf(Sf);
856  appendSfCf(Sf, Cf, mag(Sf));
857  cutSf_ += Sf;
858  }
859 
860  // Clear fields to prepare for next cut
861  cutPoints_.clear();
862  cutEdges_.clear();
863  }
864  return 1;
865 }
866 
867 
868 bool Foam::MPLICcell::cutTetCell
869 (
870  const scalar target,
871  const label faceOrig,
872  const bool ow
873 )
874 {
875  // Clear geometry data for tet cut
876  clearOneCut();
877 
878  // Single cell cut
879  forAll(tetFaces_, facei)
880  {
881  const face& f = tetFaces_[facei];
882 
883  // Collect fully submerged faces
884  if
885  (
886  min
887  (
888  min
889  (
890  tetPointsAlpha_[f[0]],
891  tetPointsAlpha_[f[1]]
892  ),
893  tetPointsAlpha_[f[2]]
894  ) >= target
895  )
896  {
897  const vector& Sf = tetSf_[facei];
898  const vector& Cf = tetCf_[facei];
899  appendSfCf(Sf, Cf, mag(Sf));
900 
901  if (unweighted_ && facei == 0)
902  {
903  subFaceMagSf_[faceOrig] += mag(Sf);
904  }
905  else if (!unweighted_ && facei == 0)
906  {
907  const face& f0 = tetFaces_[0];
908  const scalar phiU =
909  (
910  (1.0/3.0)*
911  (
912  tetPointsU_[f0[0]]
913  + tetPointsU_[f0[1]]
914  + tetPointsU_[f0[2]]
915  )
916  ) & Sf;
917  if (ow)
918  {
919  alphaPhiU_[faceOrig] += phiU;
920  }
921  else
922  {
923  alphaPhiU_[faceOrig] -= phiU;
924  }
925  }
926  continue;
927  }
928  else if
929  (
930  max
931  (
932  max
933  (
934  tetPointsAlpha_[f[0]],
935  tetPointsAlpha_[f[1]]
936  ),
937  tetPointsAlpha_[f[2]]
938  ) < target
939  )
940  {
941  continue;
942  }
943 
944  // Cut the face return label of next face and edge
945  faceCutter_.cutFace
946  (
947  tetFaces_[facei],
948  tetPoints_,
949  tetPointsAlpha_,
950  tetPointsU_,
951  target,
952  true
953  );
954 
955  // Append to the cut list of points
956  cutPoints_.append(faceCutter_.cutPoints());
957 
958  // Append area vectors and face centers
959  if (faceCutter_.subPoints().size() > 2)
960  {
961  const vector Sf = faceCutter_.Sf();
962  const vector Cf = faceCutter_.Cf(Sf);
963  const scalar magSf = mag(Sf);
964  appendSfCf(Sf, Cf, magSf);
965 
966  // For unweighted alphaf
967  if (unweighted_ && facei == 0)
968  {
969  subFaceMagSf_[faceOrig] += magSf;
970  }
971 
972  // For phiU weighted alphaf
973  else if (!unweighted_ && facei == 0)
974  {
975  if (ow)
976  {
977  alphaPhiU_[faceOrig] += faceCutter_.alphaPhiU();
978  }
979  else
980  {
981  alphaPhiU_[faceOrig] -= faceCutter_.alphaPhiU();
982  }
983  }
984  }
985  }
986 
987  // Append information from cut face
988  if (cutPoints_.size() > 2)
989  {
990  const vector Sf = calcCutSf();
991  const vector Cf = calcCutCf(Sf);
992  appendSfCf(Sf, Cf, mag(Sf));
993  cutSf_ += Sf;
994  }
995 
996  return 1;
997 }
998 
999 
1000 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1001 
1002 Foam::MPLICcell::MPLICcell(const bool unweighted, const bool multiCut)
1003 :
1004  unweighted_(unweighted),
1005  multiCut_(multiCut),
1006  faceCutter_(unweighted),
1007  cutPoints_(10),
1008  cutEdges_(10),
1009  subFaceAreas_(10),
1010  subFaceCentres_(10),
1011  tetPointsAlpha_(4),
1012  tetPointsU_(4),
1013  tetFaces_
1014  {
1015  triFace(1, 2, 3),
1016  triFace(0, 3, 2),
1017  triFace(0, 1, 3),
1018  triFace(0, 2, 1)
1019  },
1020  pointsAlpha_(8)
1021 {}
1022 
1023 
1024 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1025 
1028  const MPLICcellStorage& cellInfo
1029 )
1030 {
1031  // Addressing for multicut needs to be recomputed for each cell
1032  addressingCalculated_ = false;
1033 
1034  // Try normal cell cut matching first
1035  label status = calcMatchAlphaCutCell(cellInfo);
1036 
1037  // If volume fraction error is bigger than 10% try tet decomposition cut
1038  if (status == 0 && multiCut_)
1039  {
1040  status = calcMatchAlphaCutCell(cellInfo, true);
1041  }
1042 
1043  if (status == 0 || status == -1)
1044  {
1045  return 0;
1046  }
1047  else
1048  {
1049  return 1;
1050  }
1051 }
1052 
1053 
1054 // ************************************************************************* //
const vector Sf() const
Return subface surface area vector.
Definition: MPLICfaceI.H:115
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
const vector Cf(const vector &area) const
Return subface centre.
Definition: MPLICfaceI.H:122
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const DynamicList< point > & cutPoints() const
Access to cut points.
Definition: MPLICfaceI.H:96
const DynamicList< label > & cutEdges() const
Access to cut edges.
Definition: MPLICfaceI.H:109
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Holds information (coordinate and normal) regarding nearest wall point.
static const scalar SMALL
Definition: scalar.H:115
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Provides local cell addressing for geometry and data for MPLIC class.
const dimensionedScalar c
Speed of light in a vacuum.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
scalar cutAlpha() const
Return volume fraction corresponding to the cut.
Definition: MPLICcellI.H:287
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
face triFace(3)
void sort(UList< T > &)
Definition: UList.C:115
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
List< label > labelList
A List of labels.
Definition: labelList.H:56
MPLICcell(const bool unweighted=true, const bool multiCut=true)
Construct for given interpolation and PLIC type.
Definition: MPLICcell.C:1002
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
vector point
Point is a vector.
Definition: point.H:41
label cutFace(const labelList &f, const labelList &faceEdges, const pointField &points, const boolList &isEdgeCutOld, boolList &isEdgeCut, label &faceEdgei, const UList< scalar > &pointsAlpha, const UList< vector > &pointsU, const label facei, const scalar target, const bool ow)
Function to cut for multi cut.
Definition: MPLICface.C:40
bool matchAlpha(const MPLICcellStorage &cellInfo)
Match cell volume ratios.
Definition: MPLICcell.C:1027
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
const DynamicList< point > & subPoints() const
Access to submerged face points.
Definition: MPLICfaceI.H:103
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
scalar alphaPhiU() const
Calculate and return alphaPhiU.
Definition: MPLICfaceI.H:30