FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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"
28 
29 namespace Foam
30 {
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 defineTypeNameAndDebug(primitiveMeshGeometry, 0);
35 
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 = vector::zero;
66  scalar sumA = 0.0;
67  vector sumAc = vector::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) = vector::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) = vector::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 
233 //- Take over properties from mesh
235 {
236  faceAreas_ = mesh_.faceAreas();
237  faceCentres_ = mesh_.faceCentres();
238  cellCentres_ = mesh_.cellCentres();
239  cellVolumes_ = mesh_.cellVolumes();
240 }
241 
242 
243 //- Recalculate on selected faces
245 (
246  const pointField& p,
247  const labelList& changedFaces
248 )
249 {
250  // Update face quantities
251  updateFaceCentresAndAreas(p, changedFaces);
252  // Update cell quantities from face quantities
253  updateCellCentresAndVols(affectedCells(changedFaces), changedFaces);
254 }
255 
256 
258 (
259  const bool report,
260  const scalar orthWarn,
261  const primitiveMesh& mesh,
262  const vectorField& cellCentres,
263  const vectorField& faceAreas,
264  const labelList& checkFaces,
265  labelHashSet* setPtr
266 )
267 {
268  // for all internal faces check theat the d dot S product is positive
269 
270  const labelList& own = mesh.faceOwner();
271  const labelList& nei = mesh.faceNeighbour();
272 
273  // Severe nonorthogonality threshold
274  const scalar severeNonorthogonalityThreshold =
275  ::cos(orthWarn/180.0*mathematicalConstant::pi);
276 
277  scalar minDDotS = GREAT;
278 
279  scalar sumDDotS = 0;
280 
281  label severeNonOrth = 0;
282 
283  label errorNonOrth = 0;
284 
285  forAll(checkFaces, i)
286  {
287  label faceI = checkFaces[i];
288 
289  if (mesh.isInternalFace(faceI))
290  {
291  vector d = cellCentres[nei[faceI]] - cellCentres[own[faceI]];
292  const vector& s = faceAreas[faceI];
293 
294  scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
295 
296  if (dDotS < severeNonorthogonalityThreshold)
297  {
298  if (dDotS > SMALL)
299  {
300  if (report)
301  {
302  // Severe non-orthogonality but mesh still OK
303  Pout<< "Severe non-orthogonality for face " << faceI
304  << " between cells " << own[faceI]
305  << " and " << nei[faceI]
306  << ": Angle = "
307  << ::acos(dDotS)/mathematicalConstant::pi*180.0
308  << " deg." << endl;
309  }
310 
311  if (setPtr)
312  {
313  setPtr->insert(faceI);
314  }
315 
316  severeNonOrth++;
317  }
318  else
319  {
320  // Non-orthogonality greater than 90 deg
321  if (report)
322  {
323  WarningIn
324  (
325  "primitiveMeshGeometry::checkFaceDotProduct"
326  "(const bool, const scalar, const labelList&"
327  ", labelHashSet*)"
328  ) << "Severe non-orthogonality detected for face "
329  << faceI
330  << " between cells " << own[faceI] << " and "
331  << nei[faceI]
332  << ": Angle = "
333  << ::acos(dDotS)/mathematicalConstant::pi*180.0
334  << " deg." << endl;
335  }
336 
337  errorNonOrth++;
338 
339  if (setPtr)
340  {
341  setPtr->insert(faceI);
342  }
343  }
344  }
345 
346  if (dDotS < minDDotS)
347  {
348  minDDotS = dDotS;
349  }
350 
351  sumDDotS += dDotS;
352  }
353  }
354 
355  reduce(minDDotS, minOp<scalar>());
356  reduce(sumDDotS, sumOp<scalar>());
357  reduce(severeNonOrth, sumOp<label>());
358  reduce(errorNonOrth, sumOp<label>());
359 
360  label neiSize = nei.size();
361  reduce(neiSize, sumOp<label>());
362 
363  // Only report if there are some internal faces
364  if (neiSize > 0)
365  {
366  if (report && minDDotS < severeNonorthogonalityThreshold)
367  {
368  Info<< "Number of non-orthogonality errors: " << errorNonOrth
369  << ". Number of severely non-orthogonal faces: "
370  << severeNonOrth << "." << endl;
371  }
372  }
373 
374  if (report)
375  {
376  if (neiSize > 0)
377  {
378  Info<< "Mesh non-orthogonality Max: "
379  << ::acos(minDDotS)/mathematicalConstant::pi*180.0
380  << " average: " <<
381  ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0
382  << endl;
383  }
384  }
385 
386  if (errorNonOrth > 0)
387  {
388  if (report)
389  {
391  (
392  "primitiveMeshGeometry::checkFaceDotProduct"
393  "(const bool, const scalar, const labelList&, labelHashSet*)"
394  ) << "Error in non-orthogonality detected" << endl;
395  }
396 
397  return true;
398  }
399  else
400  {
401  if (report)
402  {
403  Info<< "Non-orthogonality check OK.\n" << endl;
404  }
405 
406  return false;
407  }
408 }
409 
410 
412 (
413  const bool report,
414  const scalar minPyrVol,
415  const primitiveMesh& mesh,
416  const vectorField& cellCentres,
417  const pointField& p,
418  const labelList& checkFaces,
419  labelHashSet* setPtr
420 )
421 {
422  // check whether face area vector points to the cell with higher label
423  const labelList& own = mesh.faceOwner();
424  const labelList& nei = mesh.faceNeighbour();
425 
426  const faceList& f = mesh.faces();
427 
428  label nErrorPyrs = 0;
429 
430  forAll(checkFaces, i)
431  {
432  label faceI = checkFaces[i];
433 
434  // Create the owner pyramid - it will have negative volume
435  scalar pyrVol = pyramidPointFaceRef
436  (
437  f[faceI],
438  cellCentres[own[faceI]]
439  ).mag(p);
440 
441  if (pyrVol > -minPyrVol)
442  {
443  if (report)
444  {
445  Pout<< "bool primitiveMeshGeometry::checkFacePyramids("
446  << "const bool, const scalar, const pointField&"
447  << ", const labelList&, labelHashSet*): "
448  << "face " << faceI << " points the wrong way. " << endl
449  << "Pyramid volume: " << -pyrVol
450  << " Face " << f[faceI] << " area: " << f[faceI].mag(p)
451  << " Owner cell: " << own[faceI] << endl
452  << "Owner cell vertex labels: "
453  << mesh.cells()[own[faceI]].labels(f)
454  << endl;
455  }
456 
457 
458  if (setPtr)
459  {
460  setPtr->insert(faceI);
461  }
462 
463  nErrorPyrs++;
464  }
465 
466  if (mesh.isInternalFace(faceI))
467  {
468  // Create the neighbour pyramid - it will have positive volume
469  scalar pyrVol =
470  pyramidPointFaceRef(f[faceI], cellCentres[nei[faceI]]).mag(p);
471 
472  if (pyrVol < minPyrVol)
473  {
474  if (report)
475  {
476  Pout<< "bool primitiveMeshGeometry::checkFacePyramids("
477  << "const bool, const scalar, const pointField&"
478  << ", const labelList&, labelHashSet*): "
479  << "face " << faceI << " points the wrong way. " << endl
480  << "Pyramid volume: " << -pyrVol
481  << " Face " << f[faceI] << " area: " << f[faceI].mag(p)
482  << " Neighbour cell: " << nei[faceI] << endl
483  << "Neighbour cell vertex labels: "
484  << mesh.cells()[nei[faceI]].labels(f)
485  << endl;
486  }
487 
488  if (setPtr)
489  {
490  setPtr->insert(faceI);
491  }
492 
493  nErrorPyrs++;
494  }
495  }
496  }
497 
498  reduce(nErrorPyrs, sumOp<label>());
499 
500  if (nErrorPyrs > 0)
501  {
502  if (report)
503  {
505  (
506  "primitiveMeshGeometry::checkFacePyramids("
507  "const bool, const scalar, const pointField&"
508  ", const labelList&, labelHashSet*)"
509  ) << "Error in face pyramids: faces pointing the wrong way!"
510  << endl;
511  }
512 
513  return true;
514  }
515  else
516  {
517  if (report)
518  {
519  Info<< "Face pyramids OK.\n" << endl;
520  }
521 
522  return false;
523  }
524 }
525 
526 
528 (
529  const bool report,
530  const scalar internalSkew,
531  const scalar boundarySkew,
532  const primitiveMesh& mesh,
533  const vectorField& cellCentres,
534  const vectorField& faceCentres,
535  const vectorField& faceAreas,
536  const labelList& checkFaces,
537  labelHashSet* setPtr
538 )
539 {
540  // Warn if the skew correction vector is more than skew times
541  // larger than the face area vector
542 
543  const labelList& own = mesh.faceOwner();
544  const labelList& nei = mesh.faceNeighbour();
545 
546  scalar maxSkew = 0;
547 
548  label nWarnSkew = 0;
549 
550  forAll(checkFaces, i)
551  {
552  label faceI = checkFaces[i];
553 
554  if (mesh.isInternalFace(faceI))
555  {
556  scalar dOwn = mag(faceCentres[faceI] - cellCentres[own[faceI]]);
557  scalar dNei = mag(faceCentres[faceI] - cellCentres[nei[faceI]]);
558 
559  point faceIntersection =
560  cellCentres[own[faceI]]*dNei/(dOwn+dNei)
561  + cellCentres[nei[faceI]]*dOwn/(dOwn+dNei);
562 
563  scalar skewness =
564  mag(faceCentres[faceI] - faceIntersection)
565  / (
566  mag(cellCentres[nei[faceI]]-cellCentres[own[faceI]])
567  + VSMALL
568  );
569 
570  // Check if the skewness vector is greater than the PN vector.
571  // This does not cause trouble but is a good indication of a poor
572  // mesh.
573  if (skewness > internalSkew)
574  {
575  if (report)
576  {
577  Pout<< "Severe skewness for face " << faceI
578  << " skewness = " << skewness << endl;
579  }
580 
581  if (setPtr)
582  {
583  setPtr->insert(faceI);
584  }
585 
586  nWarnSkew++;
587  }
588 
589  if (skewness > maxSkew)
590  {
591  maxSkew = skewness;
592  }
593  }
594  else
595  {
596  // Boundary faces: consider them to have only skewness error.
597  // (i.e. treat as if mirror cell on other side)
598 
599  vector faceNormal = faceAreas[faceI];
600  faceNormal /= mag(faceNormal) + VSMALL;
601 
602  vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
603 
604  vector dWall = faceNormal*(faceNormal & dOwn);
605 
606  point faceIntersection = cellCentres[own[faceI]] + dWall;
607 
608  scalar skewness =
609  mag(faceCentres[faceI] - faceIntersection)
610  /(2*mag(dWall) + VSMALL);
611 
612  // Check if the skewness vector is greater than the PN vector.
613  // This does not cause trouble but is a good indication of a poor
614  // mesh.
615  if (skewness > boundarySkew)
616  {
617  if (report)
618  {
619  Pout<< "Severe skewness for boundary face " << faceI
620  << " skewness = " << skewness << endl;
621  }
622 
623  if (setPtr)
624  {
625  setPtr->insert(faceI);
626  }
627 
628  nWarnSkew++;
629  }
630 
631  if (skewness > maxSkew)
632  {
633  maxSkew = skewness;
634  }
635  }
636  }
637 
638  reduce(maxSkew, maxOp<scalar>());
639  reduce(nWarnSkew, sumOp<label>());
640 
641  if (nWarnSkew > 0)
642  {
643  if (report)
644  {
645  WarningIn
646  (
647  "primitiveMeshGeometry::checkFaceSkewness"
648  "(const bool, const scalar, const labelList&, labelHashSet*)"
649  ) << "Large face skewness detected. Max skewness = "
650  << 100*maxSkew
651  << " percent.\nThis may impair the quality of the result." << nl
652  << nWarnSkew << " highly skew faces detected."
653  << endl;
654  }
655 
656  return true;
657  }
658  else
659  {
660  if (report)
661  {
662  Info<< "Max skewness = " << 100*maxSkew
663  << " percent. Face skewness OK.\n" << endl;
664  }
665 
666  return false;
667  }
668 }
669 
670 
672 (
673  const bool report,
674  const scalar warnWeight,
675  const primitiveMesh& mesh,
676  const vectorField& cellCentres,
677  const vectorField& faceCentres,
678  const vectorField& faceAreas,
679  const labelList& checkFaces,
680  labelHashSet* setPtr
681 )
682 {
683  // Warn if the delta factor (0..1) is too large.
684 
685  const labelList& own = mesh.faceOwner();
686  const labelList& nei = mesh.faceNeighbour();
687 
688  scalar minWeight = GREAT;
689 
690  label nWarnWeight = 0;
691 
692  forAll(checkFaces, i)
693  {
694  label faceI = checkFaces[i];
695 
696  if (mesh.isInternalFace(faceI))
697  {
698  const point& fc = faceCentres[faceI];
699 
700  scalar dOwn = mag(faceAreas[faceI] & (fc-cellCentres[own[faceI]]));
701  scalar dNei = mag(faceAreas[faceI] & (cellCentres[nei[faceI]]-fc));
702 
703  scalar weight = min(dNei,dOwn)/(dNei+dOwn);
704 
705  if (weight < warnWeight)
706  {
707  if (report)
708  {
709  Pout<< "Small weighting factor for face " << faceI
710  << " weight = " << weight << endl;
711  }
712 
713  if (setPtr)
714  {
715  setPtr->insert(faceI);
716  }
717 
718  nWarnWeight++;
719  }
720 
721  minWeight = min(minWeight, weight);
722  }
723  }
724 
725  reduce(minWeight, minOp<scalar>());
726  reduce(nWarnWeight, sumOp<label>());
727 
728  if (minWeight < warnWeight)
729  {
730  if (report)
731  {
732  WarningIn
733  (
734  "primitiveMeshGeometry::checkFaceWeights"
735  "(const bool, const scalar, const labelList&, labelHashSet*)"
736  ) << "Small interpolation weight detected. Min weight = "
737  << minWeight << '.' << nl
738  << nWarnWeight << " faces with small weights detected."
739  << endl;
740  }
741 
742  return true;
743  }
744  else
745  {
746  if (report)
747  {
748  Info<< "Min weight = " << minWeight
749  << " percent. Weights OK.\n" << endl;
750  }
751 
752  return false;
753  }
754 }
755 
756 
757 // Check convexity of angles in a face. Allow a slight non-convexity.
758 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
759 // (if truly concave and points not visible from face centre the face-pyramid
760 // check in checkMesh will fail)
762 (
763  const bool report,
764  const scalar maxDeg,
765  const primitiveMesh& mesh,
766  const vectorField& faceAreas,
767  const pointField& p,
768  const labelList& checkFaces,
769  labelHashSet* setPtr
770 )
771 {
772  if (maxDeg < -SMALL || maxDeg > 180+SMALL)
773  {
775  (
776  "primitiveMeshGeometry::checkFaceAngles"
777  "(const bool, const scalar, const pointField&, const labelList&"
778  ", labelHashSet*)"
779  ) << "maxDeg should be [0..180] but is now " << maxDeg
780  << abort(FatalError);
781  }
782 
783  const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
784 
785  const faceList& fcs = mesh.faces();
786 
787  scalar maxEdgeSin = 0.0;
788 
789  label nConcave = 0;
790 
791  label errorFaceI = -1;
792 
793  forAll(checkFaces, i)
794  {
795  label faceI = checkFaces[i];
796 
797  const face& f = fcs[faceI];
798 
799  vector faceNormal = faceAreas[faceI];
800  faceNormal /= mag(faceNormal) + VSMALL;
801 
802  // Get edge from f[0] to f[size-1];
803  vector ePrev(p[f[0]] - p[f[f.size()-1]]);
804  scalar magEPrev = mag(ePrev);
805  ePrev /= magEPrev + VSMALL;
806 
807  forAll(f, fp0)
808  {
809  // Get vertex after fp
810  label fp1 = f.fcIndex(fp0);
811 
812  // Normalized vector between two consecutive points
813  vector e10(p[f[fp1]] - p[f[fp0]]);
814  scalar magE10 = mag(e10);
815  e10 /= magE10 + VSMALL;
816 
817  if (magEPrev > SMALL && magE10 > SMALL)
818  {
819  vector edgeNormal = ePrev ^ e10;
820  scalar magEdgeNormal = mag(edgeNormal);
821 
822  if (magEdgeNormal < maxSin)
823  {
824  // Edges (almost) aligned -> face is ok.
825  }
826  else
827  {
828  // Check normal
829  edgeNormal /= magEdgeNormal;
830 
831  if ((edgeNormal & faceNormal) < SMALL)
832  {
833  if (faceI != errorFaceI)
834  {
835  // Count only one error per face.
836  errorFaceI = faceI;
837  nConcave++;
838  }
839 
840  if (setPtr)
841  {
842  setPtr->insert(faceI);
843  }
844 
845  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
846  }
847  }
848  }
849 
850  ePrev = e10;
851  magEPrev = magE10;
852  }
853  }
854 
855  reduce(nConcave, sumOp<label>());
856  reduce(maxEdgeSin, maxOp<scalar>());
857 
858  if (report)
859  {
860  if (maxEdgeSin > SMALL)
861  {
862  scalar maxConcaveDegr =
863  Foam::asin(Foam::min(1.0, maxEdgeSin))
864  * 180.0/mathematicalConstant::pi;
865 
866  Info<< "There are " << nConcave
867  << " faces with concave angles between consecutive"
868  << " edges. Max concave angle = " << maxConcaveDegr
869  << " degrees.\n" << endl;
870  }
871  else
872  {
873  Info<< "All angles in faces are convex or less than " << maxDeg
874  << " degrees concave.\n" << endl;
875  }
876  }
877 
878  if (nConcave > 0)
879  {
880  if (report)
881  {
882  WarningIn
883  (
884  "primitiveMeshGeometry::checkFaceAngles"
885  "(const bool, const scalar, const pointField&"
886  ", const labelList&, labelHashSet*)"
887  ) << nConcave << " face points with severe concave angle (> "
888  << maxDeg << " deg) found.\n"
889  << endl;
890  }
891 
892  return true;
893  }
894  else
895  {
896  return false;
897  }
898 }
899 
900 
904 //bool Foam::primitiveMeshGeometry::checkFaceFlatness
905 //(
906 // const bool report,
907 // const scalar warnFlatness,
908 // const primitiveMesh& mesh,
909 // const vectorField& faceAreas,
910 // const vectorField& faceCentres,
911 // const pointField& p,
912 // const labelList& checkFaces,
913 // labelHashSet* setPtr
914 //)
915 //{
916 // if (warnFlatness < 0 || warnFlatness > 1)
917 // {
918 // FatalErrorIn
919 // (
920 // "primitiveMeshGeometry::checkFaceFlatness"
921 // "(const bool, const scalar, const pointField&"
922 // ", const labelList&, labelHashSet*)"
923 // ) << "warnFlatness should be [0..1] but is now " << warnFlatness
924 // << abort(FatalError);
925 // }
926 //
927 //
928 // const faceList& fcs = mesh.faces();
929 //
930 // // Areas are calculated as the sum of areas. (see
931 // // primitiveMeshFaceCentresAndAreas.C)
932 //
933 // label nWarped = 0;
934 //
935 // scalar minFlatness = GREAT;
936 // scalar sumFlatness = 0;
937 // label nSummed = 0;
938 //
939 // forAll(checkFaces, i)
940 // {
941 // label faceI = checkFaces[i];
942 //
943 // const face& f = fcs[faceI];
944 //
945 // scalar magArea = mag(faceAreas[faceI]);
946 //
947 // if (f.size() > 3 && magArea > VSMALL)
948 // {
949 // const point& fc = faceCentres[faceI];
950 //
951 // // Calculate the sum of magnitude of areas and compare to magnitude
952 // // of sum of areas.
953 //
954 // scalar sumA = 0.0;
955 //
956 // forAll(f, fp)
957 // {
958 // const point& thisPoint = p[f[fp]];
959 // const point& nextPoint = p[f.nextLabel(fp)];
960 //
961 // // Triangle around fc.
962 // vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
963 // sumA += mag(n);
964 // }
965 //
966 // scalar flatness = magArea / (sumA+VSMALL);
967 //
968 // sumFlatness += flatness;
969 // nSummed++;
970 //
971 // minFlatness = min(minFlatness, flatness);
972 //
973 // if (flatness < warnFlatness)
974 // {
975 // nWarped++;
976 //
977 // if (setPtr)
978 // {
979 // setPtr->insert(faceI);
980 // }
981 // }
982 // }
983 // }
984 //
985 //
986 // reduce(nWarped, sumOp<label>());
987 // reduce(minFlatness, minOp<scalar>());
988 //
989 // reduce(nSummed, sumOp<label>());
990 // reduce(sumFlatness, sumOp<scalar>());
991 //
992 // if (report)
993 // {
994 // if (nSummed > 0)
995 // {
996 // Info<< "Face flatness (1 = flat, 0 = butterfly) : average = "
997 // << sumFlatness / nSummed << " min = " << minFlatness << endl;
998 // }
999 //
1000 // if (nWarped> 0)
1001 // {
1002 // Info<< "There are " << nWarped
1003 // << " faces with ratio between projected and actual area < "
1004 // << warnFlatness
1005 // << ".\nMinimum ratio (minimum flatness, maximum warpage) = "
1006 // << minFlatness << nl << endl;
1007 // }
1008 // else
1009 // {
1010 // Info<< "All faces are flat in that the ratio between projected"
1011 // << " and actual area is > " << warnFlatness << nl << endl;
1012 // }
1013 // }
1014 //
1015 // if (nWarped > 0)
1016 // {
1017 // if (report)
1018 // {
1019 // WarningIn
1020 // (
1021 // "primitiveMeshGeometry::checkFaceFlatness"
1022 // "(const bool, const scalar, const pointField&"
1023 // ", const labelList&, labelHashSet*)"
1024 // ) << nWarped << " faces with severe warpage (flatness < "
1025 // << warnFlatness << ") found.\n"
1026 // << endl;
1027 // }
1028 //
1029 // return true;
1030 // }
1031 // else
1032 // {
1033 // return false;
1034 // }
1035 //}
1036 
1037 
1038 // Check twist of faces. Is calculated as the difference between areas of
1039 // individual triangles and the overall area of the face (which ifself is
1040 // is the average of the areas of the individual triangles).
1043  const bool report,
1044  const scalar minTwist,
1045  const primitiveMesh& mesh,
1046  const vectorField& faceAreas,
1047  const vectorField& faceCentres,
1048  const pointField& p,
1049  const labelList& checkFaces,
1050  labelHashSet* setPtr
1051 )
1052 {
1053  if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1054  {
1055  FatalErrorIn
1056  (
1057  "primitiveMeshGeometry::checkFaceTwist"
1058  "(const bool, const scalar, const primitiveMesh&, const pointField&"
1059  ", const labelList&, labelHashSet*)"
1060  ) << "minTwist should be [-1..1] but is now " << minTwist
1061  << abort(FatalError);
1062  }
1063 
1064 
1065  const faceList& fcs = mesh.faces();
1066 
1067  // Areas are calculated as the sum of areas. (see
1068  // primitiveMeshFaceCentresAndAreas.C)
1069 
1070  label nWarped = 0;
1071 
1072  forAll(checkFaces, i)
1073  {
1074  label faceI = checkFaces[i];
1075 
1076  const face& f = fcs[faceI];
1077 
1078  scalar magArea = mag(faceAreas[faceI]);
1079 
1080  if (f.size() > 3 && magArea > VSMALL)
1081  {
1082  const vector nf = faceAreas[faceI] / magArea;
1083 
1084  const point& fc = faceCentres[faceI];
1085 
1086  forAll(f, fpI)
1087  {
1088  vector triArea
1089  (
1090  triPointRef
1091  (
1092  p[f[fpI]],
1093  p[f.nextLabel(fpI)],
1094  fc
1095  ).normal()
1096  );
1097 
1098  scalar magTri = mag(triArea);
1099 
1100  if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1101  {
1102  nWarped++;
1103 
1104  if (setPtr)
1105  {
1106  setPtr->insert(faceI);
1107  }
1108  }
1109  }
1110  }
1111  }
1112 
1113 
1114  reduce(nWarped, sumOp<label>());
1115 
1116  if (report)
1117  {
1118  if (nWarped> 0)
1119  {
1120  Info<< "There are " << nWarped
1121  << " faces with cosine of the angle"
1122  << " between triangle normal and face normal less than "
1123  << minTwist << nl << endl;
1124  }
1125  else
1126  {
1127  Info<< "All faces are flat in that the cosine of the angle"
1128  << " between triangle normal and face normal less than "
1129  << minTwist << nl << endl;
1130  }
1131  }
1132 
1133  if (nWarped > 0)
1134  {
1135  if (report)
1136  {
1137  WarningIn
1138  (
1139  "primitiveMeshGeometry::checkFaceTwist"
1140  "(const bool, const scalar, const primitiveMesh&"
1141  ", const pointField&, const labelList&, labelHashSet*)"
1142  ) << nWarped << " faces with severe warpage "
1143  << "(cosine of the angle between triangle normal and face normal"
1144  << " < " << minTwist << ") found.\n"
1145  << endl;
1146  }
1147 
1148  return true;
1149  }
1150  else
1151  {
1152  return false;
1153  }
1154 }
1155 
1156 
1159  const bool report,
1160  const scalar minArea,
1161  const primitiveMesh& mesh,
1162  const vectorField& faceAreas,
1163  const labelList& checkFaces,
1164  labelHashSet* setPtr
1165 )
1166 {
1167  label nZeroArea = 0;
1168 
1169  forAll(checkFaces, i)
1170  {
1171  label faceI = checkFaces[i];
1172 
1173  if (mag(faceAreas[faceI]) < minArea)
1174  {
1175  if (setPtr)
1176  {
1177  setPtr->insert(faceI);
1178  }
1179  nZeroArea++;
1180  }
1181  }
1182 
1183 
1184  reduce(nZeroArea, sumOp<label>());
1185 
1186  if (report)
1187  {
1188  if (nZeroArea > 0)
1189  {
1190  Info<< "There are " << nZeroArea
1191  << " faces with area < " << minArea << '.' << nl << endl;
1192  }
1193  else
1194  {
1195  Info<< "All faces have area > " << minArea << '.' << nl << endl;
1196  }
1197  }
1198 
1199  if (nZeroArea > 0)
1200  {
1201  if (report)
1202  {
1203  WarningIn
1204  (
1205  "primitiveMeshGeometry::checkFaceArea"
1206  "(const bool, const scalar, const primitiveMesh&"
1207  ", const pointField&, const labelList&, labelHashSet*)"
1208  ) << nZeroArea << " faces with area < " << minArea
1209  << " found.\n"
1210  << endl;
1211  }
1212 
1213  return true;
1214  }
1215  else
1216  {
1217  return false;
1218  }
1219 }
1220 
1221 
1224  const bool report,
1225  const scalar warnDet,
1226  const primitiveMesh& mesh,
1227  const vectorField& faceAreas,
1228  const labelList& checkFaces,
1229  const labelList& affectedCells,
1230  labelHashSet* setPtr
1231 )
1232 {
1233  const cellList& cells = mesh.cells();
1234 
1235  scalar minDet = GREAT;
1236  scalar sumDet = 0.0;
1237  label nSumDet = 0;
1238  label nWarnDet = 0;
1239 
1240  forAll(affectedCells, i)
1241  {
1242  const cell& cFaces = cells[affectedCells[i]];
1243 
1244  tensor areaSum(tensor::zero);
1245  scalar magAreaSum = 0;
1246 
1247  forAll(cFaces, cFaceI)
1248  {
1249  label faceI = cFaces[cFaceI];
1250 
1251  scalar magArea = mag(faceAreas[faceI]);
1252 
1253  magAreaSum += magArea;
1254  areaSum += faceAreas[faceI]*(faceAreas[faceI]/magArea);
1255  }
1256 
1257  scalar scaledDet = det(areaSum/magAreaSum)/0.037037037037037;
1258 
1259  minDet = min(minDet, scaledDet);
1260  sumDet += scaledDet;
1261  nSumDet++;
1262 
1263  if (scaledDet < warnDet)
1264  {
1265  if (setPtr)
1266  {
1267  // Insert all faces of the cell.
1268  forAll(cFaces, cFaceI)
1269  {
1270  label faceI = cFaces[cFaceI];
1271  setPtr->insert(faceI);
1272  }
1273  }
1274  nWarnDet++;
1275  }
1276  }
1277 
1278  reduce(minDet, minOp<scalar>());
1279  reduce(sumDet, sumOp<scalar>());
1280  reduce(nSumDet, sumOp<label>());
1281  reduce(nWarnDet, sumOp<label>());
1282 
1283  if (report)
1284  {
1285  if (nSumDet > 0)
1286  {
1287  Info<< "Cell determinant (1 = uniform cube) : average = "
1288  << sumDet / nSumDet << " min = " << minDet << endl;
1289  }
1290 
1291  if (nWarnDet > 0)
1292  {
1293  Info<< "There are " << nWarnDet
1294  << " cells with determinant < " << warnDet << '.' << nl
1295  << endl;
1296  }
1297  else
1298  {
1299  Info<< "All faces have determinant > " << warnDet << '.' << nl
1300  << endl;
1301  }
1302  }
1303 
1304  if (nWarnDet > 0)
1305  {
1306  if (report)
1307  {
1308  WarningIn
1309  (
1310  "primitiveMeshGeometry::checkCellDeterminant"
1311  "(const bool, const scalar, const primitiveMesh&"
1312  ", const pointField&, const labelList&, const labelList&"
1313  ", labelHashSet*)"
1314  ) << nWarnDet << " cells with determinant < " << warnDet
1315  << " found.\n"
1316  << endl;
1317  }
1318 
1319  return true;
1320  }
1321  else
1322  {
1323  return false;
1324  }
1325 }
1326 
1327 
1330  const bool report,
1331  const scalar orthWarn,
1332  const labelList& checkFaces,
1333  labelHashSet* setPtr
1334 ) const
1335 {
1336  return checkFaceDotProduct
1337  (
1338  report,
1339  orthWarn,
1340  mesh_,
1341  cellCentres_,
1342  faceAreas_,
1343  checkFaces,
1344  setPtr
1345  );
1346 }
1347 
1348 
1351  const bool report,
1352  const scalar minPyrVol,
1353  const pointField& p,
1354  const labelList& checkFaces,
1355  labelHashSet* setPtr
1356 ) const
1357 {
1358  return checkFacePyramids
1359  (
1360  report,
1361  minPyrVol,
1362  mesh_,
1363  cellCentres_,
1364  p,
1365  checkFaces,
1366  setPtr
1367  );
1368 }
1369 
1370 
1373  const bool report,
1374  const scalar internalSkew,
1375  const scalar boundarySkew,
1376  const labelList& checkFaces,
1377  labelHashSet* setPtr
1378 ) const
1379 {
1380  return checkFaceSkewness
1381  (
1382  report,
1383  internalSkew,
1384  boundarySkew,
1385  mesh_,
1386  cellCentres_,
1387  faceCentres_,
1388  faceAreas_,
1389  checkFaces,
1390  setPtr
1391  );
1392 }
1393 
1394 
1397  const bool report,
1398  const scalar warnWeight,
1399  const labelList& checkFaces,
1400  labelHashSet* setPtr
1401 ) const
1402 {
1403  return checkFaceWeights
1404  (
1405  report,
1406  warnWeight,
1407  mesh_,
1408  cellCentres_,
1409  faceCentres_,
1410  faceAreas_,
1411  checkFaces,
1412  setPtr
1413  );
1414 }
1415 
1416 
1419  const bool report,
1420  const scalar maxDeg,
1421  const pointField& p,
1422  const labelList& checkFaces,
1423  labelHashSet* setPtr
1424 ) const
1425 {
1426  return checkFaceAngles
1427  (
1428  report,
1429  maxDeg,
1430  mesh_,
1431  faceAreas_,
1432  p,
1433  checkFaces,
1434  setPtr
1435  );
1436 }
1437 
1438 
1439 //bool Foam::primitiveMeshGeometry::checkFaceFlatness
1440 //(
1441 // const bool report,
1442 // const scalar warnFlatness,
1443 // const pointField& p,
1444 // const labelList& checkFaces,
1445 // labelHashSet* setPtr
1446 //) const
1447 //{
1448 // return checkFaceFlatness
1449 // (
1450 // report,
1451 // warnFlatness,
1452 // mesh_,
1453 // faceAreas_,
1454 // faceCentres_,
1455 // p,
1456 // checkFaces,
1457 // setPtr
1458 // );
1459 //}
1460 
1461 
1464  const bool report,
1465  const scalar minTwist,
1466  const pointField& p,
1467  const labelList& checkFaces,
1468  labelHashSet* setPtr
1469 ) const
1470 {
1471  return checkFaceTwist
1472  (
1473  report,
1474  minTwist,
1475  mesh_,
1476  faceAreas_,
1477  faceCentres_,
1478  p,
1479  checkFaces,
1480  setPtr
1481  );
1482 }
1483 
1484 
1487  const bool report,
1488  const scalar minArea,
1489  const labelList& checkFaces,
1490  labelHashSet* setPtr
1491 ) const
1492 {
1493  return checkFaceArea
1494  (
1495  report,
1496  minArea,
1497  mesh_,
1498  faceAreas_,
1499  checkFaces,
1500  setPtr
1501  );
1502 }
1503 
1504 
1507  const bool report,
1508  const scalar warnDet,
1509  const labelList& checkFaces,
1510  const labelList& affectedCells,
1511  labelHashSet* setPtr
1512 ) const
1513 {
1514  return checkCellDeterminant
1515  (
1516  report,
1517  warnDet,
1518  mesh_,
1519  faceAreas_,
1520  checkFaces,
1521  affectedCells,
1522  setPtr
1523  );
1524 }
1525 
1526 
1527 // ************************ vim: set sw=4 sts=4 et: ************************ //