primitiveMeshGeometry.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 "primitiveMeshGeometry.H"
27 #include "pyramidPointFaceRef.H"
28 #include "unitConversion.H"
29 #include "triPointRef.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(primitiveMeshGeometry, 0);
36 }
37 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 void Foam::primitiveMeshGeometry::updateFaceCentresAndAreas
42 (
43  const pointField& p,
44  const labelList& changedFaces
45 )
46 {
47  const faceList& fs = mesh_.faces();
48 
49  forAll(changedFaces, i)
50  {
51  label facei = changedFaces[i];
52 
53  const labelList& f = fs[facei];
54  label nPoints = f.size();
55 
56  // If the face is a triangle, do a direct calculation for efficiency
57  // and to avoid round-off error-related problems
58  if (nPoints == 3)
59  {
60  faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
61  faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
62  }
63  else
64  {
65  vector sumN = Zero;
66  scalar sumA = 0.0;
67  vector sumAc = Zero;
68 
69  point fCentre = p[f[0]];
70  for (label pi = 1; pi < nPoints; pi++)
71  {
72  fCentre += p[f[pi]];
73  }
74 
75  fCentre /= nPoints;
76 
77  for (label pi = 0; pi < nPoints; pi++)
78  {
79  const point& nextPoint = p[f[(pi + 1) % nPoints]];
80 
81  vector c = p[f[pi]] + nextPoint + fCentre;
82  vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
83  scalar a = mag(n);
84 
85  sumN += n;
86  sumA += a;
87  sumAc += a*c;
88  }
89 
90  faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
91  faceAreas_[facei] = 0.5*sumN;
92  }
93  }
94 }
95 
96 
97 void Foam::primitiveMeshGeometry::updateCellCentresAndVols
98 (
99  const labelList& changedCells,
100  const labelList& changedFaces
101 )
102 {
103  // Clear the fields for accumulation
104  UIndirectList<vector>(cellCentres_, changedCells) = Zero;
105  UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
106 
107  const labelList& own = mesh_.faceOwner();
108  const labelList& nei = mesh_.faceNeighbour();
109 
110  // first estimate the approximate cell centre as the average of face centres
111 
112  vectorField cEst(mesh_.nCells());
113  UIndirectList<vector>(cEst, changedCells) = Zero;
114  scalarField nCellFaces(mesh_.nCells());
115  UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
116 
117  forAll(changedFaces, i)
118  {
119  label facei = changedFaces[i];
120  cEst[own[facei]] += faceCentres_[facei];
121  nCellFaces[own[facei]] += 1;
122 
123  if (mesh_.isInternalFace(facei))
124  {
125  cEst[nei[facei]] += faceCentres_[facei];
126  nCellFaces[nei[facei]] += 1;
127  }
128  }
129 
130  forAll(changedCells, i)
131  {
132  label celli = changedCells[i];
133  cEst[celli] /= nCellFaces[celli];
134  }
135 
136  forAll(changedFaces, i)
137  {
138  label facei = changedFaces[i];
139 
140  // Calculate 3*face-pyramid volume
141  scalar pyr3Vol = max
142  (
143  faceAreas_[facei] & (faceCentres_[facei] - cEst[own[facei]]),
144  VSMALL
145  );
146 
147  // Calculate face-pyramid centre
148  vector pc = (3.0/4.0)*faceCentres_[facei] + (1.0/4.0)*cEst[own[facei]];
149 
150  // Accumulate volume-weighted face-pyramid centre
151  cellCentres_[own[facei]] += pyr3Vol*pc;
152 
153  // Accumulate face-pyramid volume
154  cellVolumes_[own[facei]] += pyr3Vol;
155 
156  if (mesh_.isInternalFace(facei))
157  {
158  // Calculate 3*face-pyramid volume
159  scalar pyr3Vol = max
160  (
161  faceAreas_[facei] & (cEst[nei[facei]] - faceCentres_[facei]),
162  VSMALL
163  );
164 
165  // Calculate face-pyramid centre
166  vector pc =
167  (3.0/4.0)*faceCentres_[facei]
168  + (1.0/4.0)*cEst[nei[facei]];
169 
170  // Accumulate volume-weighted face-pyramid centre
171  cellCentres_[nei[facei]] += pyr3Vol*pc;
172 
173  // Accumulate face-pyramid volume
174  cellVolumes_[nei[facei]] += pyr3Vol;
175  }
176  }
177 
178  forAll(changedCells, i)
179  {
180  label celli = changedCells[i];
181 
182  cellCentres_[celli] /= cellVolumes_[celli];
183  cellVolumes_[celli] *= (1.0/3.0);
184  }
185 }
186 
187 
189 (
190  const labelList& changedFaces
191 ) const
192 {
193  const labelList& own = mesh_.faceOwner();
194  const labelList& nei = mesh_.faceNeighbour();
195 
196  labelHashSet affectedCells(2*changedFaces.size());
197 
198  forAll(changedFaces, i)
199  {
200  label facei = changedFaces[i];
201 
202  affectedCells.insert(own[facei]);
203 
204  if (mesh_.isInternalFace(facei))
205  {
206  affectedCells.insert(nei[facei]);
207  }
208  }
209  return affectedCells.toc();
210 }
211 
212 
213 
214 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
215 
216 // Construct from components
218 (
219  const primitiveMesh& mesh
220 )
221 :
222  mesh_(mesh)
223 {
224  correct();
225 }
226 
227 
228 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
229 
230 
231 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
232 
234 {
235  faceAreas_ = mesh_.faceAreas();
236  faceCentres_ = mesh_.faceCentres();
237  cellCentres_ = mesh_.cellCentres();
238  cellVolumes_ = mesh_.cellVolumes();
239 }
240 
241 
243 (
244  const pointField& p,
245  const labelList& changedFaces
246 )
247 {
248  // Update face quantities
249  updateFaceCentresAndAreas(p, changedFaces);
250  // Update cell quantities from face quantities
251  updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
252 }
253 
254 
256 (
257  const bool report,
258  const scalar orthWarn,
259  const primitiveMesh& mesh,
260  const vectorField& cellCentres,
261  const vectorField& faceAreas,
262  const labelList& checkFaces,
263  labelHashSet* setPtr
264 )
265 {
266  // for all internal faces check theat the d dot S product is positive
267 
268  const labelList& own = mesh.faceOwner();
269  const labelList& nei = mesh.faceNeighbour();
270 
271  // Severe nonorthogonality threshold
272  const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
273 
274  scalar minDDotS = GREAT;
275 
276  scalar sumDDotS = 0;
277 
278  label severeNonOrth = 0;
279 
280  label errorNonOrth = 0;
281 
282  forAll(checkFaces, i)
283  {
284  label facei = checkFaces[i];
285 
286  if (mesh.isInternalFace(facei))
287  {
288  vector d = cellCentres[nei[facei]] - cellCentres[own[facei]];
289  const vector& s = faceAreas[facei];
290 
291  scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
292 
293  if (dDotS < severeNonorthogonalityThreshold)
294  {
295  if (dDotS > SMALL)
296  {
297  if (report)
298  {
299  // Severe non-orthogonality but mesh still OK
300  Pout<< "Severe non-orthogonality for face " << facei
301  << " between cells " << own[facei]
302  << " and " << nei[facei]
303  << ": Angle = " << radToDeg(::acos(dDotS))
304  << " deg." << endl;
305  }
306 
307  if (setPtr)
308  {
309  setPtr->insert(facei);
310  }
311 
312  severeNonOrth++;
313  }
314  else
315  {
316  // Non-orthogonality greater than 90 deg
317  if (report)
318  {
320  << "Severe non-orthogonality detected for face "
321  << facei
322  << " between cells " << own[facei] << " and "
323  << nei[facei]
324  << ": Angle = " << radToDeg(::acos(dDotS))
325  << " deg." << endl;
326  }
327 
328  errorNonOrth++;
329 
330  if (setPtr)
331  {
332  setPtr->insert(facei);
333  }
334  }
335  }
336 
337  if (dDotS < minDDotS)
338  {
339  minDDotS = dDotS;
340  }
341 
342  sumDDotS += dDotS;
343  }
344  }
345 
346  reduce(minDDotS, minOp<scalar>());
347  reduce(sumDDotS, sumOp<scalar>());
348  reduce(severeNonOrth, sumOp<label>());
349  reduce(errorNonOrth, sumOp<label>());
350 
351  label neiSize = nei.size();
352  reduce(neiSize, sumOp<label>());
353 
354  // Only report if there are some internal faces
355  if (neiSize > 0)
356  {
357  if (report && minDDotS < severeNonorthogonalityThreshold)
358  {
359  Info<< "Number of non-orthogonality errors: " << errorNonOrth
360  << ". Number of severely non-orthogonal faces: "
361  << severeNonOrth << "." << endl;
362  }
363  }
364 
365  if (report)
366  {
367  if (neiSize > 0)
368  {
369  Info<< "Mesh non-orthogonality Max: "
370  << radToDeg(::acos(minDDotS))
371  << " average: " << radToDeg(::acos(sumDDotS/neiSize))
372  << endl;
373  }
374  }
375 
376  if (errorNonOrth > 0)
377  {
378  if (report)
379  {
381  << "Error in non-orthogonality detected" << endl;
382  }
383 
384  return true;
385  }
386  else
387  {
388  if (report)
389  {
390  Info<< "Non-orthogonality check OK.\n" << endl;
391  }
392 
393  return false;
394  }
395 }
396 
397 
399 (
400  const bool report,
401  const scalar minPyrVol,
402  const primitiveMesh& mesh,
403  const vectorField& cellCentres,
404  const pointField& p,
405  const labelList& checkFaces,
406  labelHashSet* setPtr
407 )
408 {
409  // check whether face area vector points to the cell with higher label
410  const labelList& own = mesh.faceOwner();
411  const labelList& nei = mesh.faceNeighbour();
412 
413  const faceList& f = mesh.faces();
414 
415  label nErrorPyrs = 0;
416 
417  forAll(checkFaces, i)
418  {
419  label facei = checkFaces[i];
420 
421  // Create the owner pyramid - it will have negative volume
422  scalar pyrVol = pyramidPointFaceRef
423  (
424  f[facei],
425  cellCentres[own[facei]]
426  ).mag(p);
427 
428  if (pyrVol > -minPyrVol)
429  {
430  if (report)
431  {
432  Pout<< "bool primitiveMeshGeometry::checkFacePyramids("
433  << "const bool, const scalar, const pointField&"
434  << ", const labelList&, labelHashSet*): "
435  << "face " << facei << " points the wrong way. " << endl
436  << "Pyramid volume: " << -pyrVol
437  << " Face " << f[facei] << " area: " << f[facei].mag(p)
438  << " Owner cell: " << own[facei] << endl
439  << "Owner cell vertex labels: "
440  << mesh.cells()[own[facei]].labels(f)
441  << endl;
442  }
443 
444 
445  if (setPtr)
446  {
447  setPtr->insert(facei);
448  }
449 
450  nErrorPyrs++;
451  }
452 
453  if (mesh.isInternalFace(facei))
454  {
455  // Create the neighbour pyramid - it will have positive volume
456  scalar pyrVol =
457  pyramidPointFaceRef(f[facei], cellCentres[nei[facei]]).mag(p);
458 
459  if (pyrVol < minPyrVol)
460  {
461  if (report)
462  {
463  Pout<< "bool primitiveMeshGeometry::checkFacePyramids("
464  << "const bool, const scalar, const pointField&"
465  << ", const labelList&, labelHashSet*): "
466  << "face " << facei << " points the wrong way. " << endl
467  << "Pyramid volume: " << -pyrVol
468  << " Face " << f[facei] << " area: " << f[facei].mag(p)
469  << " Neighbour cell: " << nei[facei] << endl
470  << "Neighbour cell vertex labels: "
471  << mesh.cells()[nei[facei]].labels(f)
472  << endl;
473  }
474 
475  if (setPtr)
476  {
477  setPtr->insert(facei);
478  }
479 
480  nErrorPyrs++;
481  }
482  }
483  }
484 
485  reduce(nErrorPyrs, sumOp<label>());
486 
487  if (nErrorPyrs > 0)
488  {
489  if (report)
490  {
492  << "Error in face pyramids: faces pointing the wrong way!"
493  << endl;
494  }
495 
496  return true;
497  }
498  else
499  {
500  if (report)
501  {
502  Info<< "Face pyramids OK.\n" << endl;
503  }
504 
505  return false;
506  }
507 }
508 
509 
511 (
512  const bool report,
513  const scalar internalSkew,
514  const scalar boundarySkew,
515  const primitiveMesh& mesh,
516  const vectorField& cellCentres,
517  const vectorField& faceCentres,
518  const vectorField& faceAreas,
519  const labelList& checkFaces,
520  labelHashSet* setPtr
521 )
522 {
523  // Warn if the skew correction vector is more than skew times
524  // larger than the face area vector
525 
526  const labelList& own = mesh.faceOwner();
527  const labelList& nei = mesh.faceNeighbour();
528 
529  scalar maxSkew = 0;
530 
531  label nWarnSkew = 0;
532 
533  forAll(checkFaces, i)
534  {
535  label facei = checkFaces[i];
536 
537  if (mesh.isInternalFace(facei))
538  {
539  scalar dOwn = mag(faceCentres[facei] - cellCentres[own[facei]]);
540  scalar dNei = mag(faceCentres[facei] - cellCentres[nei[facei]]);
541 
542  point faceIntersection =
543  cellCentres[own[facei]]*dNei/(dOwn+dNei)
544  + cellCentres[nei[facei]]*dOwn/(dOwn+dNei);
545 
546  scalar skewness =
547  mag(faceCentres[facei] - faceIntersection)
548  / (
549  mag(cellCentres[nei[facei]]-cellCentres[own[facei]])
550  + VSMALL
551  );
552 
553  // Check if the skewness vector is greater than the PN vector.
554  // This does not cause trouble but is a good indication of a poor
555  // mesh.
556  if (skewness > internalSkew)
557  {
558  if (report)
559  {
560  Pout<< "Severe skewness for face " << facei
561  << " skewness = " << skewness << endl;
562  }
563 
564  if (setPtr)
565  {
566  setPtr->insert(facei);
567  }
568 
569  nWarnSkew++;
570  }
571 
572  if (skewness > maxSkew)
573  {
574  maxSkew = skewness;
575  }
576  }
577  else
578  {
579  // Boundary faces: consider them to have only skewness error.
580  // (i.e. treat as if mirror cell on other side)
581 
582  vector faceNormal = faceAreas[facei];
583  faceNormal /= mag(faceNormal) + VSMALL;
584 
585  vector dOwn = faceCentres[facei] - cellCentres[own[facei]];
586 
587  vector dWall = faceNormal*(faceNormal & dOwn);
588 
589  point faceIntersection = cellCentres[own[facei]] + dWall;
590 
591  scalar skewness =
592  mag(faceCentres[facei] - faceIntersection)
593  /(2*mag(dWall) + VSMALL);
594 
595  // Check if the skewness vector is greater than the PN vector.
596  // This does not cause trouble but is a good indication of a poor
597  // mesh.
598  if (skewness > boundarySkew)
599  {
600  if (report)
601  {
602  Pout<< "Severe skewness for boundary face " << facei
603  << " skewness = " << skewness << endl;
604  }
605 
606  if (setPtr)
607  {
608  setPtr->insert(facei);
609  }
610 
611  nWarnSkew++;
612  }
613 
614  if (skewness > maxSkew)
615  {
616  maxSkew = skewness;
617  }
618  }
619  }
620 
621  reduce(maxSkew, maxOp<scalar>());
622  reduce(nWarnSkew, sumOp<label>());
623 
624  if (nWarnSkew > 0)
625  {
626  if (report)
627  {
629  << 100*maxSkew
630  << " percent.\nThis may impair the quality of the result." << nl
631  << nWarnSkew << " highly skew faces detected."
632  << endl;
633  }
634 
635  return true;
636  }
637  else
638  {
639  if (report)
640  {
641  Info<< "Max skewness = " << 100*maxSkew
642  << " percent. Face skewness OK.\n" << endl;
643  }
644 
645  return false;
646  }
647 }
648 
649 
651 (
652  const bool report,
653  const scalar warnWeight,
654  const primitiveMesh& mesh,
655  const vectorField& cellCentres,
656  const vectorField& faceCentres,
657  const vectorField& faceAreas,
658  const labelList& checkFaces,
659  labelHashSet* setPtr
660 )
661 {
662  // Warn if the delta factor (0..1) is too large.
663 
664  const labelList& own = mesh.faceOwner();
665  const labelList& nei = mesh.faceNeighbour();
666 
667  scalar minWeight = GREAT;
668 
669  label nWarnWeight = 0;
670 
671  forAll(checkFaces, i)
672  {
673  label facei = checkFaces[i];
674 
675  if (mesh.isInternalFace(facei))
676  {
677  const point& fc = faceCentres[facei];
678 
679  scalar dOwn = mag(faceAreas[facei] & (fc-cellCentres[own[facei]]));
680  scalar dNei = mag(faceAreas[facei] & (cellCentres[nei[facei]]-fc));
681 
682  scalar weight = min(dNei,dOwn)/(dNei+dOwn);
683 
684  if (weight < warnWeight)
685  {
686  if (report)
687  {
688  Pout<< "Small weighting factor for face " << facei
689  << " weight = " << weight << endl;
690  }
691 
692  if (setPtr)
693  {
694  setPtr->insert(facei);
695  }
696 
697  nWarnWeight++;
698  }
699 
700  minWeight = min(minWeight, weight);
701  }
702  }
703 
704  reduce(minWeight, minOp<scalar>());
705  reduce(nWarnWeight, sumOp<label>());
706 
707  if (minWeight < warnWeight)
708  {
709  if (report)
710  {
712  << minWeight << '.' << nl
713  << nWarnWeight << " faces with small weights detected."
714  << endl;
715  }
716 
717  return true;
718  }
719  else
720  {
721  if (report)
722  {
723  Info<< "Min weight = " << minWeight
724  << " percent. Weights OK.\n" << endl;
725  }
726 
727  return false;
728  }
729 }
730 
731 
732 // Check convexity of angles in a face. Allow a slight non-convexity.
733 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
734 // (if truly concave and points not visible from face centre the face-pyramid
735 // check in checkMesh will fail)
737 (
738  const bool report,
739  const scalar maxDeg,
740  const primitiveMesh& mesh,
741  const vectorField& faceAreas,
742  const pointField& p,
743  const labelList& checkFaces,
744  labelHashSet* setPtr
745 )
746 {
747  if (maxDeg < -SMALL || maxDeg > 180+SMALL)
748  {
750  << "maxDeg should be [0..180] but is now " << maxDeg
751  << abort(FatalError);
752  }
753 
754  const scalar maxSin = Foam::sin(degToRad(maxDeg));
755 
756  const faceList& fcs = mesh.faces();
757 
758  scalar maxEdgeSin = 0.0;
759 
760  label nConcave = 0;
761 
762  label errorFacei = -1;
763 
764  forAll(checkFaces, i)
765  {
766  label facei = checkFaces[i];
767 
768  const face& f = fcs[facei];
769 
770  vector faceNormal = faceAreas[facei];
771  faceNormal /= mag(faceNormal) + VSMALL;
772 
773  // Get edge from f[0] to f[size-1];
774  vector ePrev(p[f.first()] - p[f.last()]);
775  scalar magEPrev = mag(ePrev);
776  ePrev /= magEPrev + VSMALL;
777 
778  forAll(f, fp0)
779  {
780  // Get vertex after fp
781  label fp1 = f.fcIndex(fp0);
782 
783  // Normalized vector between two consecutive points
784  vector e10(p[f[fp1]] - p[f[fp0]]);
785  scalar magE10 = mag(e10);
786  e10 /= magE10 + VSMALL;
787 
788  if (magEPrev > SMALL && magE10 > SMALL)
789  {
790  vector edgeNormal = ePrev ^ e10;
791  scalar magEdgeNormal = mag(edgeNormal);
792 
793  if (magEdgeNormal < maxSin)
794  {
795  // Edges (almost) aligned -> face is ok.
796  }
797  else
798  {
799  // Check normal
800  edgeNormal /= magEdgeNormal;
801 
802  if ((edgeNormal & faceNormal) < SMALL)
803  {
804  if (facei != errorFacei)
805  {
806  // Count only one error per face.
807  errorFacei = facei;
808  nConcave++;
809  }
810 
811  if (setPtr)
812  {
813  setPtr->insert(facei);
814  }
815 
816  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
817  }
818  }
819  }
820 
821  ePrev = e10;
822  magEPrev = magE10;
823  }
824  }
825 
826  reduce(nConcave, sumOp<label>());
827  reduce(maxEdgeSin, maxOp<scalar>());
828 
829  if (report)
830  {
831  if (maxEdgeSin > SMALL)
832  {
833  scalar maxConcaveDegr =
834  radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
835 
836  Info<< "There are " << nConcave
837  << " faces with concave angles between consecutive"
838  << " edges. Max concave angle = " << maxConcaveDegr
839  << " degrees.\n" << endl;
840  }
841  else
842  {
843  Info<< "All angles in faces are convex or less than " << maxDeg
844  << " degrees concave.\n" << endl;
845  }
846  }
847 
848  if (nConcave > 0)
849  {
850  if (report)
851  {
853  << nConcave << " face points with severe concave angle (> "
854  << maxDeg << " deg) found.\n"
855  << endl;
856  }
857 
858  return true;
859  }
860  else
861  {
862  return false;
863  }
864 }
865 
866 
870 //bool Foam::primitiveMeshGeometry::checkFaceFlatness
871 //(
872 // const bool report,
873 // const scalar warnFlatness,
874 // const primitiveMesh& mesh,
875 // const vectorField& faceAreas,
876 // const vectorField& faceCentres,
877 // const pointField& p,
878 // const labelList& checkFaces,
879 // labelHashSet* setPtr
880 //)
881 //{
882 // if (warnFlatness < 0 || warnFlatness > 1)
883 // {
884 // FatalErrorInFunction
885 // << "warnFlatness should be [0..1] but is now " << warnFlatness
886 // << abort(FatalError);
887 // }
888 //
889 //
890 // const faceList& fcs = mesh.faces();
891 //
892 // // Areas are calculated as the sum of areas. (see
893 // // primitiveMeshFaceCentresAndAreas.C)
894 //
895 // label nWarped = 0;
896 //
897 // scalar minFlatness = GREAT;
898 // scalar sumFlatness = 0;
899 // label nSummed = 0;
900 //
901 // forAll(checkFaces, i)
902 // {
903 // label facei = checkFaces[i];
904 //
905 // const face& f = fcs[facei];
906 //
907 // scalar magArea = mag(faceAreas[facei]);
908 //
909 // if (f.size() > 3 && magArea > VSMALL)
910 // {
911 // const point& fc = faceCentres[facei];
912 //
913 // // Calculate the sum of magnitude of areas and compare to
914 // // magnitude of sum of areas.
915 //
916 // scalar sumA = 0.0;
917 //
918 // forAll(f, fp)
919 // {
920 // const point& thisPoint = p[f[fp]];
921 // const point& nextPoint = p[f.nextLabel(fp)];
922 //
923 // // Triangle around fc.
924 // vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
925 // sumA += mag(n);
926 // }
927 //
928 // scalar flatness = magArea / (sumA+VSMALL);
929 //
930 // sumFlatness += flatness;
931 // nSummed++;
932 //
933 // minFlatness = min(minFlatness, flatness);
934 //
935 // if (flatness < warnFlatness)
936 // {
937 // nWarped++;
938 //
939 // if (setPtr)
940 // {
941 // setPtr->insert(facei);
942 // }
943 // }
944 // }
945 // }
946 //
947 //
948 // reduce(nWarped, sumOp<label>());
949 // reduce(minFlatness, minOp<scalar>());
950 //
951 // reduce(nSummed, sumOp<label>());
952 // reduce(sumFlatness, sumOp<scalar>());
953 //
954 // if (report)
955 // {
956 // if (nSummed > 0)
957 // {
958 // Info<< "Face flatness (1 = flat, 0 = butterfly) : average = "
959 // << sumFlatness / nSummed << " min = " << minFlatness << endl;
960 // }
961 //
962 // if (nWarped> 0)
963 // {
964 // Info<< "There are " << nWarped
965 // << " faces with ratio between projected and actual area < "
966 // << warnFlatness
967 // << ".\nMinimum ratio (minimum flatness, maximum warpage) = "
968 // << minFlatness << nl << endl;
969 // }
970 // else
971 // {
972 // Info<< "All faces are flat in that the ratio between projected"
973 // << " and actual area is > " << warnFlatness << nl << endl;
974 // }
975 // }
976 //
977 // if (nWarped > 0)
978 // {
979 // if (report)
980 // {
981 // WarningInFunction
982 // << nWarped << " faces with severe warpage (flatness < "
983 // << warnFlatness << ") found.\n"
984 // << endl;
985 // }
986 //
987 // return true;
988 // }
989 // else
990 // {
991 // return false;
992 // }
993 //}
994 
995 
996 // Check twist of faces. Is calculated as the difference between areas of
997 // individual triangles and the overall area of the face (which ifself is
998 // is the average of the areas of the individual triangles).
1001  const bool report,
1002  const scalar minTwist,
1003  const primitiveMesh& mesh,
1004  const vectorField& faceAreas,
1005  const vectorField& faceCentres,
1006  const pointField& p,
1007  const labelList& checkFaces,
1008  labelHashSet* setPtr
1009 )
1010 {
1011  if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1012  {
1014  << "minTwist should be [-1..1] but is now " << minTwist
1015  << abort(FatalError);
1016  }
1017 
1018 
1019  const faceList& fcs = mesh.faces();
1020 
1021  // Areas are calculated as the sum of areas. (see
1022  // primitiveMeshFaceCentresAndAreas.C)
1023 
1024  label nWarped = 0;
1025 
1026  forAll(checkFaces, i)
1027  {
1028  label facei = checkFaces[i];
1029 
1030  const face& f = fcs[facei];
1031 
1032  scalar magArea = mag(faceAreas[facei]);
1033 
1034  if (f.size() > 3 && magArea > VSMALL)
1035  {
1036  const vector nf = faceAreas[facei] / magArea;
1037 
1038  const point& fc = faceCentres[facei];
1039 
1040  forAll(f, fpI)
1041  {
1042  vector triArea
1043  (
1044  triPointRef
1045  (
1046  p[f[fpI]],
1047  p[f.nextLabel(fpI)],
1048  fc
1049  ).normal()
1050  );
1051 
1052  scalar magTri = mag(triArea);
1053 
1054  if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1055  {
1056  nWarped++;
1057 
1058  if (setPtr)
1059  {
1060  setPtr->insert(facei);
1061  }
1062  }
1063  }
1064  }
1065  }
1066 
1067 
1068  reduce(nWarped, sumOp<label>());
1069 
1070  if (report)
1071  {
1072  if (nWarped> 0)
1073  {
1074  Info<< "There are " << nWarped
1075  << " faces with cosine of the angle"
1076  << " between triangle normal and face normal less than "
1077  << minTwist << nl << endl;
1078  }
1079  else
1080  {
1081  Info<< "All faces are flat in that the cosine of the angle"
1082  << " between triangle normal and face normal less than "
1083  << minTwist << nl << endl;
1084  }
1085  }
1086 
1087  if (nWarped > 0)
1088  {
1089  if (report)
1090  {
1092  << nWarped << " faces with severe warpage "
1093  << "(cosine of the angle between triangle normal and "
1094  << "face normal < " << minTwist << ") found.\n"
1095  << endl;
1096  }
1097 
1098  return true;
1099  }
1100  else
1101  {
1102  return false;
1103  }
1104 }
1105 
1106 
1109  const bool report,
1110  const scalar minArea,
1111  const primitiveMesh& mesh,
1112  const vectorField& faceAreas,
1113  const labelList& checkFaces,
1114  labelHashSet* setPtr
1115 )
1116 {
1117  label nZeroArea = 0;
1118 
1119  forAll(checkFaces, i)
1120  {
1121  label facei = checkFaces[i];
1122 
1123  if (mag(faceAreas[facei]) < minArea)
1124  {
1125  if (setPtr)
1126  {
1127  setPtr->insert(facei);
1128  }
1129  nZeroArea++;
1130  }
1131  }
1132 
1133 
1134  reduce(nZeroArea, sumOp<label>());
1135 
1136  if (report)
1137  {
1138  if (nZeroArea > 0)
1139  {
1140  Info<< "There are " << nZeroArea
1141  << " faces with area < " << minArea << '.' << nl << endl;
1142  }
1143  else
1144  {
1145  Info<< "All faces have area > " << minArea << '.' << nl << endl;
1146  }
1147  }
1148 
1149  if (nZeroArea > 0)
1150  {
1151  if (report)
1152  {
1154  << nZeroArea << " faces with area < " << minArea
1155  << " found.\n"
1156  << endl;
1157  }
1158 
1159  return true;
1160  }
1161  else
1162  {
1163  return false;
1164  }
1165 }
1166 
1167 
1170  const bool report,
1171  const scalar warnDet,
1172  const primitiveMesh& mesh,
1173  const vectorField& faceAreas,
1174  const labelList& checkFaces,
1175  const labelList& affectedCells,
1176  labelHashSet* setPtr
1177 )
1178 {
1179  const cellList& cells = mesh.cells();
1180 
1181  scalar minDet = GREAT;
1182  scalar sumDet = 0.0;
1183  label nSumDet = 0;
1184  label nWarnDet = 0;
1185 
1186  forAll(affectedCells, i)
1187  {
1188  const cell& cFaces = cells[affectedCells[i]];
1189 
1190  tensor areaSum(Zero);
1191  scalar magAreaSum = 0;
1192 
1193  forAll(cFaces, cFacei)
1194  {
1195  label facei = cFaces[cFacei];
1196 
1197  scalar magArea = mag(faceAreas[facei]);
1198 
1199  magAreaSum += magArea;
1200  areaSum += faceAreas[facei]*(faceAreas[facei]/magArea);
1201  }
1202 
1203  scalar scaledDet = det(areaSum/magAreaSum)/0.037037037037037;
1204 
1205  minDet = min(minDet, scaledDet);
1206  sumDet += scaledDet;
1207  nSumDet++;
1208 
1209  if (scaledDet < warnDet)
1210  {
1211  if (setPtr)
1212  {
1213  // Insert all faces of the cell.
1214  forAll(cFaces, cFacei)
1215  {
1216  label facei = cFaces[cFacei];
1217  setPtr->insert(facei);
1218  }
1219  }
1220  nWarnDet++;
1221  }
1222  }
1223 
1224  reduce(minDet, minOp<scalar>());
1225  reduce(sumDet, sumOp<scalar>());
1226  reduce(nSumDet, sumOp<label>());
1227  reduce(nWarnDet, sumOp<label>());
1228 
1229  if (report)
1230  {
1231  if (nSumDet > 0)
1232  {
1233  Info<< "Cell determinant (1 = uniform cube) : average = "
1234  << sumDet / nSumDet << " min = " << minDet << endl;
1235  }
1236 
1237  if (nWarnDet > 0)
1238  {
1239  Info<< "There are " << nWarnDet
1240  << " cells with determinant < " << warnDet << '.' << nl
1241  << endl;
1242  }
1243  else
1244  {
1245  Info<< "All faces have determinant > " << warnDet << '.' << nl
1246  << endl;
1247  }
1248  }
1249 
1250  if (nWarnDet > 0)
1251  {
1252  if (report)
1253  {
1255  << nWarnDet << " cells with determinant < " << warnDet
1256  << " found.\n"
1257  << endl;
1258  }
1259 
1260  return true;
1261  }
1262  else
1263  {
1264  return false;
1265  }
1266 }
1267 
1268 
1271  const bool report,
1272  const scalar orthWarn,
1273  const labelList& checkFaces,
1274  labelHashSet* setPtr
1275 ) const
1276 {
1277  return checkFaceDotProduct
1278  (
1279  report,
1280  orthWarn,
1281  mesh_,
1282  cellCentres_,
1283  faceAreas_,
1284  checkFaces,
1285  setPtr
1286  );
1287 }
1288 
1289 
1292  const bool report,
1293  const scalar minPyrVol,
1294  const pointField& p,
1295  const labelList& checkFaces,
1296  labelHashSet* setPtr
1297 ) const
1298 {
1299  return checkFacePyramids
1300  (
1301  report,
1302  minPyrVol,
1303  mesh_,
1304  cellCentres_,
1305  p,
1306  checkFaces,
1307  setPtr
1308  );
1309 }
1310 
1311 
1314  const bool report,
1315  const scalar internalSkew,
1316  const scalar boundarySkew,
1317  const labelList& checkFaces,
1318  labelHashSet* setPtr
1319 ) const
1320 {
1321  return checkFaceSkewness
1322  (
1323  report,
1324  internalSkew,
1325  boundarySkew,
1326  mesh_,
1327  cellCentres_,
1328  faceCentres_,
1329  faceAreas_,
1330  checkFaces,
1331  setPtr
1332  );
1333 }
1334 
1335 
1338  const bool report,
1339  const scalar warnWeight,
1340  const labelList& checkFaces,
1341  labelHashSet* setPtr
1342 ) const
1343 {
1344  return checkFaceWeights
1345  (
1346  report,
1347  warnWeight,
1348  mesh_,
1349  cellCentres_,
1350  faceCentres_,
1351  faceAreas_,
1352  checkFaces,
1353  setPtr
1354  );
1355 }
1356 
1357 
1360  const bool report,
1361  const scalar maxDeg,
1362  const pointField& p,
1363  const labelList& checkFaces,
1364  labelHashSet* setPtr
1365 ) const
1366 {
1367  return checkFaceAngles
1368  (
1369  report,
1370  maxDeg,
1371  mesh_,
1372  faceAreas_,
1373  p,
1374  checkFaces,
1375  setPtr
1376  );
1377 }
1378 
1379 
1380 //bool Foam::primitiveMeshGeometry::checkFaceFlatness
1381 //(
1382 // const bool report,
1383 // const scalar warnFlatness,
1384 // const pointField& p,
1385 // const labelList& checkFaces,
1386 // labelHashSet* setPtr
1387 //) const
1388 //{
1389 // return checkFaceFlatness
1390 // (
1391 // report,
1392 // warnFlatness,
1393 // mesh_,
1394 // faceAreas_,
1395 // faceCentres_,
1396 // p,
1397 // checkFaces,
1398 // setPtr
1399 // );
1400 //}
1401 
1402 
1405  const bool report,
1406  const scalar minTwist,
1407  const pointField& p,
1408  const labelList& checkFaces,
1409  labelHashSet* setPtr
1410 ) const
1411 {
1412  return checkFaceTwist
1413  (
1414  report,
1415  minTwist,
1416  mesh_,
1417  faceAreas_,
1418  faceCentres_,
1419  p,
1420  checkFaces,
1421  setPtr
1422  );
1423 }
1424 
1425 
1428  const bool report,
1429  const scalar minArea,
1430  const labelList& checkFaces,
1431  labelHashSet* setPtr
1432 ) const
1433 {
1434  return checkFaceArea
1435  (
1436  report,
1437  minArea,
1438  mesh_,
1439  faceAreas_,
1440  checkFaces,
1441  setPtr
1442  );
1443 }
1444 
1445 
1448  const bool report,
1449  const scalar warnDet,
1450  const labelList& checkFaces,
1451  const labelList& affectedCells,
1452  labelHashSet* setPtr
1453 ) const
1454 {
1455  return checkCellDeterminant
1456  (
1457  report,
1458  warnDet,
1459  mesh_,
1460  faceAreas_,
1461  checkFaces,
1462  affectedCells,
1463  setPtr
1464  );
1465 }
1466 
1467 
1468 // ************************************************************************* //
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:58
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
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
const cellList & cells() const
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensionedScalar det(const dimensionedSphericalTensor &dt)
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
T & first()
Return the first element of the list.
Definition: UListI.H:114
static bool checkFacePyramids(const bool report, const scalar minPyrVol, const primitiveMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, labelHashSet *)
dimensionedScalar asin(const dimensionedScalar &ds)
labelList affectedCells(const labelList &changedFaces) const
Helper function: get affected cells from faces.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
static bool checkFaceArea(const bool report, const scalar minArea, const primitiveMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
virtual const labelList & faceNeighbour() const =0
Face face-neighbour addressing.
static bool checkCellDeterminant(const bool report, const scalar minDet, const primitiveMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
static bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const primitiveMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
label nPoints
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static bool checkFaceAngles(const bool report, const scalar maxDeg, const primitiveMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
primitiveMeshGeometry(const primitiveMesh &)
Construct from mesh.
pyramid< point, const point &, const face & > pyramidPointFaceRef
errorManip< error > abort(error &err)
Definition: errorManip.H:131
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
void correct()
Take over properties from mesh.
dimensionedScalar sin(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:262
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/ubuntu/OpenFOAM-5.0/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:72
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
A normal distribution model.
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const primitiveMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
static bool checkFaceWeights(const bool report, const scalar warnWeight, const primitiveMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
virtual const faceList & faces() const =0
Return faces.
const dimensionedScalar c
Speed of light in a vacuum.
messageStream Info
virtual const labelList & faceOwner() const =0
Face face-owner addresing.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
T & last()
Return the last element of the list.
Definition: UListI.H:128
static bool checkFaceTwist(const bool report, const scalar minTwist, const primitiveMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Namespace for OpenFOAM.