FreeFOAM The Cross-Platform CFD Toolkit
polyMeshGeometry.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 "polyMeshGeometry.H"
28 #include <OpenFOAM/syncTools.H>
29 
30 namespace Foam
31 {
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 defineTypeNameAndDebug(polyMeshGeometry, 0);
36 
37 }
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::polyMeshGeometry::updateFaceCentresAndAreas
43 (
44  const pointField& p,
45  const labelList& changedFaces
46 )
47 {
48  const faceList& fs = mesh_.faces();
49 
50  forAll(changedFaces, i)
51  {
52  label facei = changedFaces[i];
53 
54  const labelList& f = fs[facei];
55  label nPoints = f.size();
56 
57  // If the face is a triangle, do a direct calculation for efficiency
58  // and to avoid round-off error-related problems
59  if (nPoints == 3)
60  {
61  faceCentres_[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
62  faceAreas_[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
63  }
64  else
65  {
66  vector sumN = vector::zero;
67  scalar sumA = 0.0;
68  vector sumAc = vector::zero;
69 
70  point fCentre = p[f[0]];
71  for (label pi = 1; pi < nPoints; pi++)
72  {
73  fCentre += p[f[pi]];
74  }
75 
76  fCentre /= nPoints;
77 
78  for (label pi = 0; pi < nPoints; pi++)
79  {
80  const point& nextPoint = p[f[(pi + 1) % nPoints]];
81 
82  vector c = p[f[pi]] + nextPoint + fCentre;
83  vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]);
84  scalar a = mag(n);
85 
86  sumN += n;
87  sumA += a;
88  sumAc += a*c;
89  }
90 
91  faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL);
92  faceAreas_[facei] = 0.5*sumN;
93  }
94  }
95 }
96 
97 
98 void Foam::polyMeshGeometry::updateCellCentresAndVols
99 (
100  const labelList& changedCells,
101  const labelList& changedFaces
102 )
103 {
104  // Clear the fields for accumulation
105  UIndirectList<vector>(cellCentres_, changedCells) = vector::zero;
106  UIndirectList<scalar>(cellVolumes_, changedCells) = 0.0;
107 
108  const labelList& own = mesh_.faceOwner();
109  const labelList& nei = mesh_.faceNeighbour();
110 
111  // first estimate the approximate cell centre as the average of face centres
112 
113  vectorField cEst(mesh_.nCells());
114  UIndirectList<vector>(cEst, changedCells) = vector::zero;
115  scalarField nCellFaces(mesh_.nCells());
116  UIndirectList<scalar>(nCellFaces, changedCells) = 0.0;
117 
118  forAll(changedFaces, i)
119  {
120  label faceI = changedFaces[i];
121  cEst[own[faceI]] += faceCentres_[faceI];
122  nCellFaces[own[faceI]] += 1;
123 
124  if (mesh_.isInternalFace(faceI))
125  {
126  cEst[nei[faceI]] += faceCentres_[faceI];
127  nCellFaces[nei[faceI]] += 1;
128  }
129  }
130 
131  forAll(changedCells, i)
132  {
133  label cellI = changedCells[i];
134  cEst[cellI] /= nCellFaces[cellI];
135  }
136 
137  forAll(changedFaces, i)
138  {
139  label faceI = changedFaces[i];
140 
141  // Calculate 3*face-pyramid volume
142  scalar pyr3Vol = max
143  (
144  faceAreas_[faceI] & (faceCentres_[faceI] - cEst[own[faceI]]),
145  VSMALL
146  );
147 
148  // Calculate face-pyramid centre
149  vector pc = (3.0/4.0)*faceCentres_[faceI] + (1.0/4.0)*cEst[own[faceI]];
150 
151  // Accumulate volume-weighted face-pyramid centre
152  cellCentres_[own[faceI]] += pyr3Vol*pc;
153 
154  // Accumulate face-pyramid volume
155  cellVolumes_[own[faceI]] += pyr3Vol;
156 
157  if (mesh_.isInternalFace(faceI))
158  {
159  // Calculate 3*face-pyramid volume
160  scalar pyr3Vol = max
161  (
162  faceAreas_[faceI] & (cEst[nei[faceI]] - faceCentres_[faceI]),
163  VSMALL
164  );
165 
166  // Calculate face-pyramid centre
167  vector pc =
168  (3.0/4.0)*faceCentres_[faceI]
169  + (1.0/4.0)*cEst[nei[faceI]];
170 
171  // Accumulate volume-weighted face-pyramid centre
172  cellCentres_[nei[faceI]] += pyr3Vol*pc;
173 
174  // Accumulate face-pyramid volume
175  cellVolumes_[nei[faceI]] += pyr3Vol;
176  }
177  }
178 
179  forAll(changedCells, i)
180  {
181  label cellI = changedCells[i];
182 
183  cellCentres_[cellI] /= cellVolumes_[cellI] + VSMALL;
184  cellVolumes_[cellI] *= (1.0/3.0);
185  }
186 }
187 
188 
190 (
191  const polyMesh& mesh,
192  const labelList& changedFaces
193 )
194 {
195  const labelList& own = mesh.faceOwner();
196  const labelList& nei = mesh.faceNeighbour();
197 
198  labelHashSet affectedCells(2*changedFaces.size());
199 
200  forAll(changedFaces, i)
201  {
202  label faceI = changedFaces[i];
203 
204  affectedCells.insert(own[faceI]);
205 
206  if (mesh.isInternalFace(faceI))
207  {
208  affectedCells.insert(nei[faceI]);
209  }
210  }
211  return affectedCells.toc();
212 }
213 
214 
215 Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
216 (
217  const polyMesh& mesh,
218  const bool report,
219  const scalar severeNonorthogonalityThreshold,
220  const label faceI,
221  const vector& s, // face area vector
222  const vector& d, // cc-cc vector
223 
224  label& severeNonOrth,
225  label& errorNonOrth,
226  labelHashSet* setPtr
227 )
228 {
229  scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
230 
231  if (dDotS < severeNonorthogonalityThreshold)
232  {
233  label nei = -1;
234 
235  if (mesh.isInternalFace(faceI))
236  {
237  nei = mesh.faceNeighbour()[faceI];
238  }
239 
240  if (dDotS > SMALL)
241  {
242  if (report)
243  {
244  // Severe non-orthogonality but mesh still OK
245  Pout<< "Severe non-orthogonality for face " << faceI
246  << " between cells " << mesh.faceOwner()[faceI]
247  << " and " << nei
248  << ": Angle = "
249  << ::acos(dDotS)/mathematicalConstant::pi*180.0
250  << " deg." << endl;
251  }
252 
253  severeNonOrth++;
254  }
255  else
256  {
257  // Non-orthogonality greater than 90 deg
258  if (report)
259  {
260  WarningIn
261  (
262  "polyMeshGeometry::checkFaceDotProduct"
263  "(const bool, const scalar, const labelList&"
264  ", labelHashSet*)"
265  ) << "Severe non-orthogonality detected for face "
266  << faceI
267  << " between cells " << mesh.faceOwner()[faceI]
268  << " and " << nei
269  << ": Angle = "
270  << ::acos(dDotS)/mathematicalConstant::pi*180.0
271  << " deg." << endl;
272  }
273 
274  errorNonOrth++;
275  }
276 
277  if (setPtr)
278  {
279  setPtr->insert(faceI);
280  }
281  }
282  return dDotS;
283 }
284 
285 
286 Foam::scalar Foam::polyMeshGeometry::calcSkewness
287 (
288  const point& ownCc,
289  const point& neiCc,
290  const point& fc
291 )
292 {
293  scalar dOwn = mag(fc - ownCc);
294  scalar dNei = mag(fc - neiCc);
295 
296  point faceIntersection =
297  ownCc*dNei/(dOwn+dNei+VSMALL)
298  + neiCc*dOwn/(dOwn+dNei+VSMALL);
299 
300  return
301  mag(fc - faceIntersection)
302  / (
303  mag(neiCc-ownCc)
304  + VSMALL
305  );
306 }
307 
308 
309 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
310 
311 // Construct from components
313 :
314  mesh_(mesh)
315 {
316  correct();
317 }
318 
319 
320 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
321 
322 
323 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
324 
325 //- Take over properties from mesh
327 {
328  faceAreas_ = mesh_.faceAreas();
329  faceCentres_ = mesh_.faceCentres();
330  cellCentres_ = mesh_.cellCentres();
331  cellVolumes_ = mesh_.cellVolumes();
332 }
333 
334 
335 //- Recalculate on selected faces
337 (
338  const pointField& p,
339  const labelList& changedFaces
340 )
341 {
342  // Update face quantities
343  updateFaceCentresAndAreas(p, changedFaces);
344  // Update cell quantities from face quantities
345  updateCellCentresAndVols(affectedCells(mesh_, changedFaces), changedFaces);
346 }
347 
348 
350 (
351  const bool report,
352  const scalar orthWarn,
353  const polyMesh& mesh,
354  const vectorField& cellCentres,
355  const vectorField& faceAreas,
356  const labelList& checkFaces,
357  const List<labelPair>& baffles,
358  labelHashSet* setPtr
359 )
360 {
361  // for all internal and coupled faces check theat the d dot S product
362  // is positive
363 
364  const labelList& own = mesh.faceOwner();
365  const labelList& nei = mesh.faceNeighbour();
366  const polyBoundaryMesh& patches = mesh.boundaryMesh();
367 
368  // Severe nonorthogonality threshold
369  const scalar severeNonorthogonalityThreshold =
370  ::cos(orthWarn/180.0*mathematicalConstant::pi);
371 
372 
373  // Calculate coupled cell centre
374  vectorField neiCc(mesh.nFaces()-mesh.nInternalFaces());
375 
376  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
377  {
378  neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
379  }
380  syncTools::swapBoundaryFaceList(mesh, neiCc, true);
381 
382 
383  scalar minDDotS = GREAT;
384 
385  scalar sumDDotS = 0;
386  label nDDotS = 0;
387 
388  label severeNonOrth = 0;
389 
390  label errorNonOrth = 0;
391 
392  forAll(checkFaces, i)
393  {
394  label faceI = checkFaces[i];
395 
396  const point& ownCc = cellCentres[own[faceI]];
397 
398  if (mesh.isInternalFace(faceI))
399  {
400  scalar dDotS = checkNonOrtho
401  (
402  mesh,
403  report,
404  severeNonorthogonalityThreshold,
405  faceI,
406  faceAreas[faceI],
407  cellCentres[nei[faceI]] - ownCc,
408 
409  severeNonOrth,
410  errorNonOrth,
411  setPtr
412  );
413 
414  if (dDotS < minDDotS)
415  {
416  minDDotS = dDotS;
417  }
418 
419  sumDDotS += dDotS;
420  nDDotS++;
421  }
422  else
423  {
424  label patchI = patches.whichPatch(faceI);
425 
426  if (patches[patchI].coupled())
427  {
428  scalar dDotS = checkNonOrtho
429  (
430  mesh,
431  report,
432  severeNonorthogonalityThreshold,
433  faceI,
434  faceAreas[faceI],
435  neiCc[faceI-mesh.nInternalFaces()] - ownCc,
436 
437  severeNonOrth,
438  errorNonOrth,
439  setPtr
440  );
441 
442  if (dDotS < minDDotS)
443  {
444  minDDotS = dDotS;
445  }
446 
447  sumDDotS += dDotS;
448  nDDotS++;
449  }
450  }
451  }
452 
453  forAll(baffles, i)
454  {
455  label face0 = baffles[i].first();
456  label face1 = baffles[i].second();
457 
458  const point& ownCc = cellCentres[own[face0]];
459 
460  scalar dDotS = checkNonOrtho
461  (
462  mesh,
463  report,
464  severeNonorthogonalityThreshold,
465  face0,
466  faceAreas[face0],
467  cellCentres[own[face1]] - ownCc,
468 
469  severeNonOrth,
470  errorNonOrth,
471  setPtr
472  );
473 
474  if (dDotS < minDDotS)
475  {
476  minDDotS = dDotS;
477  }
478 
479  sumDDotS += dDotS;
480  nDDotS++;
481  }
482 
483  reduce(minDDotS, minOp<scalar>());
484  reduce(sumDDotS, sumOp<scalar>());
485  reduce(nDDotS, sumOp<label>());
486  reduce(severeNonOrth, sumOp<label>());
487  reduce(errorNonOrth, sumOp<label>());
488 
489  // Only report if there are some internal faces
490  if (nDDotS > 0)
491  {
492  if (report && minDDotS < severeNonorthogonalityThreshold)
493  {
494  Info<< "Number of non-orthogonality errors: " << errorNonOrth
495  << ". Number of severely non-orthogonal faces: "
496  << severeNonOrth << "." << endl;
497  }
498  }
499 
500  if (report)
501  {
502  if (nDDotS > 0)
503  {
504  Info<< "Mesh non-orthogonality Max: "
505  << ::acos(minDDotS)/mathematicalConstant::pi*180.0
506  << " average: " <<
507  ::acos(sumDDotS/nDDotS)/mathematicalConstant::pi*180.0
508  << endl;
509  }
510  }
511 
512  if (errorNonOrth > 0)
513  {
514  if (report)
515  {
517  (
518  "polyMeshGeometry::checkFaceDotProduct"
519  "(const bool, const scalar, const labelList&, labelHashSet*)"
520  ) << "Error in non-orthogonality detected" << endl;
521  }
522 
523  return true;
524  }
525  else
526  {
527  if (report)
528  {
529  Info<< "Non-orthogonality check OK.\n" << endl;
530  }
531 
532  return false;
533  }
534 }
535 
536 
538 (
539  const bool report,
540  const scalar minPyrVol,
541  const polyMesh& mesh,
542  const vectorField& cellCentres,
543  const pointField& p,
544  const labelList& checkFaces,
545  const List<labelPair>& baffles,
546  labelHashSet* setPtr
547 )
548 {
549  // check whether face area vector points to the cell with higher label
550  const labelList& own = mesh.faceOwner();
551  const labelList& nei = mesh.faceNeighbour();
552 
553  const faceList& f = mesh.faces();
554 
555  label nErrorPyrs = 0;
556 
557  forAll(checkFaces, i)
558  {
559  label faceI = checkFaces[i];
560 
561  // Create the owner pyramid - it will have negative volume
562  scalar pyrVol = pyramidPointFaceRef
563  (
564  f[faceI],
565  cellCentres[own[faceI]]
566  ).mag(p);
567 
568  if (pyrVol > -minPyrVol)
569  {
570  if (report)
571  {
572  Pout<< "bool polyMeshGeometry::checkFacePyramids("
573  << "const bool, const scalar, const pointField&"
574  << ", const labelList&, labelHashSet*): "
575  << "face " << faceI << " points the wrong way. " << endl
576  << "Pyramid volume: " << -pyrVol
577  << " Face " << f[faceI] << " area: " << f[faceI].mag(p)
578  << " Owner cell: " << own[faceI] << endl
579  << "Owner cell vertex labels: "
580  << mesh.cells()[own[faceI]].labels(f)
581  << endl;
582  }
583 
584 
585  if (setPtr)
586  {
587  setPtr->insert(faceI);
588  }
589 
590  nErrorPyrs++;
591  }
592 
593  if (mesh.isInternalFace(faceI))
594  {
595  // Create the neighbour pyramid - it will have positive volume
596  scalar pyrVol =
597  pyramidPointFaceRef(f[faceI], cellCentres[nei[faceI]]).mag(p);
598 
599  if (pyrVol < minPyrVol)
600  {
601  if (report)
602  {
603  Pout<< "bool polyMeshGeometry::checkFacePyramids("
604  << "const bool, const scalar, const pointField&"
605  << ", const labelList&, labelHashSet*): "
606  << "face " << faceI << " points the wrong way. " << endl
607  << "Pyramid volume: " << -pyrVol
608  << " Face " << f[faceI] << " area: " << f[faceI].mag(p)
609  << " Neighbour cell: " << nei[faceI] << endl
610  << "Neighbour cell vertex labels: "
611  << mesh.cells()[nei[faceI]].labels(f)
612  << endl;
613  }
614 
615  if (setPtr)
616  {
617  setPtr->insert(faceI);
618  }
619 
620  nErrorPyrs++;
621  }
622  }
623  }
624 
625  forAll(baffles, i)
626  {
627  label face0 = baffles[i].first();
628  label face1 = baffles[i].second();
629 
630  const point& ownCc = cellCentres[own[face0]];
631 
632  // Create the owner pyramid - it will have negative volume
633  scalar pyrVolOwn = pyramidPointFaceRef
634  (
635  f[face0],
636  ownCc
637  ).mag(p);
638 
639  if (pyrVolOwn > -minPyrVol)
640  {
641  if (report)
642  {
643  Pout<< "bool polyMeshGeometry::checkFacePyramids("
644  << "const bool, const scalar, const pointField&"
645  << ", const labelList&, labelHashSet*): "
646  << "face " << face0 << " points the wrong way. " << endl
647  << "Pyramid volume: " << -pyrVolOwn
648  << " Face " << f[face0] << " area: " << f[face0].mag(p)
649  << " Owner cell: " << own[face0] << endl
650  << "Owner cell vertex labels: "
651  << mesh.cells()[own[face0]].labels(f)
652  << endl;
653  }
654 
655 
656  if (setPtr)
657  {
658  setPtr->insert(face0);
659  }
660 
661  nErrorPyrs++;
662  }
663 
664  // Create the neighbour pyramid - it will have positive volume
665  scalar pyrVolNbr =
666  pyramidPointFaceRef(f[face0], cellCentres[own[face1]]).mag(p);
667 
668  if (pyrVolNbr < minPyrVol)
669  {
670  if (report)
671  {
672  Pout<< "bool polyMeshGeometry::checkFacePyramids("
673  << "const bool, const scalar, const pointField&"
674  << ", const labelList&, labelHashSet*): "
675  << "face " << face0 << " points the wrong way. " << endl
676  << "Pyramid volume: " << -pyrVolNbr
677  << " Face " << f[face0] << " area: " << f[face0].mag(p)
678  << " Neighbour cell: " << own[face1] << endl
679  << "Neighbour cell vertex labels: "
680  << mesh.cells()[own[face1]].labels(f)
681  << endl;
682  }
683 
684  if (setPtr)
685  {
686  setPtr->insert(face0);
687  }
688 
689  nErrorPyrs++;
690  }
691  }
692 
693  reduce(nErrorPyrs, sumOp<label>());
694 
695  if (nErrorPyrs > 0)
696  {
697  if (report)
698  {
700  (
701  "polyMeshGeometry::checkFacePyramids("
702  "const bool, const scalar, const pointField&"
703  ", const labelList&, labelHashSet*)"
704  ) << "Error in face pyramids: faces pointing the wrong way!"
705  << endl;
706  }
707 
708  return true;
709  }
710  else
711  {
712  if (report)
713  {
714  Info<< "Face pyramids OK.\n" << endl;
715  }
716 
717  return false;
718  }
719 }
720 
721 
723 (
724  const bool report,
725  const scalar internalSkew,
726  const scalar boundarySkew,
727  const polyMesh& mesh,
728  const vectorField& cellCentres,
729  const vectorField& faceCentres,
730  const vectorField& faceAreas,
731  const labelList& checkFaces,
732  const List<labelPair>& baffles,
733  labelHashSet* setPtr
734 )
735 {
736  // Warn if the skew correction vector is more than skew times
737  // larger than the face area vector
738 
739  const labelList& own = mesh.faceOwner();
740  const labelList& nei = mesh.faceNeighbour();
741  const polyBoundaryMesh& patches = mesh.boundaryMesh();
742 
743  // Calculate coupled cell centre
744  vectorField neiCc(mesh.nFaces()-mesh.nInternalFaces());
745 
746  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
747  {
748  neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
749  }
750  syncTools::swapBoundaryFaceList(mesh, neiCc, true);
751 
752 
753  scalar maxSkew = 0;
754 
755  label nWarnSkew = 0;
756 
757  forAll(checkFaces, i)
758  {
759  label faceI = checkFaces[i];
760 
761  if (mesh.isInternalFace(faceI))
762  {
763  scalar skewness = calcSkewness
764  (
765  cellCentres[own[faceI]],
766  cellCentres[nei[faceI]],
767  faceCentres[faceI]
768  );
769 
770  // Check if the skewness vector is greater than the PN vector.
771  // This does not cause trouble but is a good indication of a poor
772  // mesh.
773  if (skewness > internalSkew)
774  {
775  if (report)
776  {
777  Pout<< "Severe skewness for face " << faceI
778  << " skewness = " << skewness << endl;
779  }
780 
781  if (setPtr)
782  {
783  setPtr->insert(faceI);
784  }
785 
786  nWarnSkew++;
787  }
788 
789  maxSkew = max(maxSkew, skewness);
790  }
791  else if (patches[patches.whichPatch(faceI)].coupled())
792  {
793  scalar skewness = calcSkewness
794  (
795  cellCentres[own[faceI]],
796  neiCc[faceI-mesh.nInternalFaces()],
797  faceCentres[faceI]
798  );
799 
800  // Check if the skewness vector is greater than the PN vector.
801  // This does not cause trouble but is a good indication of a poor
802  // mesh.
803  if (skewness > internalSkew)
804  {
805  if (report)
806  {
807  Pout<< "Severe skewness for coupled face " << faceI
808  << " skewness = " << skewness << endl;
809  }
810 
811  if (setPtr)
812  {
813  setPtr->insert(faceI);
814  }
815 
816  nWarnSkew++;
817  }
818 
819  maxSkew = max(maxSkew, skewness);
820  }
821  else
822  {
823  // Boundary faces: consider them to have only skewness error.
824  // (i.e. treat as if mirror cell on other side)
825 
826  vector faceNormal = faceAreas[faceI];
827  faceNormal /= mag(faceNormal) + ROOTVSMALL;
828 
829  vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]];
830 
831  vector dWall = faceNormal*(faceNormal & dOwn);
832 
833  point faceIntersection = cellCentres[own[faceI]] + dWall;
834 
835  scalar skewness =
836  mag(faceCentres[faceI] - faceIntersection)
837  /(2*mag(dWall) + ROOTVSMALL);
838 
839  // Check if the skewness vector is greater than the PN vector.
840  // This does not cause trouble but is a good indication of a poor
841  // mesh.
842  if (skewness > boundarySkew)
843  {
844  if (report)
845  {
846  Pout<< "Severe skewness for boundary face " << faceI
847  << " skewness = " << skewness << endl;
848  }
849 
850  if (setPtr)
851  {
852  setPtr->insert(faceI);
853  }
854 
855  nWarnSkew++;
856  }
857 
858  maxSkew = max(maxSkew, skewness);
859  }
860  }
861 
862  forAll(baffles, i)
863  {
864  label face0 = baffles[i].first();
865  label face1 = baffles[i].second();
866 
867  const point& ownCc = cellCentres[own[face0]];
868 
869  scalar skewness = calcSkewness
870  (
871  ownCc,
872  cellCentres[own[face1]],
873  faceCentres[face0]
874  );
875 
876  // Check if the skewness vector is greater than the PN vector.
877  // This does not cause trouble but is a good indication of a poor
878  // mesh.
879  if (skewness > internalSkew)
880  {
881  if (report)
882  {
883  Pout<< "Severe skewness for face " << face0
884  << " skewness = " << skewness << endl;
885  }
886 
887  if (setPtr)
888  {
889  setPtr->insert(face0);
890  }
891 
892  nWarnSkew++;
893  }
894 
895  maxSkew = max(maxSkew, skewness);
896  }
897 
898 
899  reduce(maxSkew, maxOp<scalar>());
900  reduce(nWarnSkew, sumOp<label>());
901 
902  if (nWarnSkew > 0)
903  {
904  if (report)
905  {
906  WarningIn
907  (
908  "polyMeshGeometry::checkFaceSkewness"
909  "(const bool, const scalar, const labelList&, labelHashSet*)"
910  ) << "Large face skewness detected. Max skewness = "
911  << 100*maxSkew
912  << " percent.\nThis may impair the quality of the result." << nl
913  << nWarnSkew << " highly skew faces detected."
914  << endl;
915  }
916 
917  return true;
918  }
919  else
920  {
921  if (report)
922  {
923  Info<< "Max skewness = " << 100*maxSkew
924  << " percent. Face skewness OK.\n" << endl;
925  }
926 
927  return false;
928  }
929 }
930 
931 
933 (
934  const bool report,
935  const scalar warnWeight,
936  const polyMesh& mesh,
937  const vectorField& cellCentres,
938  const vectorField& faceCentres,
939  const vectorField& faceAreas,
940  const labelList& checkFaces,
941  const List<labelPair>& baffles,
942  labelHashSet* setPtr
943 )
944 {
945  // Warn if the delta factor (0..1) is too large.
946 
947  const labelList& own = mesh.faceOwner();
948  const labelList& nei = mesh.faceNeighbour();
949  const polyBoundaryMesh& patches = mesh.boundaryMesh();
950 
951  // Calculate coupled cell centre
952  vectorField neiCc(mesh.nFaces()-mesh.nInternalFaces());
953 
954  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
955  {
956  neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
957  }
958  syncTools::swapBoundaryFaceList(mesh, neiCc, true);
959 
960 
961  scalar minWeight = GREAT;
962 
963  label nWarnWeight = 0;
964 
965  forAll(checkFaces, i)
966  {
967  label faceI = checkFaces[i];
968 
969  const point& fc = faceCentres[faceI];
970  const vector& fa = faceAreas[faceI];
971 
972  scalar dOwn = mag(fa & (fc-cellCentres[own[faceI]]));
973 
974  if (mesh.isInternalFace(faceI))
975  {
976  scalar dNei = mag(fa & (cellCentres[nei[faceI]]-fc));
977  scalar weight = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
978 
979  if (weight < warnWeight)
980  {
981  if (report)
982  {
983  Pout<< "Small weighting factor for face " << faceI
984  << " weight = " << weight << endl;
985  }
986 
987  if (setPtr)
988  {
989  setPtr->insert(faceI);
990  }
991 
992  nWarnWeight++;
993  }
994 
995  minWeight = min(minWeight, weight);
996  }
997  else
998  {
999  label patchI = patches.whichPatch(faceI);
1000 
1001  if (patches[patchI].coupled())
1002  {
1003  scalar dNei = mag(fa & (neiCc[faceI-mesh.nInternalFaces()]-fc));
1004  scalar weight = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1005 
1006  if (weight < warnWeight)
1007  {
1008  if (report)
1009  {
1010  Pout<< "Small weighting factor for face " << faceI
1011  << " weight = " << weight << endl;
1012  }
1013 
1014  if (setPtr)
1015  {
1016  setPtr->insert(faceI);
1017  }
1018 
1019  nWarnWeight++;
1020  }
1021 
1022  minWeight = min(minWeight, weight);
1023  }
1024  }
1025  }
1026 
1027  forAll(baffles, i)
1028  {
1029  label face0 = baffles[i].first();
1030  label face1 = baffles[i].second();
1031 
1032  const point& ownCc = cellCentres[own[face0]];
1033  const point& fc = faceCentres[face0];
1034  const vector& fa = faceAreas[face0];
1035 
1036  scalar dOwn = mag(fa & (fc-ownCc));
1037  scalar dNei = mag(fa & (cellCentres[own[face1]]-fc));
1038  scalar weight = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
1039 
1040  if (weight < warnWeight)
1041  {
1042  if (report)
1043  {
1044  Pout<< "Small weighting factor for face " << face0
1045  << " weight = " << weight << endl;
1046  }
1047 
1048  if (setPtr)
1049  {
1050  setPtr->insert(face0);
1051  }
1052 
1053  nWarnWeight++;
1054  }
1055 
1056  minWeight = min(minWeight, weight);
1057  }
1058 
1059  reduce(minWeight, minOp<scalar>());
1060  reduce(nWarnWeight, sumOp<label>());
1061 
1062  if (minWeight < warnWeight)
1063  {
1064  if (report)
1065  {
1066  WarningIn
1067  (
1068  "polyMeshGeometry::checkFaceWeights"
1069  "(const bool, const scalar, const labelList&, labelHashSet*)"
1070  ) << "Small interpolation weight detected. Min weight = "
1071  << minWeight << '.' << nl
1072  << nWarnWeight << " faces with small weights detected."
1073  << endl;
1074  }
1075 
1076  return true;
1077  }
1078  else
1079  {
1080  if (report)
1081  {
1082  Info<< "Min weight = " << minWeight
1083  << ". Weights OK.\n" << endl;
1084  }
1085 
1086  return false;
1087  }
1088 }
1089 
1090 
1093  const bool report,
1094  const scalar warnRatio,
1095  const polyMesh& mesh,
1096  const scalarField& cellVolumes,
1097  const labelList& checkFaces,
1098  const List<labelPair>& baffles,
1099  labelHashSet* setPtr
1100 )
1101 {
1102  // Warn if the volume ratio between neighbouring cells is too large
1103 
1104  const labelList& own = mesh.faceOwner();
1105  const labelList& nei = mesh.faceNeighbour();
1106  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1107 
1108  // Calculate coupled cell vol
1109  scalarField neiVols(mesh.nFaces()-mesh.nInternalFaces());
1110 
1111  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
1112  {
1113  neiVols[faceI-mesh.nInternalFaces()] = cellVolumes[own[faceI]];
1114  }
1115  syncTools::swapBoundaryFaceList(mesh, neiVols, true);
1116 
1117 
1118  scalar minRatio = GREAT;
1119 
1120  label nWarnRatio = 0;
1121 
1122  forAll(checkFaces, i)
1123  {
1124  label faceI = checkFaces[i];
1125 
1126  scalar ownVol = mag(cellVolumes[own[faceI]]);
1127 
1128  scalar neiVol = -GREAT;
1129 
1130  if (mesh.isInternalFace(faceI))
1131  {
1132  neiVol = mag(cellVolumes[nei[faceI]]);
1133  }
1134  else
1135  {
1136  label patchI = patches.whichPatch(faceI);
1137 
1138  if (patches[patchI].coupled())
1139  {
1140  neiVol = mag(neiVols[faceI-mesh.nInternalFaces()]);
1141  }
1142  }
1143 
1144  if (neiVol >= 0)
1145  {
1146  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + VSMALL);
1147 
1148  if (ratio < warnRatio)
1149  {
1150  if (report)
1151  {
1152  Pout<< "Small ratio for face " << faceI
1153  << " ratio = " << ratio << endl;
1154  }
1155 
1156  if (setPtr)
1157  {
1158  setPtr->insert(faceI);
1159  }
1160 
1161  nWarnRatio++;
1162  }
1163 
1164  minRatio = min(minRatio, ratio);
1165  }
1166  }
1167 
1168  forAll(baffles, i)
1169  {
1170  label face0 = baffles[i].first();
1171  label face1 = baffles[i].second();
1172 
1173  scalar ownVol = mag(cellVolumes[own[face0]]);
1174 
1175  scalar neiVol = mag(cellVolumes[own[face1]]);
1176 
1177  if (neiVol >= 0)
1178  {
1179  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + VSMALL);
1180 
1181  if (ratio < warnRatio)
1182  {
1183  if (report)
1184  {
1185  Pout<< "Small ratio for face " << face0
1186  << " ratio = " << ratio << endl;
1187  }
1188 
1189  if (setPtr)
1190  {
1191  setPtr->insert(face0);
1192  }
1193 
1194  nWarnRatio++;
1195  }
1196 
1197  minRatio = min(minRatio, ratio);
1198  }
1199  }
1200 
1201  reduce(minRatio, minOp<scalar>());
1202  reduce(nWarnRatio, sumOp<label>());
1203 
1204  if (minRatio < warnRatio)
1205  {
1206  if (report)
1207  {
1208  WarningIn
1209  (
1210  "polyMeshGeometry::checkVolRatio"
1211  "(const bool, const scalar, const labelList&, labelHashSet*)"
1212  ) << "Small volume ratio detected. Min ratio = "
1213  << minRatio << '.' << nl
1214  << nWarnRatio << " faces with small ratios detected."
1215  << endl;
1216  }
1217 
1218  return true;
1219  }
1220  else
1221  {
1222  if (report)
1223  {
1224  Info<< "Min ratio = " << minRatio
1225  << ". Ratios OK.\n" << endl;
1226  }
1227 
1228  return false;
1229  }
1230 }
1231 
1232 
1233 // Check convexity of angles in a face. Allow a slight non-convexity.
1234 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
1235 // (if truly concave and points not visible from face centre the face-pyramid
1236 // check in checkMesh will fail)
1239  const bool report,
1240  const scalar maxDeg,
1241  const polyMesh& mesh,
1242  const vectorField& faceAreas,
1243  const pointField& p,
1244  const labelList& checkFaces,
1245  labelHashSet* setPtr
1246 )
1247 {
1248  if (maxDeg < -SMALL || maxDeg > 180+SMALL)
1249  {
1250  FatalErrorIn
1251  (
1252  "polyMeshGeometry::checkFaceAngles"
1253  "(const bool, const scalar, const pointField&, const labelList&"
1254  ", labelHashSet*)"
1255  ) << "maxDeg should be [0..180] but is now " << maxDeg
1256  << abort(FatalError);
1257  }
1258 
1259  const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
1260 
1261  const faceList& fcs = mesh.faces();
1262 
1263  scalar maxEdgeSin = 0.0;
1264 
1265  label nConcave = 0;
1266 
1267  label errorFaceI = -1;
1268 
1269  forAll(checkFaces, i)
1270  {
1271  label faceI = checkFaces[i];
1272 
1273  const face& f = fcs[faceI];
1274 
1275  vector faceNormal = faceAreas[faceI];
1276  faceNormal /= mag(faceNormal) + VSMALL;
1277 
1278  // Get edge from f[0] to f[size-1];
1279  vector ePrev(p[f[0]] - p[f[f.size()-1]]);
1280  scalar magEPrev = mag(ePrev);
1281  ePrev /= magEPrev + VSMALL;
1282 
1283  forAll(f, fp0)
1284  {
1285  // Get vertex after fp
1286  label fp1 = f.fcIndex(fp0);
1287 
1288  // Normalized vector between two consecutive points
1289  vector e10(p[f[fp1]] - p[f[fp0]]);
1290  scalar magE10 = mag(e10);
1291  e10 /= magE10 + VSMALL;
1292 
1293  if (magEPrev > SMALL && magE10 > SMALL)
1294  {
1295  vector edgeNormal = ePrev ^ e10;
1296  scalar magEdgeNormal = mag(edgeNormal);
1297 
1298  if (magEdgeNormal < maxSin)
1299  {
1300  // Edges (almost) aligned -> face is ok.
1301  }
1302  else
1303  {
1304  // Check normal
1305  edgeNormal /= magEdgeNormal;
1306 
1307  if ((edgeNormal & faceNormal) < SMALL)
1308  {
1309  if (faceI != errorFaceI)
1310  {
1311  // Count only one error per face.
1312  errorFaceI = faceI;
1313  nConcave++;
1314  }
1315 
1316  if (setPtr)
1317  {
1318  setPtr->insert(faceI);
1319  }
1320 
1321  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
1322  }
1323  }
1324  }
1325 
1326  ePrev = e10;
1327  magEPrev = magE10;
1328  }
1329  }
1330 
1331  reduce(nConcave, sumOp<label>());
1332  reduce(maxEdgeSin, maxOp<scalar>());
1333 
1334  if (report)
1335  {
1336  if (maxEdgeSin > SMALL)
1337  {
1338  scalar maxConcaveDegr =
1339  Foam::asin(Foam::min(1.0, maxEdgeSin))
1340  * 180.0/mathematicalConstant::pi;
1341 
1342  Info<< "There are " << nConcave
1343  << " faces with concave angles between consecutive"
1344  << " edges. Max concave angle = " << maxConcaveDegr
1345  << " degrees.\n" << endl;
1346  }
1347  else
1348  {
1349  Info<< "All angles in faces are convex or less than " << maxDeg
1350  << " degrees concave.\n" << endl;
1351  }
1352  }
1353 
1354  if (nConcave > 0)
1355  {
1356  if (report)
1357  {
1358  WarningIn
1359  (
1360  "polyMeshGeometry::checkFaceAngles"
1361  "(const bool, const scalar, const pointField&"
1362  ", const labelList&, labelHashSet*)"
1363  ) << nConcave << " face points with severe concave angle (> "
1364  << maxDeg << " deg) found.\n"
1365  << endl;
1366  }
1367 
1368  return true;
1369  }
1370  else
1371  {
1372  return false;
1373  }
1374 }
1375 
1376 
1377 // Check twist of faces. Is calculated as the difference between normals of
1378 // individual triangles and the cell-cell centre edge.
1381  const bool report,
1382  const scalar minTwist,
1383  const polyMesh& mesh,
1384  const vectorField& cellCentres,
1385  const vectorField& faceAreas,
1386  const vectorField& faceCentres,
1387  const pointField& p,
1388  const labelList& checkFaces,
1389  labelHashSet* setPtr
1390 )
1391 {
1392  if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1393  {
1394  FatalErrorIn
1395  (
1396  "polyMeshGeometry::checkFaceTwist"
1397  "(const bool, const scalar, const polyMesh&, const pointField&"
1398  ", const pointField&, const pointField&, const pointField&"
1399  ", const labelList&, labelHashSet*)"
1400  ) << "minTwist should be [-1..1] but is now " << minTwist
1401  << abort(FatalError);
1402  }
1403 
1404 
1405  const faceList& fcs = mesh.faces();
1406 
1407 
1408  label nWarped = 0;
1409 
1410 // forAll(checkFaces, i)
1411 // {
1412 // label faceI = checkFaces[i];
1413 //
1414 // const face& f = fcs[faceI];
1415 //
1416 // scalar magArea = mag(faceAreas[faceI]);
1417 //
1418 // if (f.size() > 3 && magArea > VSMALL)
1419 // {
1420 // const vector nf = faceAreas[faceI] / magArea;
1421 //
1422 // const point& fc = faceCentres[faceI];
1423 //
1424 // forAll(f, fpI)
1425 // {
1426 // vector triArea
1427 // (
1428 // triPointRef
1429 // (
1430 // p[f[fpI]],
1431 // p[f.nextLabel(fpI)],
1432 // fc
1433 // ).normal()
1434 // );
1435 //
1436 // scalar magTri = mag(triArea);
1437 //
1438 // if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1439 // {
1440 // nWarped++;
1441 //
1442 // if (setPtr)
1443 // {
1444 // setPtr->insert(faceI);
1445 // }
1446 //
1447 // break;
1448 // }
1449 // }
1450 // }
1451 // }
1452 
1453 
1454  const labelList& own = mesh.faceOwner();
1455  const labelList& nei = mesh.faceNeighbour();
1456  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1457 
1458  // Calculate coupled cell centre
1459  vectorField neiCc(mesh.nFaces()-mesh.nInternalFaces());
1460 
1461  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
1462  {
1463  neiCc[faceI-mesh.nInternalFaces()] = cellCentres[own[faceI]];
1464  }
1465  syncTools::swapBoundaryFaceList(mesh, neiCc, true);
1466 
1467  forAll(checkFaces, i)
1468  {
1469  label faceI = checkFaces[i];
1470 
1471  const face& f = fcs[faceI];
1472 
1473  if (f.size() > 3)
1474  {
1475  vector nf(vector::zero);
1476 
1477  if (mesh.isInternalFace(faceI))
1478  {
1479  nf = cellCentres[nei[faceI]] - cellCentres[own[faceI]];
1480  nf /= mag(nf) + VSMALL;
1481  }
1482  else if (patches[patches.whichPatch(faceI)].coupled())
1483  {
1484  nf =
1485  neiCc[faceI-mesh.nInternalFaces()]
1486  - cellCentres[own[faceI]];
1487  nf /= mag(nf) + VSMALL;
1488  }
1489  else
1490  {
1491  nf = faceCentres[faceI] - cellCentres[own[faceI]];
1492  nf /= mag(nf) + VSMALL;
1493  }
1494 
1495  if (nf != vector::zero)
1496  {
1497  const point& fc = faceCentres[faceI];
1498 
1499  forAll(f, fpI)
1500  {
1501  vector triArea
1502  (
1503  triPointRef
1504  (
1505  p[f[fpI]],
1506  p[f.nextLabel(fpI)],
1507  fc
1508  ).normal()
1509  );
1510 
1511  scalar magTri = mag(triArea);
1512 
1513  if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
1514  {
1515  nWarped++;
1516 
1517  if (setPtr)
1518  {
1519  setPtr->insert(faceI);
1520  }
1521 
1522  break;
1523  }
1524  }
1525  }
1526  }
1527  }
1528 
1529  reduce(nWarped, sumOp<label>());
1530 
1531  if (report)
1532  {
1533  if (nWarped> 0)
1534  {
1535  Info<< "There are " << nWarped
1536  << " faces with cosine of the angle"
1537  << " between triangle normal and face normal less than "
1538  << minTwist << nl << endl;
1539  }
1540  else
1541  {
1542  Info<< "All faces are flat in that the cosine of the angle"
1543  << " between triangle normal and face normal less than "
1544  << minTwist << nl << endl;
1545  }
1546  }
1547 
1548  if (nWarped > 0)
1549  {
1550  if (report)
1551  {
1552  WarningIn
1553  (
1554  "polyMeshGeometry::checkFaceTwist"
1555  "(const bool, const scalar, const polyMesh&, const pointField&"
1556  ", const pointField&, const pointField&, const pointField&"
1557  ", const labelList&, labelHashSet*)"
1558  ) << nWarped << " faces with severe warpage "
1559  << "(cosine of the angle between triangle normal and face normal"
1560  << " < " << minTwist << ") found.\n"
1561  << endl;
1562  }
1563 
1564  return true;
1565  }
1566  else
1567  {
1568  return false;
1569  }
1570 }
1571 
1572 
1573 // Like checkFaceTwist but compares normals of consecutive triangles.
1576  const bool report,
1577  const scalar minTwist,
1578  const polyMesh& mesh,
1579  const vectorField& faceAreas,
1580  const vectorField& faceCentres,
1581  const pointField& p,
1582  const labelList& checkFaces,
1583  labelHashSet* setPtr
1584 )
1585 {
1586  if (minTwist < -1-SMALL || minTwist > 1+SMALL)
1587  {
1588  FatalErrorIn
1589  (
1590  "polyMeshGeometry::checkTriangleTwist"
1591  "(const bool, const scalar, const polyMesh&, const pointField&"
1592  ", const labelList&, labelHashSet*)"
1593  ) << "minTwist should be [-1..1] but is now " << minTwist
1594  << abort(FatalError);
1595  }
1596 
1597  const faceList& fcs = mesh.faces();
1598 
1599  label nWarped = 0;
1600 
1601  forAll(checkFaces, i)
1602  {
1603  label faceI = checkFaces[i];
1604 
1605  const face& f = fcs[faceI];
1606 
1607  if (f.size() > 3)
1608  {
1609  const point& fc = faceCentres[faceI];
1610 
1611  // Find starting triangle (at startFp) with non-zero area
1612  label startFp = -1;
1613  vector prevN;
1614 
1615  forAll(f, fp)
1616  {
1617  prevN = triPointRef
1618  (
1619  p[f[fp]],
1620  p[f.nextLabel(fp)],
1621  fc
1622  ).normal();
1623 
1624  scalar magTri = mag(prevN);
1625 
1626  if (magTri > VSMALL)
1627  {
1628  startFp = fp;
1629  prevN /= magTri;
1630  break;
1631  }
1632  }
1633 
1634  if (startFp != -1)
1635  {
1636  label fp = startFp;
1637 
1638  do
1639  {
1640  fp = f.fcIndex(fp);
1641 
1642  vector triN
1643  (
1644  triPointRef
1645  (
1646  p[f[fp]],
1647  p[f.nextLabel(fp)],
1648  fc
1649  ).normal()
1650  );
1651  scalar magTri = mag(triN);
1652 
1653  if (magTri > VSMALL)
1654  {
1655  triN /= magTri;
1656 
1657  if ((prevN & triN) < minTwist)
1658  {
1659  nWarped++;
1660 
1661  if (setPtr)
1662  {
1663  setPtr->insert(faceI);
1664  }
1665 
1666  break;
1667  }
1668 
1669  prevN = triN;
1670  }
1671  else if (minTwist > 0)
1672  {
1673  nWarped++;
1674 
1675  if (setPtr)
1676  {
1677  setPtr->insert(faceI);
1678  }
1679 
1680  break;
1681  }
1682  }
1683  while (fp != startFp);
1684  }
1685  }
1686  }
1687 
1688 
1689  reduce(nWarped, sumOp<label>());
1690 
1691  if (report)
1692  {
1693  if (nWarped> 0)
1694  {
1695  Info<< "There are " << nWarped
1696  << " faces with cosine of the angle"
1697  << " between consecutive triangle normals less than "
1698  << minTwist << nl << endl;
1699  }
1700  else
1701  {
1702  Info<< "All faces are flat in that the cosine of the angle"
1703  << " between consecutive triangle normals is less than "
1704  << minTwist << nl << endl;
1705  }
1706  }
1707 
1708  if (nWarped > 0)
1709  {
1710  if (report)
1711  {
1712  WarningIn
1713  (
1714  "polyMeshGeometry::checkTriangleTwist"
1715  "(const bool, const scalar, const polyMesh&"
1716  ", const pointField&, const labelList&, labelHashSet*)"
1717  ) << nWarped << " faces with severe warpage "
1718  << "(cosine of the angle between consecutive triangle normals"
1719  << " < " << minTwist << ") found.\n"
1720  << endl;
1721  }
1722 
1723  return true;
1724  }
1725  else
1726  {
1727  return false;
1728  }
1729 }
1730 
1731 
1734  const bool report,
1735  const scalar minArea,
1736  const polyMesh& mesh,
1737  const vectorField& faceAreas,
1738  const labelList& checkFaces,
1739  labelHashSet* setPtr
1740 )
1741 {
1742  label nZeroArea = 0;
1743 
1744  forAll(checkFaces, i)
1745  {
1746  label faceI = checkFaces[i];
1747 
1748  if (mag(faceAreas[faceI]) < minArea)
1749  {
1750  if (setPtr)
1751  {
1752  setPtr->insert(faceI);
1753  }
1754  nZeroArea++;
1755  }
1756  }
1757 
1758 
1759  reduce(nZeroArea, sumOp<label>());
1760 
1761  if (report)
1762  {
1763  if (nZeroArea > 0)
1764  {
1765  Info<< "There are " << nZeroArea
1766  << " faces with area < " << minArea << '.' << nl << endl;
1767  }
1768  else
1769  {
1770  Info<< "All faces have area > " << minArea << '.' << nl << endl;
1771  }
1772  }
1773 
1774  if (nZeroArea > 0)
1775  {
1776  if (report)
1777  {
1778  WarningIn
1779  (
1780  "polyMeshGeometry::checkFaceArea"
1781  "(const bool, const scalar, const polyMesh&"
1782  ", const pointField&, const labelList&, labelHashSet*)"
1783  ) << nZeroArea << " faces with area < " << minArea
1784  << " found.\n"
1785  << endl;
1786  }
1787 
1788  return true;
1789  }
1790  else
1791  {
1792  return false;
1793  }
1794 }
1795 
1796 
1799  const bool report,
1800  const scalar warnDet,
1801  const polyMesh& mesh,
1802  const vectorField& faceAreas,
1803  const labelList& checkFaces,
1804  const labelList& affectedCells,
1805  labelHashSet* setPtr
1806 )
1807 {
1808  const cellList& cells = mesh.cells();
1809 
1810  scalar minDet = GREAT;
1811  scalar sumDet = 0.0;
1812  label nSumDet = 0;
1813  label nWarnDet = 0;
1814 
1815  forAll(affectedCells, i)
1816  {
1817  const cell& cFaces = cells[affectedCells[i]];
1818 
1819  tensor areaSum(tensor::zero);
1820  scalar magAreaSum = 0;
1821 
1822  forAll(cFaces, cFaceI)
1823  {
1824  label faceI = cFaces[cFaceI];
1825 
1826  scalar magArea = mag(faceAreas[faceI]);
1827 
1828  magAreaSum += magArea;
1829  areaSum += faceAreas[faceI]*(faceAreas[faceI]/(magArea+VSMALL));
1830  }
1831 
1832  scalar scaledDet = det(areaSum/(magAreaSum+VSMALL))/0.037037037037037;
1833 
1834  minDet = min(minDet, scaledDet);
1835  sumDet += scaledDet;
1836  nSumDet++;
1837 
1838  if (scaledDet < warnDet)
1839  {
1840  if (setPtr)
1841  {
1842  // Insert all faces of the cell.
1843  forAll(cFaces, cFaceI)
1844  {
1845  label faceI = cFaces[cFaceI];
1846  setPtr->insert(faceI);
1847  }
1848  }
1849  nWarnDet++;
1850  }
1851  }
1852 
1853  reduce(minDet, minOp<scalar>());
1854  reduce(sumDet, sumOp<scalar>());
1855  reduce(nSumDet, sumOp<label>());
1856  reduce(nWarnDet, sumOp<label>());
1857 
1858  if (report)
1859  {
1860  if (nSumDet > 0)
1861  {
1862  Info<< "Cell determinant (1 = uniform cube) : average = "
1863  << sumDet / nSumDet << " min = " << minDet << endl;
1864  }
1865 
1866  if (nWarnDet > 0)
1867  {
1868  Info<< "There are " << nWarnDet
1869  << " cells with determinant < " << warnDet << '.' << nl
1870  << endl;
1871  }
1872  else
1873  {
1874  Info<< "All faces have determinant > " << warnDet << '.' << nl
1875  << endl;
1876  }
1877  }
1878 
1879  if (nWarnDet > 0)
1880  {
1881  if (report)
1882  {
1883  WarningIn
1884  (
1885  "polyMeshGeometry::checkCellDeterminant"
1886  "(const bool, const scalar, const polyMesh&"
1887  ", const pointField&, const labelList&, const labelList&"
1888  ", labelHashSet*)"
1889  ) << nWarnDet << " cells with determinant < " << warnDet
1890  << " found.\n"
1891  << endl;
1892  }
1893 
1894  return true;
1895  }
1896  else
1897  {
1898  return false;
1899  }
1900 }
1901 
1902 
1905  const bool report,
1906  const scalar orthWarn,
1907  const labelList& checkFaces,
1908  const List<labelPair>& baffles,
1909  labelHashSet* setPtr
1910 ) const
1911 {
1912  return checkFaceDotProduct
1913  (
1914  report,
1915  orthWarn,
1916  mesh_,
1917  cellCentres_,
1918  faceAreas_,
1919  checkFaces,
1920  baffles,
1921  setPtr
1922  );
1923 }
1924 
1925 
1928  const bool report,
1929  const scalar minPyrVol,
1930  const pointField& p,
1931  const labelList& checkFaces,
1932  const List<labelPair>& baffles,
1933  labelHashSet* setPtr
1934 ) const
1935 {
1936  return checkFacePyramids
1937  (
1938  report,
1939  minPyrVol,
1940  mesh_,
1941  cellCentres_,
1942  p,
1943  checkFaces,
1944  baffles,
1945  setPtr
1946  );
1947 }
1948 
1949 
1952  const bool report,
1953  const scalar internalSkew,
1954  const scalar boundarySkew,
1955  const labelList& checkFaces,
1956  const List<labelPair>& baffles,
1957  labelHashSet* setPtr
1958 ) const
1959 {
1960  return checkFaceSkewness
1961  (
1962  report,
1963  internalSkew,
1964  boundarySkew,
1965  mesh_,
1966  cellCentres_,
1967  faceCentres_,
1968  faceAreas_,
1969  checkFaces,
1970  baffles,
1971  setPtr
1972  );
1973 }
1974 
1975 
1978  const bool report,
1979  const scalar warnWeight,
1980  const labelList& checkFaces,
1981  const List<labelPair>& baffles,
1982  labelHashSet* setPtr
1983 ) const
1984 {
1985  return checkFaceWeights
1986  (
1987  report,
1988  warnWeight,
1989  mesh_,
1990  cellCentres_,
1991  faceCentres_,
1992  faceAreas_,
1993  checkFaces,
1994  baffles,
1995  setPtr
1996  );
1997 }
1998 
1999 
2002  const bool report,
2003  const scalar warnRatio,
2004  const labelList& checkFaces,
2005  const List<labelPair>& baffles,
2006  labelHashSet* setPtr
2007 ) const
2008 {
2009  return checkVolRatio
2010  (
2011  report,
2012  warnRatio,
2013  mesh_,
2014  cellVolumes_,
2015  checkFaces,
2016  baffles,
2017  setPtr
2018  );
2019 }
2020 
2021 
2024  const bool report,
2025  const scalar maxDeg,
2026  const pointField& p,
2027  const labelList& checkFaces,
2028  labelHashSet* setPtr
2029 ) const
2030 {
2031  return checkFaceAngles
2032  (
2033  report,
2034  maxDeg,
2035  mesh_,
2036  faceAreas_,
2037  p,
2038  checkFaces,
2039  setPtr
2040  );
2041 }
2042 
2043 
2046  const bool report,
2047  const scalar minTwist,
2048  const pointField& p,
2049  const labelList& checkFaces,
2050  labelHashSet* setPtr
2051 ) const
2052 {
2053  return checkFaceTwist
2054  (
2055  report,
2056  minTwist,
2057  mesh_,
2058  cellCentres_,
2059  faceAreas_,
2060  faceCentres_,
2061  p,
2062  checkFaces,
2063  setPtr
2064  );
2065 }
2066 
2067 
2070  const bool report,
2071  const scalar minTwist,
2072  const pointField& p,
2073  const labelList& checkFaces,
2074  labelHashSet* setPtr
2075 ) const
2076 {
2077  return checkTriangleTwist
2078  (
2079  report,
2080  minTwist,
2081  mesh_,
2082  faceAreas_,
2083  faceCentres_,
2084  p,
2085  checkFaces,
2086  setPtr
2087  );
2088 }
2089 
2090 
2093  const bool report,
2094  const scalar minArea,
2095  const labelList& checkFaces,
2096  labelHashSet* setPtr
2097 ) const
2098 {
2099  return checkFaceArea
2100  (
2101  report,
2102  minArea,
2103  mesh_,
2104  faceAreas_,
2105  checkFaces,
2106  setPtr
2107  );
2108 }
2109 
2110 
2113  const bool report,
2114  const scalar warnDet,
2115  const labelList& checkFaces,
2116  const labelList& affectedCells,
2117  labelHashSet* setPtr
2118 ) const
2119 {
2120  return checkCellDeterminant
2121  (
2122  report,
2123  warnDet,
2124  mesh_,
2125  faceAreas_,
2126  checkFaces,
2127  affectedCells,
2128  setPtr
2129  );
2130 }
2131 
2132 
2133 // ************************ vim: set sw=4 sts=4 et: ************************ //