polyMeshGeometry.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "polyMeshGeometry.H"
28 #include "pyramidPointFaceRef.H"
29 #include "tetPointRef.H"
30 #include "syncTools.H"
31 #include "unitConversion.H"
32 #include "primitiveMeshTools.H"
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(polyMeshGeometry, 0);
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 = Zero;
67  scalar sumA = 0.0;
68  vector sumAc = 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) = 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) = 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  << radToDeg(::acos(dDotS))
250  << " deg." << endl;
251  }
252 
253  severeNonOrth++;
254  }
255  else
256  {
257  // Non-orthogonality greater than 90 deg
258  if (report)
259  {
261  << "Severe non-orthogonality detected for face "
262  << facei
263  << " between cells " << mesh.faceOwner()[facei]
264  << " and " << nei
265  << ": Angle = "
266  << radToDeg(::acos(dDotS))
267  << " deg." << endl;
268  }
269 
270  errorNonOrth++;
271  }
272 
273  if (setPtr)
274  {
275  setPtr->insert(facei);
276  }
277  }
278  return dDotS;
279 }
280 
281 
282 // Create the neighbour pyramid - it will have positive volume
283 bool Foam::polyMeshGeometry::checkFaceTet
284 (
285  const polyMesh& mesh,
286  const bool report,
287  const scalar minTetQuality,
288  const pointField& p,
289  const label facei,
290  const point& fc, // face centre
291  const point& cc, // cell centre
292 
293  labelHashSet* setPtr
294 )
295 {
296  const face& f = mesh.faces()[facei];
297 
298  forAll(f, fp)
299  {
300  scalar tetQual = tetPointRef
301  (
302  p[f[fp]],
303  p[f.nextLabel(fp)],
304  fc,
305  cc
306  ).quality();
307 
308  if (tetQual < minTetQuality)
309  {
310  if (report)
311  {
312  Pout<< "bool polyMeshGeometry::checkFaceTets("
313  << "const bool, const scalar, const pointField&"
314  << ", const pointField&"
315  << ", const labelList&, labelHashSet*) : "
316  << "face " << facei
317  << " has a triangle that points the wrong way."
318  << endl
319  << "Tet quality: " << tetQual
320  << " Face " << facei
321  << endl;
322  }
323 
324  if (setPtr)
325  {
326  setPtr->insert(facei);
327  }
328  return true;
329  }
330  }
331  return false;
332 }
333 
334 
335 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
336 
337 // Construct from components
339 :
340  mesh_(mesh)
341 {
342  correct();
343 }
344 
345 
346 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
347 
349 {
350  faceAreas_ = mesh_.faceAreas();
351  faceCentres_ = mesh_.faceCentres();
352  cellCentres_ = mesh_.cellCentres();
353  cellVolumes_ = mesh_.cellVolumes();
354 }
355 
356 
358 (
359  const pointField& p,
360  const labelList& changedFaces
361 )
362 {
363  // Update face quantities
364  updateFaceCentresAndAreas(p, changedFaces);
365  // Update cell quantities from face quantities
366  updateCellCentresAndVols(affectedCells(mesh_, changedFaces), changedFaces);
367 }
368 
369 
371 (
372  const bool report,
373  const scalar orthWarn,
374  const polyMesh& mesh,
375  const vectorField& cellCentres,
376  const vectorField& faceAreas,
377  const labelList& checkFaces,
378  const List<labelPair>& baffles,
379  labelHashSet* setPtr
380 )
381 {
382  // for all internal and coupled faces check theat the d dot S product
383  // is positive
384 
385  const labelList& own = mesh.faceOwner();
386  const labelList& nei = mesh.faceNeighbour();
387  const polyBoundaryMesh& patches = mesh.boundaryMesh();
388 
389  // Severe nonorthogonality threshold
390  const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
391 
392  // Calculate coupled cell centre
393  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
394 
395  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
396  {
397  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
398  }
399 
401 
402  scalar minDDotS = great;
403 
404  scalar sumDDotS = 0;
405  label nDDotS = 0;
406 
407  label severeNonOrth = 0;
408 
409  label errorNonOrth = 0;
410 
411  forAll(checkFaces, i)
412  {
413  label facei = checkFaces[i];
414 
415  const point& ownCc = cellCentres[own[facei]];
416 
417  if (mesh.isInternalFace(facei))
418  {
419  scalar dDotS = checkNonOrtho
420  (
421  mesh,
422  report,
423  severeNonorthogonalityThreshold,
424  facei,
425  faceAreas[facei],
426  cellCentres[nei[facei]] - ownCc,
427 
428  severeNonOrth,
429  errorNonOrth,
430  setPtr
431  );
432 
433  if (dDotS < minDDotS)
434  {
435  minDDotS = dDotS;
436  }
437 
438  sumDDotS += dDotS;
439  nDDotS++;
440  }
441  else
442  {
443  label patchi = patches.whichPatch(facei);
444 
445  if (patches[patchi].coupled())
446  {
447  scalar dDotS = checkNonOrtho
448  (
449  mesh,
450  report,
451  severeNonorthogonalityThreshold,
452  facei,
453  faceAreas[facei],
454  neiCc[facei-mesh.nInternalFaces()] - ownCc,
455 
456  severeNonOrth,
457  errorNonOrth,
458  setPtr
459  );
460 
461  if (dDotS < minDDotS)
462  {
463  minDDotS = dDotS;
464  }
465 
466  sumDDotS += dDotS;
467  nDDotS++;
468  }
469  }
470  }
471 
472  forAll(baffles, i)
473  {
474  label face0 = baffles[i].first();
475  label face1 = baffles[i].second();
476 
477  const point& ownCc = cellCentres[own[face0]];
478 
479  scalar dDotS = checkNonOrtho
480  (
481  mesh,
482  report,
483  severeNonorthogonalityThreshold,
484  face0,
485  faceAreas[face0],
486  cellCentres[own[face1]] - ownCc,
487 
488  severeNonOrth,
489  errorNonOrth,
490  setPtr
491  );
492 
493  if (dDotS < minDDotS)
494  {
495  minDDotS = dDotS;
496  }
497 
498  sumDDotS += dDotS;
499  nDDotS++;
500  }
501 
502  reduce(minDDotS, minOp<scalar>());
503  reduce(sumDDotS, sumOp<scalar>());
504  reduce(nDDotS, sumOp<label>());
505  reduce(severeNonOrth, sumOp<label>());
506  reduce(errorNonOrth, sumOp<label>());
507 
508  // Only report if there are some internal faces
509  if (nDDotS > 0)
510  {
511  if (report && minDDotS < severeNonorthogonalityThreshold)
512  {
513  Info<< "Number of non-orthogonality errors: " << errorNonOrth
514  << ". Number of severely non-orthogonal faces: "
515  << severeNonOrth << "." << endl;
516  }
517  }
518 
519  if (report)
520  {
521  if (nDDotS > 0)
522  {
523  Info<< "Mesh non-orthogonality Max: "
524  << radToDeg(::acos(minDDotS))
525  << " average: " << radToDeg(::acos(sumDDotS/nDDotS))
526  << endl;
527  }
528  }
529 
530  if (errorNonOrth > 0)
531  {
532  if (report)
533  {
535  << "Error in non-orthogonality detected" << endl;
536  }
537 
538  return true;
539  }
540  else
541  {
542  if (report)
543  {
544  Info<< "Non-orthogonality check OK.\n" << endl;
545  }
546 
547  return false;
548  }
549 }
550 
551 
553 (
554  const bool report,
555  const scalar minPyrVol,
556  const polyMesh& mesh,
557  const vectorField& cellCentres,
558  const pointField& p,
559  const labelList& checkFaces,
560  const List<labelPair>& baffles,
561  labelHashSet* setPtr
562 )
563 {
564  // check whether face area vector points to the cell with higher label
565  const labelList& own = mesh.faceOwner();
566  const labelList& nei = mesh.faceNeighbour();
567 
568  const faceList& f = mesh.faces();
569 
570  label nErrorPyrs = 0;
571 
572  forAll(checkFaces, i)
573  {
574  label facei = checkFaces[i];
575 
576  // Create the owner pyramid - it will have negative volume
577  scalar pyrVol = pyramidPointFaceRef
578  (
579  f[facei],
580  cellCentres[own[facei]]
581  ).mag(p);
582 
583  if (pyrVol > -minPyrVol)
584  {
585  if (report)
586  {
587  Pout<< "bool polyMeshGeometry::checkFacePyramids("
588  << "const bool, const scalar, const pointField&"
589  << ", const labelList&, labelHashSet*): "
590  << "face " << facei << " points the wrong way. " << endl
591  << "Pyramid volume: " << -pyrVol
592  << " Face " << f[facei] << " area: " << f[facei].mag(p)
593  << " Owner cell: " << own[facei] << endl
594  << "Owner cell vertex labels: "
595  << mesh.cells()[own[facei]].labels(f)
596  << endl;
597  }
598 
599 
600  if (setPtr)
601  {
602  setPtr->insert(facei);
603  }
604 
605  nErrorPyrs++;
606  }
607 
608  if (mesh.isInternalFace(facei))
609  {
610  // Create the neighbour pyramid - it will have positive volume
611  scalar pyrVol =
612  pyramidPointFaceRef(f[facei], cellCentres[nei[facei]]).mag(p);
613 
614  if (pyrVol < minPyrVol)
615  {
616  if (report)
617  {
618  Pout<< "bool polyMeshGeometry::checkFacePyramids("
619  << "const bool, const scalar, const pointField&"
620  << ", const labelList&, labelHashSet*): "
621  << "face " << facei << " points the wrong way. " << endl
622  << "Pyramid volume: " << -pyrVol
623  << " Face " << f[facei] << " area: " << f[facei].mag(p)
624  << " Neighbour cell: " << nei[facei] << endl
625  << "Neighbour cell vertex labels: "
626  << mesh.cells()[nei[facei]].labels(f)
627  << endl;
628  }
629 
630  if (setPtr)
631  {
632  setPtr->insert(facei);
633  }
634 
635  nErrorPyrs++;
636  }
637  }
638  }
639 
640  forAll(baffles, i)
641  {
642  label face0 = baffles[i].first();
643  label face1 = baffles[i].second();
644 
645  const point& ownCc = cellCentres[own[face0]];
646 
647  // Create the owner pyramid - it will have negative volume
648  scalar pyrVolOwn = pyramidPointFaceRef
649  (
650  f[face0],
651  ownCc
652  ).mag(p);
653 
654  if (pyrVolOwn > -minPyrVol)
655  {
656  if (report)
657  {
658  Pout<< "bool polyMeshGeometry::checkFacePyramids("
659  << "const bool, const scalar, const pointField&"
660  << ", const labelList&, labelHashSet*): "
661  << "face " << face0 << " points the wrong way. " << endl
662  << "Pyramid volume: " << -pyrVolOwn
663  << " Face " << f[face0] << " area: " << f[face0].mag(p)
664  << " Owner cell: " << own[face0] << endl
665  << "Owner cell vertex labels: "
666  << mesh.cells()[own[face0]].labels(f)
667  << endl;
668  }
669 
670 
671  if (setPtr)
672  {
673  setPtr->insert(face0);
674  }
675 
676  nErrorPyrs++;
677  }
678 
679  // Create the neighbour pyramid - it will have positive volume
680  scalar pyrVolNbr =
681  pyramidPointFaceRef(f[face0], cellCentres[own[face1]]).mag(p);
682 
683  if (pyrVolNbr < minPyrVol)
684  {
685  if (report)
686  {
687  Pout<< "bool polyMeshGeometry::checkFacePyramids("
688  << "const bool, const scalar, const pointField&"
689  << ", const labelList&, labelHashSet*): "
690  << "face " << face0 << " points the wrong way. " << endl
691  << "Pyramid volume: " << -pyrVolNbr
692  << " Face " << f[face0] << " area: " << f[face0].mag(p)
693  << " Neighbour cell: " << own[face1] << endl
694  << "Neighbour cell vertex labels: "
695  << mesh.cells()[own[face1]].labels(f)
696  << endl;
697  }
698 
699  if (setPtr)
700  {
701  setPtr->insert(face0);
702  }
703 
704  nErrorPyrs++;
705  }
706  }
707 
708  reduce(nErrorPyrs, sumOp<label>());
709 
710  if (nErrorPyrs > 0)
711  {
712  if (report)
713  {
715  << "Error in face pyramids: faces pointing the wrong way."
716  << endl;
717  }
718 
719  return true;
720  }
721  else
722  {
723  if (report)
724  {
725  Info<< "Face pyramids OK.\n" << endl;
726  }
727 
728  return false;
729  }
730 }
731 
732 
734 (
735  const bool report,
736  const scalar minTetQuality,
737  const polyMesh& mesh,
738  const vectorField& cellCentres,
739  const vectorField& faceCentres,
740  const pointField& p,
741  const labelList& checkFaces,
742  const List<labelPair>& baffles,
743  labelHashSet* setPtr
744 )
745 {
746  // check whether decomposing each cell into tets results in
747  // positive volume, non-flat tets
748  const labelList& own = mesh.faceOwner();
749  const labelList& nei = mesh.faceNeighbour();
750  const polyBoundaryMesh& patches = mesh.boundaryMesh();
751 
752  // Calculate coupled cell centre
753  pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
754 
755  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
756  {
757  neiCc[facei - mesh.nInternalFaces()] = cellCentres[own[facei]];
758  }
759 
761 
762  label nErrorTets = 0;
763 
764  forAll(checkFaces, i)
765  {
766  label facei = checkFaces[i];
767 
768  // Create the owner pyramid - note: exchange cell and face centre
769  // to get positive volume.
770  bool tetError = checkFaceTet
771  (
772  mesh,
773  report,
774  minTetQuality,
775  p,
776  facei,
777  cellCentres[own[facei]], // face centre
778  faceCentres[facei], // cell centre
779  setPtr
780  );
781 
782  if (tetError)
783  {
784  nErrorTets++;
785  }
786 
787  if (mesh.isInternalFace(facei))
788  {
789  // Create the neighbour tets - they will have positive volume
790  bool tetError = checkFaceTet
791  (
792  mesh,
793  report,
794  minTetQuality,
795  p,
796  facei,
797  faceCentres[facei], // face centre
798  cellCentres[nei[facei]], // cell centre
799  setPtr
800  );
801 
802  if (tetError)
803  {
804  nErrorTets++;
805  }
806 
807  if
808  (
810  (
811  mesh,
812  facei,
813  minTetQuality,
814  report
815  ) == -1
816  )
817  {
818  if (setPtr)
819  {
820  setPtr->insert(facei);
821  }
822 
823  nErrorTets++;
824  }
825  }
826  else
827  {
828  label patchi = patches.whichPatch(facei);
829 
830  if (patches[patchi].coupled())
831  {
832  if
833  (
835  (
836  mesh,
837  facei,
838  neiCc[facei - mesh.nInternalFaces()],
839  minTetQuality,
840  report
841  ) == -1
842  )
843  {
844  if (setPtr)
845  {
846  setPtr->insert(facei);
847  }
848 
849  nErrorTets++;
850  }
851  }
852  else
853  {
854  if
855  (
857  (
858  mesh,
859  facei,
860  minTetQuality,
861  report
862  ) == -1
863  )
864  {
865  if (setPtr)
866  {
867  setPtr->insert(facei);
868  }
869 
870  nErrorTets++;
871  }
872  }
873  }
874  }
875 
876  forAll(baffles, i)
877  {
878  label face0 = baffles[i].first();
879  label face1 = baffles[i].second();
880 
881  bool tetError = checkFaceTet
882  (
883  mesh,
884  report,
885  minTetQuality,
886  p,
887  face0,
888  cellCentres[own[face0]], // face centre
889  faceCentres[face0], // cell centre
890  setPtr
891  );
892 
893  if (tetError)
894  {
895  nErrorTets++;
896  }
897 
898  // Create the neighbour tets - they will have positive volume
899  tetError = checkFaceTet
900  (
901  mesh,
902  report,
903  minTetQuality,
904  p,
905  face0,
906  faceCentres[face0], // face centre
907  cellCentres[own[face1]], // cell centre
908  setPtr
909  );
910 
911  if (tetError)
912  {
913  nErrorTets++;
914  }
915 
916  if
917  (
919  (
920  mesh,
921  face0,
922  cellCentres[own[face1]],
923  minTetQuality,
924  report
925  ) == -1
926  )
927  {
928  if (setPtr)
929  {
930  setPtr->insert(face0);
931  }
932 
933  nErrorTets++;
934  }
935  }
936 
937  reduce(nErrorTets, sumOp<label>());
938 
939  if (nErrorTets > 0)
940  {
941  if (report)
942  {
944  << "Error in face decomposition: negative tets."
945  << endl;
946  }
947 
948  return true;
949  }
950  else
951  {
952  if (report)
953  {
954  Info<< "Face tets OK.\n" << endl;
955  }
956 
957  return false;
958  }
959 }
960 
961 
963 (
964  const bool report,
965  const scalar internalSkew,
966  const scalar boundarySkew,
967  const polyMesh& mesh,
968  const pointField& points,
969  const vectorField& cellCentres,
970  const vectorField& faceCentres,
971  const vectorField& faceAreas,
972  const labelList& checkFaces,
973  const List<labelPair>& baffles,
974  labelHashSet* setPtr
975 )
976 {
977  // Warn if the skew correction vector is more than skew times
978  // larger than the face area vector
979 
980  const labelList& own = mesh.faceOwner();
981  const labelList& nei = mesh.faceNeighbour();
982  const polyBoundaryMesh& patches = mesh.boundaryMesh();
983 
984  // Calculate coupled cell centre
985  pointField neiCc;
986  syncTools::swapBoundaryCellPositions(mesh, cellCentres, neiCc);
987 
988  scalar maxSkew = 0;
989 
990  label nWarnSkew = 0;
991 
992  forAll(checkFaces, i)
993  {
994  label facei = checkFaces[i];
995 
996  if (mesh.isInternalFace(facei))
997  {
998  scalar skewness = primitiveMeshTools::faceSkewness
999  (
1000  mesh,
1001  points,
1002  faceCentres,
1003  faceAreas,
1004 
1005  facei,
1006  cellCentres[own[facei]],
1007  cellCentres[nei[facei]]
1008  );
1009 
1010  // Check if the skewness vector is greater than the PN vector.
1011  // This does not cause trouble but is a good indication of a poor
1012  // mesh.
1013  if (skewness > internalSkew)
1014  {
1015  if (report)
1016  {
1017  Pout<< "Severe skewness for face " << facei
1018  << " skewness = " << skewness << endl;
1019  }
1020 
1021  if (setPtr)
1022  {
1023  setPtr->insert(facei);
1024  }
1025 
1026  nWarnSkew++;
1027  }
1028 
1029  maxSkew = max(maxSkew, skewness);
1030  }
1031  else if (patches[patches.whichPatch(facei)].coupled())
1032  {
1033  scalar skewness = primitiveMeshTools::faceSkewness
1034  (
1035  mesh,
1036  points,
1037  faceCentres,
1038  faceAreas,
1039 
1040  facei,
1041  cellCentres[own[facei]],
1042  neiCc[facei-mesh.nInternalFaces()]
1043  );
1044 
1045  // Check if the skewness vector is greater than the PN vector.
1046  // This does not cause trouble but is a good indication of a poor
1047  // mesh.
1048  if (skewness > internalSkew)
1049  {
1050  if (report)
1051  {
1052  Pout<< "Severe skewness for coupled face " << facei
1053  << " skewness = " << skewness << endl;
1054  }
1055 
1056  if (setPtr)
1057  {
1058  setPtr->insert(facei);
1059  }
1060 
1061  nWarnSkew++;
1062  }
1063 
1064  maxSkew = max(maxSkew, skewness);
1065  }
1066  else
1067  {
1069  (
1070  mesh,
1071  points,
1072  faceCentres,
1073  faceAreas,
1074 
1075  facei,
1076  cellCentres[own[facei]]
1077  );
1078 
1079 
1080  // Check if the skewness vector is greater than the PN vector.
1081  // This does not cause trouble but is a good indication of a poor
1082  // mesh.
1083  if (skewness > boundarySkew)
1084  {
1085  if (report)
1086  {
1087  Pout<< "Severe skewness for boundary face " << facei
1088  << " skewness = " << skewness << endl;
1089  }
1090 
1091  if (setPtr)
1092  {
1093  setPtr->insert(facei);
1094  }
1095 
1096  nWarnSkew++;
1097  }
1098 
1099  maxSkew = max(maxSkew, skewness);
1100  }
1101  }
1102 
1103  forAll(baffles, i)
1104  {
1105  label face0 = baffles[i].first();
1106  label face1 = baffles[i].second();
1107 
1108  const point& ownCc = cellCentres[own[face0]];
1109  const point& neiCc = cellCentres[own[face1]];
1110 
1111  scalar skewness = primitiveMeshTools::faceSkewness
1112  (
1113  mesh,
1114  points,
1115  faceCentres,
1116  faceAreas,
1117 
1118  face0,
1119  ownCc,
1120  neiCc
1121  );
1122 
1123  // Check if the skewness vector is greater than the PN vector.
1124  // This does not cause trouble but is a good indication of a poor
1125  // mesh.
1126  if (skewness > internalSkew)
1127  {
1128  if (report)
1129  {
1130  Pout<< "Severe skewness for face " << face0
1131  << " skewness = " << skewness << endl;
1132  }
1133 
1134  if (setPtr)
1135  {
1136  setPtr->insert(face0);
1137  }
1138 
1139  nWarnSkew++;
1140  }
1141 
1142  maxSkew = max(maxSkew, skewness);
1143  }
1144 
1145 
1146  reduce(maxSkew, maxOp<scalar>());
1147  reduce(nWarnSkew, sumOp<label>());
1148 
1149  if (nWarnSkew > 0)
1150  {
1151  if (report)
1152  {
1154  << 100*maxSkew
1155  << " percent.\nThis may impair the quality of the result." << nl
1156  << nWarnSkew << " highly skew faces detected."
1157  << endl;
1158  }
1159 
1160  return true;
1161  }
1162  else
1163  {
1164  if (report)
1165  {
1166  Info<< "Max skewness = " << 100*maxSkew
1167  << " percent. Face skewness OK.\n" << endl;
1168  }
1169 
1170  return false;
1171  }
1172 }
1173 
1174 
1177  const bool report,
1178  const scalar warnWeight,
1179  const polyMesh& mesh,
1180  const vectorField& cellCentres,
1181  const vectorField& faceCentres,
1182  const vectorField& faceAreas,
1183  const labelList& checkFaces,
1184  const List<labelPair>& baffles,
1185  labelHashSet* setPtr
1186 )
1187 {
1188  // Warn if the delta factor (0..1) is too large.
1189 
1190  const labelList& own = mesh.faceOwner();
1191  const labelList& nei = mesh.faceNeighbour();
1192  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1193 
1194  // Calculate coupled cell centre
1195  pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
1196 
1197  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
1198  {
1199  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
1200  }
1202 
1203 
1204  scalar minWeight = great;
1205 
1206  label nWarnWeight = 0;
1207 
1208  forAll(checkFaces, i)
1209  {
1210  label facei = checkFaces[i];
1211 
1212  const point& fc = faceCentres[facei];
1213  const vector& fa = faceAreas[facei];
1214 
1215  scalar dOwn = mag(fa & (fc-cellCentres[own[facei]]));
1216 
1217  if (mesh.isInternalFace(facei))
1218  {
1219  scalar dNei = mag(fa & (cellCentres[nei[facei]]-fc));
1220  scalar weight = min(dNei,dOwn)/(dNei+dOwn+vSmall);
1221 
1222  if (weight < warnWeight)
1223  {
1224  if (report)
1225  {
1226  Pout<< "Small weighting factor for face " << facei
1227  << " weight = " << weight << endl;
1228  }
1229 
1230  if (setPtr)
1231  {
1232  setPtr->insert(facei);
1233  }
1234 
1235  nWarnWeight++;
1236  }
1237 
1238  minWeight = min(minWeight, weight);
1239  }
1240  else
1241  {
1242  label patchi = patches.whichPatch(facei);
1243 
1244  if (patches[patchi].coupled())
1245  {
1246  scalar dNei = mag(fa & (neiCc[facei-mesh.nInternalFaces()]-fc));
1247  scalar weight = min(dNei,dOwn)/(dNei+dOwn+vSmall);
1248 
1249  if (weight < warnWeight)
1250  {
1251  if (report)
1252  {
1253  Pout<< "Small weighting factor for face " << facei
1254  << " weight = " << weight << endl;
1255  }
1256 
1257  if (setPtr)
1258  {
1259  setPtr->insert(facei);
1260  }
1261 
1262  nWarnWeight++;
1263  }
1264 
1265  minWeight = min(minWeight, weight);
1266  }
1267  }
1268  }
1269 
1270  forAll(baffles, i)
1271  {
1272  label face0 = baffles[i].first();
1273  label face1 = baffles[i].second();
1274 
1275  const point& ownCc = cellCentres[own[face0]];
1276  const point& fc = faceCentres[face0];
1277  const vector& fa = faceAreas[face0];
1278 
1279  scalar dOwn = mag(fa & (fc-ownCc));
1280  scalar dNei = mag(fa & (cellCentres[own[face1]]-fc));
1281  scalar weight = min(dNei,dOwn)/(dNei+dOwn+vSmall);
1282 
1283  if (weight < warnWeight)
1284  {
1285  if (report)
1286  {
1287  Pout<< "Small weighting factor for face " << face0
1288  << " weight = " << weight << endl;
1289  }
1290 
1291  if (setPtr)
1292  {
1293  setPtr->insert(face0);
1294  }
1295 
1296  nWarnWeight++;
1297  }
1298 
1299  minWeight = min(minWeight, weight);
1300  }
1301 
1302  reduce(minWeight, minOp<scalar>());
1303  reduce(nWarnWeight, sumOp<label>());
1304 
1305  if (minWeight < warnWeight)
1306  {
1307  if (report)
1308  {
1310  << minWeight << '.' << nl
1311  << nWarnWeight << " faces with small weights detected."
1312  << endl;
1313  }
1314 
1315  return true;
1316  }
1317  else
1318  {
1319  if (report)
1320  {
1321  Info<< "Min weight = " << minWeight
1322  << ". Weights OK.\n" << endl;
1323  }
1324 
1325  return false;
1326  }
1327 }
1328 
1329 
1332  const bool report,
1333  const scalar warnRatio,
1334  const polyMesh& mesh,
1335  const scalarField& cellVolumes,
1336  const labelList& checkFaces,
1337  const List<labelPair>& baffles,
1338  labelHashSet* setPtr
1339 )
1340 {
1341  // Warn if the volume ratio between neighbouring cells is too large
1342 
1343  const labelList& own = mesh.faceOwner();
1344  const labelList& nei = mesh.faceNeighbour();
1345  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1346 
1347  // Calculate coupled cell vol
1348  scalarField neiVols(mesh.nFaces()-mesh.nInternalFaces());
1349 
1350  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
1351  {
1352  neiVols[facei-mesh.nInternalFaces()] = cellVolumes[own[facei]];
1353  }
1354  syncTools::swapBoundaryFaceList(mesh, neiVols);
1355 
1356 
1357  scalar minRatio = great;
1358 
1359  label nWarnRatio = 0;
1360 
1361  forAll(checkFaces, i)
1362  {
1363  label facei = checkFaces[i];
1364 
1365  scalar ownVol = mag(cellVolumes[own[facei]]);
1366 
1367  scalar neiVol = -great;
1368 
1369  if (mesh.isInternalFace(facei))
1370  {
1371  neiVol = mag(cellVolumes[nei[facei]]);
1372  }
1373  else
1374  {
1375  label patchi = patches.whichPatch(facei);
1376 
1377  if (patches[patchi].coupled())
1378  {
1379  neiVol = mag(neiVols[facei-mesh.nInternalFaces()]);
1380  }
1381  }
1382 
1383  if (neiVol >= 0)
1384  {
1385  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + vSmall);
1386 
1387  if (ratio < warnRatio)
1388  {
1389  if (report)
1390  {
1391  Pout<< "Small ratio for face " << facei
1392  << " ratio = " << ratio << endl;
1393  }
1394 
1395  if (setPtr)
1396  {
1397  setPtr->insert(facei);
1398  }
1399 
1400  nWarnRatio++;
1401  }
1402 
1403  minRatio = min(minRatio, ratio);
1404  }
1405  }
1406 
1407  forAll(baffles, i)
1408  {
1409  label face0 = baffles[i].first();
1410  label face1 = baffles[i].second();
1411 
1412  scalar ownVol = mag(cellVolumes[own[face0]]);
1413 
1414  scalar neiVol = mag(cellVolumes[own[face1]]);
1415 
1416  if (neiVol >= 0)
1417  {
1418  scalar ratio = min(ownVol, neiVol) / (max(ownVol, neiVol) + vSmall);
1419 
1420  if (ratio < warnRatio)
1421  {
1422  if (report)
1423  {
1424  Pout<< "Small ratio for face " << face0
1425  << " ratio = " << ratio << endl;
1426  }
1427 
1428  if (setPtr)
1429  {
1430  setPtr->insert(face0);
1431  }
1432 
1433  nWarnRatio++;
1434  }
1435 
1436  minRatio = min(minRatio, ratio);
1437  }
1438  }
1439 
1440  reduce(minRatio, minOp<scalar>());
1441  reduce(nWarnRatio, sumOp<label>());
1442 
1443  if (minRatio < warnRatio)
1444  {
1445  if (report)
1446  {
1448  << minRatio << '.' << nl
1449  << nWarnRatio << " faces with small ratios detected."
1450  << endl;
1451  }
1452 
1453  return true;
1454  }
1455  else
1456  {
1457  if (report)
1458  {
1459  Info<< "Min ratio = " << minRatio
1460  << ". Ratios OK.\n" << endl;
1461  }
1462 
1463  return false;
1464  }
1465 }
1466 
1467 
1468 // Check convexity of angles in a face. Allow a slight non-convexity.
1469 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
1470 // (if truly concave and points not visible from face centre the face-pyramid
1471 // check in checkMesh will fail)
1474  const bool report,
1475  const scalar maxDeg,
1476  const polyMesh& mesh,
1477  const vectorField& faceAreas,
1478  const pointField& p,
1479  const labelList& checkFaces,
1480  labelHashSet* setPtr
1481 )
1482 {
1483  if (maxDeg < -small || maxDeg > 180+small)
1484  {
1486  << "maxDeg should be [0..180] but is now " << maxDeg
1487  << abort(FatalError);
1488  }
1489 
1490  const scalar maxSin = Foam::sin(degToRad(maxDeg));
1491 
1492  const faceList& fcs = mesh.faces();
1493 
1494  scalar maxEdgeSin = 0.0;
1495 
1496  label nConcave = 0;
1497 
1498  label errorFacei = -1;
1499 
1500  forAll(checkFaces, i)
1501  {
1502  label facei = checkFaces[i];
1503 
1504  const face& f = fcs[facei];
1505 
1506  vector faceNormal = faceAreas[facei];
1507  faceNormal /= mag(faceNormal) + vSmall;
1508 
1509  // Get edge from f[0] to f[size-1];
1510  vector ePrev(p[f.first()] - p[f.last()]);
1511  scalar magEPrev = mag(ePrev);
1512  ePrev /= magEPrev + vSmall;
1513 
1514  forAll(f, fp0)
1515  {
1516  // Get vertex after fp
1517  label fp1 = f.fcIndex(fp0);
1518 
1519  // Normalized vector between two consecutive points
1520  vector e10(p[f[fp1]] - p[f[fp0]]);
1521  scalar magE10 = mag(e10);
1522  e10 /= magE10 + vSmall;
1523 
1524  if (magEPrev > small && magE10 > small)
1525  {
1526  vector edgeNormal = ePrev ^ e10;
1527  scalar magEdgeNormal = mag(edgeNormal);
1528 
1529  if (magEdgeNormal < maxSin)
1530  {
1531  // Edges (almost) aligned -> face is ok.
1532  }
1533  else
1534  {
1535  // Check normal
1536  edgeNormal /= magEdgeNormal;
1537 
1538  if ((edgeNormal & faceNormal) < small)
1539  {
1540  if (facei != errorFacei)
1541  {
1542  // Count only one error per face.
1543  errorFacei = facei;
1544  nConcave++;
1545  }
1546 
1547  if (setPtr)
1548  {
1549  setPtr->insert(facei);
1550  }
1551 
1552  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
1553  }
1554  }
1555  }
1556 
1557  ePrev = e10;
1558  magEPrev = magE10;
1559  }
1560  }
1561 
1562  reduce(nConcave, sumOp<label>());
1563  reduce(maxEdgeSin, maxOp<scalar>());
1564 
1565  if (report)
1566  {
1567  if (maxEdgeSin > small)
1568  {
1569  scalar maxConcaveDegr =
1570  radToDeg(Foam::asin(Foam::min(1.0, maxEdgeSin)));
1571 
1572  Info<< "There are " << nConcave
1573  << " faces with concave angles between consecutive"
1574  << " edges. Max concave angle = " << maxConcaveDegr
1575  << " degrees.\n" << endl;
1576  }
1577  else
1578  {
1579  Info<< "All angles in faces are convex or less than " << maxDeg
1580  << " degrees concave.\n" << endl;
1581  }
1582  }
1583 
1584  if (nConcave > 0)
1585  {
1586  if (report)
1587  {
1589  << nConcave << " face points with severe concave angle (> "
1590  << maxDeg << " deg) found.\n"
1591  << endl;
1592  }
1593 
1594  return true;
1595  }
1596  else
1597  {
1598  return false;
1599  }
1600 }
1601 
1602 
1603 // Check twist of faces. Is calculated as the difference between normals of
1604 // individual triangles and the cell-cell centre edge.
1607  const bool report,
1608  const scalar minTwist,
1609  const polyMesh& mesh,
1610  const vectorField& cellCentres,
1611  const vectorField& faceAreas,
1612  const vectorField& faceCentres,
1613  const pointField& p,
1614  const labelList& checkFaces,
1615  labelHashSet* setPtr
1616 )
1617 {
1618  if (minTwist < -1-small || minTwist > 1+small)
1619  {
1621  << "minTwist should be [-1..1] but is now " << minTwist
1622  << abort(FatalError);
1623  }
1624 
1625 
1626  const faceList& fcs = mesh.faces();
1627 
1628 
1629  label nWarped = 0;
1630 
1631 // forAll(checkFaces, i)
1632 // {
1633 // label facei = checkFaces[i];
1634 //
1635 // const face& f = fcs[facei];
1636 //
1637 // scalar magArea = mag(faceAreas[facei]);
1638 //
1639 // if (f.size() > 3 && magArea > vSmall)
1640 // {
1641 // const vector nf = faceAreas[facei] / magArea;
1642 //
1643 // const point& fc = faceCentres[facei];
1644 //
1645 // forAll(f, fpI)
1646 // {
1647 // vector triArea
1648 // (
1649 // triPointRef
1650 // (
1651 // p[f[fpI]],
1652 // p[f.nextLabel(fpI)],
1653 // fc
1654 // ).area()
1655 // );
1656 //
1657 // scalar magTri = mag(triArea);
1658 //
1659 // if (magTri > vSmall && ((nf & triArea/magTri) < minTwist))
1660 // {
1661 // nWarped++;
1662 //
1663 // if (setPtr)
1664 // {
1665 // setPtr->insert(facei);
1666 // }
1667 //
1668 // break;
1669 // }
1670 // }
1671 // }
1672 // }
1673 
1674 
1675  const labelList& own = mesh.faceOwner();
1676  const labelList& nei = mesh.faceNeighbour();
1677  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1678 
1679  // Calculate coupled cell centre
1680  pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
1681 
1682  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
1683  {
1684  neiCc[facei-mesh.nInternalFaces()] = cellCentres[own[facei]];
1685  }
1687 
1688  forAll(checkFaces, i)
1689  {
1690  label facei = checkFaces[i];
1691 
1692  const face& f = fcs[facei];
1693 
1694  if (f.size() > 3)
1695  {
1696  vector nf(Zero);
1697 
1698  if (mesh.isInternalFace(facei))
1699  {
1700  nf = cellCentres[nei[facei]] - cellCentres[own[facei]];
1701  nf /= mag(nf) + vSmall;
1702  }
1703  else if (patches[patches.whichPatch(facei)].coupled())
1704  {
1705  nf =
1706  neiCc[facei-mesh.nInternalFaces()]
1707  - cellCentres[own[facei]];
1708  nf /= mag(nf) + vSmall;
1709  }
1710  else
1711  {
1712  nf = faceCentres[facei] - cellCentres[own[facei]];
1713  nf /= mag(nf) + vSmall;
1714  }
1715 
1716  if (nf != vector::zero)
1717  {
1718  const point& fc = faceCentres[facei];
1719 
1720  forAll(f, fpI)
1721  {
1722  vector triArea
1723  (
1724  triPointRef
1725  (
1726  p[f[fpI]],
1727  p[f.nextLabel(fpI)],
1728  fc
1729  ).area()
1730  );
1731 
1732  scalar magTri = mag(triArea);
1733 
1734  if (magTri > vSmall && ((nf & triArea/magTri) < minTwist))
1735  {
1736  nWarped++;
1737 
1738  if (setPtr)
1739  {
1740  setPtr->insert(facei);
1741  }
1742 
1743  break;
1744  }
1745  }
1746  }
1747  }
1748  }
1749 
1750  reduce(nWarped, sumOp<label>());
1751 
1752  if (report)
1753  {
1754  if (nWarped> 0)
1755  {
1756  Info<< "There are " << nWarped
1757  << " faces with cosine of the angle"
1758  << " between triangle normal and face normal less than "
1759  << minTwist << nl << endl;
1760  }
1761  else
1762  {
1763  Info<< "All faces are flat in that the cosine of the angle"
1764  << " between triangle normal and face normal less than "
1765  << minTwist << nl << endl;
1766  }
1767  }
1768 
1769  if (nWarped > 0)
1770  {
1771  if (report)
1772  {
1774  << nWarped << " faces with severe warpage "
1775  << "(cosine of the angle between triangle normal and "
1776  << "face normal < " << minTwist << ") found.\n"
1777  << endl;
1778  }
1779 
1780  return true;
1781  }
1782  else
1783  {
1784  return false;
1785  }
1786 }
1787 
1788 
1789 // Like checkFaceTwist but compares normals of consecutive triangles.
1792  const bool report,
1793  const scalar minTwist,
1794  const polyMesh& mesh,
1795  const vectorField& faceAreas,
1796  const vectorField& faceCentres,
1797  const pointField& p,
1798  const labelList& checkFaces,
1799  labelHashSet* setPtr
1800 )
1801 {
1802  if (minTwist < -1-small || minTwist > 1+small)
1803  {
1805  << "minTwist should be [-1..1] but is now " << minTwist
1806  << abort(FatalError);
1807  }
1808 
1809  const faceList& fcs = mesh.faces();
1810 
1811  label nWarped = 0;
1812 
1813  forAll(checkFaces, i)
1814  {
1815  label facei = checkFaces[i];
1816 
1817  const face& f = fcs[facei];
1818 
1819  if (f.size() > 3)
1820  {
1821  const point& fc = faceCentres[facei];
1822 
1823  // Find starting triangle (at startFp) with non-zero area
1824  label startFp = -1;
1825  vector prevN;
1826 
1827  forAll(f, fp)
1828  {
1829  prevN = triPointRef
1830  (
1831  p[f[fp]],
1832  p[f.nextLabel(fp)],
1833  fc
1834  ).area();
1835 
1836  scalar magTri = mag(prevN);
1837 
1838  if (magTri > vSmall)
1839  {
1840  startFp = fp;
1841  prevN /= magTri;
1842  break;
1843  }
1844  }
1845 
1846  if (startFp != -1)
1847  {
1848  label fp = startFp;
1849 
1850  do
1851  {
1852  fp = f.fcIndex(fp);
1853 
1854  vector triN
1855  (
1856  triPointRef
1857  (
1858  p[f[fp]],
1859  p[f.nextLabel(fp)],
1860  fc
1861  ).area()
1862  );
1863  scalar magTri = mag(triN);
1864 
1865  if (magTri > vSmall)
1866  {
1867  triN /= magTri;
1868 
1869  if ((prevN & triN) < minTwist)
1870  {
1871  nWarped++;
1872 
1873  if (setPtr)
1874  {
1875  setPtr->insert(facei);
1876  }
1877 
1878  break;
1879  }
1880 
1881  prevN = triN;
1882  }
1883  else if (minTwist > 0)
1884  {
1885  nWarped++;
1886 
1887  if (setPtr)
1888  {
1889  setPtr->insert(facei);
1890  }
1891 
1892  break;
1893  }
1894  }
1895  while (fp != startFp);
1896  }
1897  }
1898  }
1899 
1900 
1901  reduce(nWarped, sumOp<label>());
1902 
1903  if (report)
1904  {
1905  if (nWarped> 0)
1906  {
1907  Info<< "There are " << nWarped
1908  << " faces with cosine of the angle"
1909  << " between consecutive triangle normals less than "
1910  << minTwist << nl << endl;
1911  }
1912  else
1913  {
1914  Info<< "All faces are flat in that the cosine of the angle"
1915  << " between consecutive triangle normals is less than "
1916  << minTwist << nl << endl;
1917  }
1918  }
1919 
1920  if (nWarped > 0)
1921  {
1922  if (report)
1923  {
1925  << nWarped << " faces with severe warpage "
1926  << "(cosine of the angle between consecutive triangle normals"
1927  << " < " << minTwist << ") found.\n"
1928  << endl;
1929  }
1930 
1931  return true;
1932  }
1933  else
1934  {
1935  return false;
1936  }
1937 }
1938 
1939 
1942  const bool report,
1943  const scalar minFlatness,
1944  const polyMesh& mesh,
1945  const vectorField& faceAreas,
1946  const vectorField& faceCentres,
1947  const pointField& p,
1948  const labelList& checkFaces,
1949  labelHashSet* setPtr
1950 )
1951 {
1952  if (minFlatness < -small || minFlatness > 1+small)
1953  {
1955  << "minFlatness should be [0..1] but is now " << minFlatness
1956  << abort(FatalError);
1957  }
1958 
1959  const faceList& fcs = mesh.faces();
1960 
1961  label nWarped = 0;
1962 
1963  forAll(checkFaces, i)
1964  {
1965  label facei = checkFaces[i];
1966 
1967  const face& f = fcs[facei];
1968 
1969  if (f.size() > 3)
1970  {
1971  const point& fc = faceCentres[facei];
1972 
1973  // Sum triangle areas
1974  scalar sumArea = 0.0;
1975 
1976  forAll(f, fp)
1977  {
1978  sumArea += triPointRef
1979  (
1980  p[f[fp]],
1981  p[f.nextLabel(fp)],
1982  fc
1983  ).mag();
1984  }
1985 
1986  if (sumArea/mag(faceAreas[facei]) < minFlatness)
1987  {
1988  nWarped++;
1989 
1990  if (setPtr)
1991  {
1992  setPtr->insert(facei);
1993  }
1994  }
1995  }
1996  }
1997 
1998  reduce(nWarped, sumOp<label>());
1999 
2000  if (report)
2001  {
2002  if (nWarped> 0)
2003  {
2004  Info<< "There are " << nWarped
2005  << " faces with area of individual triangles"
2006  << " compared to overall area less than "
2007  << minFlatness << nl << endl;
2008  }
2009  else
2010  {
2011  Info<< "All faces are flat in that the area of individual triangles"
2012  << " compared to overall area is less than "
2013  << minFlatness << nl << endl;
2014  }
2015  }
2016 
2017  if (nWarped > 0)
2018  {
2019  if (report)
2020  {
2022  << nWarped << " non-flat faces "
2023  << "(area of individual triangles"
2024  << " compared to overall area"
2025  << " < " << minFlatness << ") found.\n"
2026  << endl;
2027  }
2028 
2029  return true;
2030  }
2031  else
2032  {
2033  return false;
2034  }
2035 }
2036 
2037 
2040  const bool report,
2041  const scalar minArea,
2042  const polyMesh& mesh,
2043  const vectorField& faceAreas,
2044  const labelList& checkFaces,
2045  labelHashSet* setPtr
2046 )
2047 {
2048  label nZeroArea = 0;
2049 
2050  forAll(checkFaces, i)
2051  {
2052  label facei = checkFaces[i];
2053 
2054  if (mag(faceAreas[facei]) < minArea)
2055  {
2056  if (setPtr)
2057  {
2058  setPtr->insert(facei);
2059  }
2060  nZeroArea++;
2061  }
2062  }
2063 
2064 
2065  reduce(nZeroArea, sumOp<label>());
2066 
2067  if (report)
2068  {
2069  if (nZeroArea > 0)
2070  {
2071  Info<< "There are " << nZeroArea
2072  << " faces with area < " << minArea << '.' << nl << endl;
2073  }
2074  else
2075  {
2076  Info<< "All faces have area > " << minArea << '.' << nl << endl;
2077  }
2078  }
2079 
2080  if (nZeroArea > 0)
2081  {
2082  if (report)
2083  {
2085  << nZeroArea << " faces with area < " << minArea
2086  << " found.\n"
2087  << endl;
2088  }
2089 
2090  return true;
2091  }
2092  else
2093  {
2094  return false;
2095  }
2096 }
2097 
2098 
2101  const bool report,
2102  const scalar warnDet,
2103  const polyMesh& mesh,
2104  const vectorField& faceAreas,
2105  const labelList& checkFaces,
2106  const labelList& affectedCells,
2107  labelHashSet* setPtr
2108 )
2109 {
2110  const cellList& cells = mesh.cells();
2111 
2112  scalar minDet = great;
2113  scalar sumDet = 0.0;
2114  label nSumDet = 0;
2115  label nWarnDet = 0;
2116 
2117  forAll(affectedCells, i)
2118  {
2119  const cell& cFaces = cells[affectedCells[i]];
2120 
2121  tensor areaSum(Zero);
2122  scalar magAreaSum = 0;
2123 
2124  forAll(cFaces, cFacei)
2125  {
2126  label facei = cFaces[cFacei];
2127 
2128  scalar magArea = mag(faceAreas[facei]);
2129 
2130  magAreaSum += magArea;
2131  areaSum += faceAreas[facei]*(faceAreas[facei]/(magArea+vSmall));
2132  }
2133 
2134  scalar scaledDet = det(areaSum/(magAreaSum+vSmall))/0.037037037037037;
2135 
2136  minDet = min(minDet, scaledDet);
2137  sumDet += scaledDet;
2138  nSumDet++;
2139 
2140  if (scaledDet < warnDet)
2141  {
2142  if (setPtr)
2143  {
2144  // Insert all faces of the cell.
2145  forAll(cFaces, cFacei)
2146  {
2147  label facei = cFaces[cFacei];
2148  setPtr->insert(facei);
2149  }
2150  }
2151  nWarnDet++;
2152  }
2153  }
2154 
2155  reduce(minDet, minOp<scalar>());
2156  reduce(sumDet, sumOp<scalar>());
2157  reduce(nSumDet, sumOp<label>());
2158  reduce(nWarnDet, sumOp<label>());
2159 
2160  if (report)
2161  {
2162  if (nSumDet > 0)
2163  {
2164  Info<< "Cell determinant (1 = uniform cube) : average = "
2165  << sumDet / nSumDet << " min = " << minDet << endl;
2166  }
2167 
2168  if (nWarnDet > 0)
2169  {
2170  Info<< "There are " << nWarnDet
2171  << " cells with determinant < " << warnDet << '.' << nl
2172  << endl;
2173  }
2174  else
2175  {
2176  Info<< "All faces have determinant > " << warnDet << '.' << nl
2177  << endl;
2178  }
2179  }
2180 
2181  if (nWarnDet > 0)
2182  {
2183  if (report)
2184  {
2186  << nWarnDet << " cells with determinant < " << warnDet
2187  << " found.\n"
2188  << endl;
2189  }
2190 
2191  return true;
2192  }
2193  else
2194  {
2195  return false;
2196  }
2197 }
2198 
2199 
2202  const bool report,
2203  const scalar orthWarn,
2204  const labelList& checkFaces,
2205  const List<labelPair>& baffles,
2206  labelHashSet* setPtr
2207 ) const
2208 {
2209  return checkFaceDotProduct
2210  (
2211  report,
2212  orthWarn,
2213  mesh_,
2214  cellCentres_,
2215  faceAreas_,
2216  checkFaces,
2217  baffles,
2218  setPtr
2219  );
2220 }
2221 
2222 
2225  const bool report,
2226  const scalar minPyrVol,
2227  const pointField& p,
2228  const labelList& checkFaces,
2229  const List<labelPair>& baffles,
2230  labelHashSet* setPtr
2231 ) const
2232 {
2233  return checkFacePyramids
2234  (
2235  report,
2236  minPyrVol,
2237  mesh_,
2238  cellCentres_,
2239  p,
2240  checkFaces,
2241  baffles,
2242  setPtr
2243  );
2244 }
2245 
2246 
2249  const bool report,
2250  const scalar minTetQuality,
2251  const pointField& p,
2252  const labelList& checkFaces,
2253  const List<labelPair>& baffles,
2254  labelHashSet* setPtr
2255 ) const
2256 {
2257  return checkFaceTets
2258  (
2259  report,
2260  minTetQuality,
2261  mesh_,
2262  cellCentres_,
2263  faceCentres_,
2264  p,
2265  checkFaces,
2266  baffles,
2267  setPtr
2268  );
2269 }
2270 
2271 
2272 //bool Foam::polyMeshGeometry::checkFaceSkewness
2273 //(
2274 // const bool report,
2275 // const scalar internalSkew,
2276 // const scalar boundarySkew,
2277 // const labelList& checkFaces,
2278 // const List<labelPair>& baffles,
2279 // labelHashSet* setPtr
2280 //) const
2281 //{
2282 // return checkFaceSkewness
2283 // (
2284 // report,
2285 // internalSkew,
2286 // boundarySkew,
2287 // mesh_,
2288 // cellCentres_,
2289 // faceCentres_,
2290 // faceAreas_,
2291 // checkFaces,
2292 // baffles,
2293 // setPtr
2294 // );
2295 //}
2296 
2297 
2300  const bool report,
2301  const scalar warnWeight,
2302  const labelList& checkFaces,
2303  const List<labelPair>& baffles,
2304  labelHashSet* setPtr
2305 ) const
2306 {
2307  return checkFaceWeights
2308  (
2309  report,
2310  warnWeight,
2311  mesh_,
2312  cellCentres_,
2313  faceCentres_,
2314  faceAreas_,
2315  checkFaces,
2316  baffles,
2317  setPtr
2318  );
2319 }
2320 
2321 
2324  const bool report,
2325  const scalar warnRatio,
2326  const labelList& checkFaces,
2327  const List<labelPair>& baffles,
2328  labelHashSet* setPtr
2329 ) const
2330 {
2331  return checkVolRatio
2332  (
2333  report,
2334  warnRatio,
2335  mesh_,
2336  cellVolumes_,
2337  checkFaces,
2338  baffles,
2339  setPtr
2340  );
2341 }
2342 
2343 
2346  const bool report,
2347  const scalar maxDeg,
2348  const pointField& p,
2349  const labelList& checkFaces,
2350  labelHashSet* setPtr
2351 ) const
2352 {
2353  return checkFaceAngles
2354  (
2355  report,
2356  maxDeg,
2357  mesh_,
2358  faceAreas_,
2359  p,
2360  checkFaces,
2361  setPtr
2362  );
2363 }
2364 
2365 
2368  const bool report,
2369  const scalar minTwist,
2370  const pointField& p,
2371  const labelList& checkFaces,
2372  labelHashSet* setPtr
2373 ) const
2374 {
2375  return checkFaceTwist
2376  (
2377  report,
2378  minTwist,
2379  mesh_,
2380  cellCentres_,
2381  faceAreas_,
2382  faceCentres_,
2383  p,
2384  checkFaces,
2385  setPtr
2386  );
2387 }
2388 
2389 
2392  const bool report,
2393  const scalar minTwist,
2394  const pointField& p,
2395  const labelList& checkFaces,
2396  labelHashSet* setPtr
2397 ) const
2398 {
2399  return checkTriangleTwist
2400  (
2401  report,
2402  minTwist,
2403  mesh_,
2404  faceAreas_,
2405  faceCentres_,
2406  p,
2407  checkFaces,
2408  setPtr
2409  );
2410 }
2411 
2412 
2415  const bool report,
2416  const scalar minFlatness,
2417  const pointField& p,
2418  const labelList& checkFaces,
2419  labelHashSet* setPtr
2420 ) const
2421 {
2422  return checkFaceFlatness
2423  (
2424  report,
2425  minFlatness,
2426  mesh_,
2427  faceAreas_,
2428  faceCentres_,
2429  p,
2430  checkFaces,
2431  setPtr
2432  );
2433 }
2434 
2435 
2438  const bool report,
2439  const scalar minArea,
2440  const labelList& checkFaces,
2441  labelHashSet* setPtr
2442 ) const
2443 {
2444  return checkFaceArea
2445  (
2446  report,
2447  minArea,
2448  mesh_,
2449  faceAreas_,
2450  checkFaces,
2451  setPtr
2452  );
2453 }
2454 
2455 
2458  const bool report,
2459  const scalar warnDet,
2460  const labelList& checkFaces,
2461  const labelList& affectedCells,
2462  labelHashSet* setPtr
2463 ) const
2464 {
2465  return checkCellDeterminant
2466  (
2467  report,
2468  warnDet,
2469  mesh_,
2470  faceAreas_,
2471  checkFaces,
2472  affectedCells,
2473  setPtr
2474  );
2475 }
2476 
2477 
2478 // ************************************************************************* //
static bool checkFacePyramids(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
See primitiveMesh.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
static bool checkFaceTets(const bool report, const scalar minPyrVol, const polyMesh &, const vectorField &cellCentres, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *)
See primitiveMesh.
dimensionedScalar acos(const dimensionedScalar &ds)
tetrahedron< point, const point & > tetPointRef
Definition: tetPointRef.H:44
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A triangle primitive used to calculate face areas and swept volumes.
Definition: triangle.H:57
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.
static bool checkTriangleTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Consecutive triangle (from face-centre decomposition) normals.
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
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
static void swapBoundaryFacePositions(const polyMesh &mesh, UList< point > &l)
Swap coupled positions.
Definition: syncTools.H:446
label nFaces() const
Unit conversion functions.
static bool checkFaceTwist(const bool report, const scalar minTwist, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Triangle (from face-centre decomposition) normal v.s.
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:256
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)
patches[0]
static void swapBoundaryCellPositions(const polyMesh &mesh, const UList< point > &cellData, List< point > &neighbourCellData)
Swap to obtain neighbour cell positions for all boundary faces.
Definition: syncTools.C:31
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:108
T & first()
Return the first element of the list.
Definition: UListI.H:114
static bool checkVolRatio(const bool report, const scalar warnRatio, const polyMesh &mesh, const scalarField &cellVolumes, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Cell volume ratio of neighbouring cells (1 for regular mesh)
dimensionedScalar asin(const dimensionedScalar &ds)
static bool checkFaceFlatness(const bool report, const scalar minFlatness, const polyMesh &, const vectorField &faceAreas, const vectorField &faceCentres, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
Area of faces v.s. sum of triangle areas.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
static label findBasePoint(const polyMesh &mesh, label fI, scalar tol, bool report=false)
Find the base point to use for a minimum triangle.
#define SeriousErrorInFunction
Report an error message using Foam::SeriousError.
static tmp< scalarField > faceSkewness(const primitiveMesh &mesh, const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate skewness field.
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
static scalar boundaryFaceSkewness(const primitiveMesh &mesh, const pointField &p, const vectorField &fCtrs, const vectorField &fAreas, const label facei, const point &ownCc)
Skewness of single boundary face.
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
label nPoints
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static label findSharedBasePoint(const polyMesh &mesh, label fI, const point &nCc, scalar tol, bool report=false)
Find the first suitable base point to use for a minimum.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
const vectorField & cellCentres() const
pyramid< point, const point &, const face & > pyramidPointFaceRef
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
dimensionedScalar sin(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:265
static bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Small faces.
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)
const vectorField & faceCentres() const
static bool checkFaceAngles(const bool report, const scalar maxDeg, const polyMesh &mesh, const vectorField &faceAreas, const pointField &p, const labelList &checkFaces, labelHashSet *setPtr)
See primitiveMesh.
label patchi
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
static bool checkFaceDotProduct(const bool report, const scalar orthWarn, const polyMesh &, const vectorField &cellCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
const vectorField & faceAreas() const
const dimensionedScalar c
Speed of light in a vacuum.
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
void correct()
Take over properties from mesh.
static labelList affectedCells(const polyMesh &, const labelList &changedFaces)
Helper function: get affected cells from faces.
messageStream Info
static bool checkFaceWeights(const bool report, const scalar warnWeight, const polyMesh &mesh, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
Interpolation weights (0.5 for regular mesh)
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
static bool checkFaceSkewness(const bool report, const scalar internalSkew, const scalar boundarySkew, const polyMesh &mesh, const pointField &points, const vectorField &cellCentres, const vectorField &faceCentres, const vectorField &faceAreas, const labelList &checkFaces, const List< labelPair > &baffles, labelHashSet *setPtr)
See primitiveMesh.
T & last()
Return the last element of the list.
Definition: UListI.H:128
static bool checkCellDeterminant(const bool report, const scalar minDet, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
Area of internal faces v.s. boundary faces.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Namespace for OpenFOAM.
polyMeshGeometry(const polyMesh &)
Construct from mesh.
const scalarField & cellVolumes() const