FreeFOAM The Cross-Platform CFD Toolkit
primitiveMeshCheck.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 <OpenFOAM/primitiveMesh.H>
28 #include <OpenFOAM/ListOps.H>
30 #include <OpenFOAM/SortableList.H>
31 
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 Foam::scalar Foam::primitiveMesh::closedThreshold_ = 1.0e-6;
36 Foam::scalar Foam::primitiveMesh::aspectThreshold_ = 1000;
37 Foam::scalar Foam::primitiveMesh::nonOrthThreshold_ = 70; // deg
38 Foam::scalar Foam::primitiveMesh::skewThreshold_ = 4;
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
44 {
45  if (debug)
46  {
47  Info<< "bool primitiveMesh::checkClosedBoundary("
48  << "const bool) const: "
49  << "checking whether the boundary is closed" << endl;
50  }
51 
52  // Loop through all boundary faces and sum up the face area vectors.
53  // For a closed boundary, this should be zero in all vector components
54 
55  vector sumClosed(vector::zero);
56  scalar sumMagClosedBoundary = 0;
57 
58  const vectorField& areas = faceAreas();
59 
60  for (label faceI = nInternalFaces(); faceI < areas.size(); faceI++)
61  {
62  sumClosed += areas[faceI];
63  sumMagClosedBoundary += mag(areas[faceI]);
64  }
65 
66  reduce(sumClosed, sumOp<vector>());
67  reduce(sumMagClosedBoundary, sumOp<scalar>());
68 
69  vector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
70 
71  if (cmptMax(cmptMag(openness)) > closedThreshold_)
72  {
73  if (debug || report)
74  {
75  Info<< " ***Boundary openness " << openness
76  << " possible hole in boundary description."
77  << endl;
78  }
79 
80  return true;
81  }
82  else
83  {
84  if (debug || report)
85  {
86  Info<< " Boundary openness " << openness << " OK."
87  << endl;
88  }
89 
90  return false;
91  }
92 }
93 
94 
96 (
97  const bool report,
98  labelHashSet* setPtr,
99  labelHashSet* aspectSetPtr,
100  const Vector<label>& meshD
101 ) const
102 {
103  if (debug)
104  {
105  Info<< "bool primitiveMesh::checkClosedCells("
106  << "const bool, labelHashSet*, labelHashSet*"
107  << ", const Vector<label>&) const: "
108  << "checking whether cells are closed" << endl;
109  }
110 
111  // Check that all cells labels are valid
112  const cellList& c = cells();
113 
114  label nErrorClosed = 0;
115 
116  forAll (c, cI)
117  {
118  const cell& curCell = c[cI];
119 
120  if (min(curCell) < 0 || max(curCell) > nFaces())
121  {
122  if (setPtr)
123  {
124  setPtr->insert(cI);
125  }
126 
127  nErrorClosed++;
128  }
129  }
130 
131  if (nErrorClosed > 0)
132  {
133  if (debug || report)
134  {
135  Info<< " ***Cells with invalid face labels found, number of cells "
136  << nErrorClosed << endl;
137  }
138 
139  return true;
140  }
141 
142  // Loop through cell faces and sum up the face area vectors for each cell.
143  // This should be zero in all vector components
144 
145  vectorField sumClosed(nCells(), vector::zero);
146  vectorField sumMagClosed(nCells(), vector::zero);
147 
148  const labelList& own = faceOwner();
149  const labelList& nei = faceNeighbour();
150  const vectorField& areas = faceAreas();
151 
152  forAll (own, faceI)
153  {
154  // Add to owner
155  sumClosed[own[faceI]] += areas[faceI];
156  sumMagClosed[own[faceI]] += cmptMag(areas[faceI]);
157  }
158 
159  forAll (nei, faceI)
160  {
161  // Subtract from neighbour
162  sumClosed[nei[faceI]] -= areas[faceI];
163  sumMagClosed[nei[faceI]] += cmptMag(areas[faceI]);
164  }
165 
166  label nOpen = 0;
167  scalar maxOpennessCell = 0;
168 
169  label nAspect = 0;
170  scalar maxAspectRatio = 0;
171 
172  const scalarField& vols = cellVolumes();
173 
174  label nDims = 0;
175  for (direction dir = 0; dir < vector::nComponents; dir++)
176  {
177  if (meshD[dir] == 1)
178  {
179  nDims++;
180  }
181  }
182 
183 
184  // Check the sums
185  forAll(sumClosed, cellI)
186  {
187  scalar maxOpenness = 0;
188 
189  for(direction cmpt=0; cmpt<vector::nComponents; cmpt++)
190  {
191  maxOpenness = max
192  (
193  maxOpenness,
194  mag(sumClosed[cellI][cmpt])
195  /(sumMagClosed[cellI][cmpt] + VSMALL)
196  );
197  }
198 
199  maxOpennessCell = max(maxOpennessCell, maxOpenness);
200 
201  if (maxOpenness > closedThreshold_)
202  {
203  if (setPtr)
204  {
205  setPtr->insert(cellI);
206  }
207 
208  nOpen++;
209  }
210 
211  // Calculate the aspect ration as the maximum of Cartesian component
212  // aspect ratio to the total area hydraulic area aspect ratio
213  scalar minCmpt = VGREAT;
214  scalar maxCmpt = -VGREAT;
215  for (direction dir = 0; dir < vector::nComponents; dir++)
216  {
217  if (meshD[dir] == 1)
218  {
219  minCmpt = min(minCmpt, sumMagClosed[cellI][dir]);
220  maxCmpt = max(maxCmpt, sumMagClosed[cellI][dir]);
221  }
222  }
223 
224  scalar aspectRatio = maxCmpt/(minCmpt + VSMALL);
225  if (nDims == 3)
226  {
227  aspectRatio = max
228  (
229  aspectRatio,
230  1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
231  );
232  }
233 
234  maxAspectRatio = max(maxAspectRatio, aspectRatio);
235 
236  if (aspectRatio > aspectThreshold_)
237  {
238  if (aspectSetPtr)
239  {
240  aspectSetPtr->insert(cellI);
241  }
242 
243  nAspect++;
244  }
245  }
246 
247  reduce(nOpen, sumOp<label>());
248  reduce(maxOpennessCell, maxOp<scalar>());
249 
250  reduce(nAspect, sumOp<label>());
251  reduce(maxAspectRatio, maxOp<scalar>());
252 
253 
254  if (nOpen > 0)
255  {
256  if (debug || report)
257  {
258  Info<< " ***Open cells found, max cell openness: "
259  << maxOpennessCell << ", number of open cells " << nOpen
260  << endl;
261  }
262 
263  return true;
264  }
265 
266  if (nAspect > 0)
267  {
268  if (debug || report)
269  {
270  Info<< " ***High aspect ratio cells found, Max aspect ratio: "
271  << maxAspectRatio
272  << ", number of cells " << nAspect
273  << endl;
274  }
275 
276  return true;
277  }
278 
279  if (debug || report)
280  {
281  Info<< " Max cell openness = " << maxOpennessCell << " OK." << nl
282  << " Max aspect ratio = " << maxAspectRatio << " OK."
283  << endl;
284  }
285 
286  return false;
287 }
288 
289 
291 (
292  const bool report,
293  labelHashSet* setPtr
294 ) const
295 {
296  if (debug)
297  {
298  Info<< "bool primitiveMesh::checkFaceAreas("
299  << "const bool, labelHashSet*) const: "
300  << "checking face area magnitudes" << endl;
301  }
302 
303  const scalarField magFaceAreas = mag(faceAreas());
304 
305  scalar minArea = GREAT;
306  scalar maxArea = -GREAT;
307 
308  forAll (magFaceAreas, faceI)
309  {
310  if (magFaceAreas[faceI] < VSMALL)
311  {
312  if (setPtr)
313  {
314  setPtr->insert(faceI);
315  }
316  }
317 
318  minArea = min(minArea, magFaceAreas[faceI]);
319  maxArea = max(maxArea, magFaceAreas[faceI]);
320  }
321 
322  reduce(minArea, minOp<scalar>());
323  reduce(maxArea, maxOp<scalar>());
324 
325  if (minArea < VSMALL)
326  {
327  if (debug || report)
328  {
329  Info<< " ***Zero or negative face area detected. "
330  "Minimum area: " << minArea << endl;
331  }
332 
333  return true;
334  }
335  else
336  {
337  if (debug || report)
338  {
339  Info<< " Minumum face area = " << minArea
340  << ". Maximum face area = " << maxArea
341  << ". Face area magnitudes OK." << endl;
342  }
343 
344  return false;
345  }
346 }
347 
348 
350 (
351  const bool report,
352  labelHashSet* setPtr
353 ) const
354 {
355  if (debug)
356  {
357  Info<< "bool primitiveMesh::checkCellVolumes("
358  << "const bool, labelHashSet*) const: "
359  << "checking cell volumes" << endl;
360  }
361 
362  const scalarField& vols = cellVolumes();
363 
364  scalar minVolume = GREAT;
365  scalar maxVolume = -GREAT;
366 
367  label nNegVolCells = 0;
368 
369  forAll (vols, cellI)
370  {
371  if (vols[cellI] < VSMALL)
372  {
373  if (setPtr)
374  {
375  setPtr->insert(cellI);
376  }
377 
378  nNegVolCells++;
379  }
380 
381  minVolume = min(minVolume, vols[cellI]);
382  maxVolume = max(maxVolume, vols[cellI]);
383  }
384 
385  reduce(minVolume, minOp<scalar>());
386  reduce(maxVolume, maxOp<scalar>());
387  reduce(nNegVolCells, sumOp<label>());
388 
389  if (minVolume < VSMALL)
390  {
391  if (debug || report)
392  {
393  Info<< " ***Zero or negative cell volume detected. "
394  << "Minimum negative volume: " << minVolume
395  << ", Number of negative volume cells: " << nNegVolCells
396  << endl;
397  }
398 
399  return true;
400  }
401  else
402  {
403  if (debug || report)
404  {
405  Info<< " Min volume = " << minVolume
406  << ". Max volume = " << maxVolume
407  << ". Total volume = " << gSum(vols)
408  << ". Cell volumes OK." << endl;
409  }
410 
411  return false;
412  }
413 }
414 
415 
417 (
418  const bool report,
419  labelHashSet* setPtr
420 ) const
421 {
422  if (debug)
423  {
424  Info<< "bool primitiveMesh::checkFaceOrthogonality("
425  << "const bool, labelHashSet*) const: "
426  << "checking mesh non-orthogonality" << endl;
427  }
428 
429  // for all internal faces check theat the d dot S product is positive
430  const vectorField& centres = cellCentres();
431  const vectorField& areas = faceAreas();
432 
433  const labelList& own = faceOwner();
434  const labelList& nei = faceNeighbour();
435 
436  // Severe nonorthogonality threshold
437  const scalar severeNonorthogonalityThreshold =
438  ::cos(nonOrthThreshold_/180.0*mathematicalConstant::pi);
439 
440  scalar minDDotS = GREAT;
441 
442  scalar sumDDotS = 0;
443 
444  label severeNonOrth = 0;
445 
446  label errorNonOrth = 0;
447 
448  forAll (nei, faceI)
449  {
450  vector d = centres[nei[faceI]] - centres[own[faceI]];
451  const vector& s = areas[faceI];
452 
453  scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
454 
455  if (dDotS < severeNonorthogonalityThreshold)
456  {
457  if (dDotS > SMALL)
458  {
459  if (setPtr)
460  {
461  setPtr->insert(faceI);
462  }
463 
464  severeNonOrth++;
465  }
466  else
467  {
468  if (setPtr)
469  {
470  setPtr->insert(faceI);
471  }
472 
473  errorNonOrth++;
474  }
475  }
476 
477  if (dDotS < minDDotS)
478  {
479  minDDotS = dDotS;
480  }
481 
482  sumDDotS += dDotS;
483  }
484 
485  reduce(minDDotS, minOp<scalar>());
486  reduce(sumDDotS, sumOp<scalar>());
487  reduce(severeNonOrth, sumOp<label>());
488  reduce(errorNonOrth, sumOp<label>());
489 
490  if (debug || report)
491  {
492  label neiSize = nei.size();
493  reduce(neiSize, sumOp<label>());
494 
495  if (neiSize > 0)
496  {
497  if (debug || report)
498  {
499  Info<< " Mesh non-orthogonality Max: "
500  << ::acos(minDDotS)/mathematicalConstant::pi*180.0
501  << " average: " <<
502  ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0
503  << endl;
504  }
505  }
506 
507  if (severeNonOrth > 0)
508  {
509  Info<< " *Number of severely non-orthogonal faces: "
510  << severeNonOrth << "." << endl;
511  }
512  }
513 
514  if (errorNonOrth > 0)
515  {
516  if (debug || report)
517  {
518  Info<< " ***Number of non-orthogonality errors: "
519  << errorNonOrth << "." << endl;
520  }
521 
522  return true;
523  }
524  else
525  {
526  if (debug || report)
527  {
528  Info<< " Non-orthogonality check OK." << endl;
529  }
530 
531  return false;
532  }
533 }
534 
535 
537 (
538  const bool report,
539  const scalar minPyrVol,
540  labelHashSet* setPtr
541 ) const
542 {
543  if (debug)
544  {
545  Info<< "bool primitiveMesh::checkFacePyramids("
546  << "const bool, const scalar, labelHashSet*) const: "
547  << "checking face orientation" << endl;
548  }
549 
550  // check whether face area vector points to the cell with higher label
551  const vectorField& ctrs = cellCentres();
552 
553  const labelList& own = faceOwner();
554  const labelList& nei = faceNeighbour();
555 
556  const faceList& f = faces();
557 
558  const pointField& p = points();
559 
560  label nErrorPyrs = 0;
561 
562  forAll (f, faceI)
563  {
564  // Create the owner pyramid - it will have negative volume
565  scalar pyrVol = pyramidPointFaceRef(f[faceI], ctrs[own[faceI]]).mag(p);
566 
567  if (pyrVol > -minPyrVol)
568  {
569  if (setPtr)
570  {
571  setPtr->insert(faceI);
572  }
573 
574  nErrorPyrs++;
575  }
576 
577  if (isInternalFace(faceI))
578  {
579  // Create the neighbour pyramid - it will have positive volume
580  scalar pyrVol =
581  pyramidPointFaceRef(f[faceI], ctrs[nei[faceI]]).mag(p);
582 
583  if (pyrVol < minPyrVol)
584  {
585  if (setPtr)
586  {
587  setPtr->insert(faceI);
588  }
589 
590  nErrorPyrs++;
591  }
592  }
593  }
594 
595  reduce(nErrorPyrs, sumOp<label>());
596 
597  if (nErrorPyrs > 0)
598  {
599  if (debug || report)
600  {
601  Info<< " ***Error in face pyramids: "
602  << nErrorPyrs << " faces are incorrectly oriented."
603  << endl;
604  }
605 
606  return true;
607  }
608  else
609  {
610  if (debug || report)
611  {
612  Info<< " Face pyramids OK." << endl;
613  }
614 
615  return false;
616  }
617 }
618 
619 
621 (
622  const bool report,
623  labelHashSet* setPtr
624 ) const
625 {
626  if (debug)
627  {
628  Info<< "bool primitiveMesh::checkFaceSkewnesss("
629  << "const bool, labelHashSet*) const: "
630  << "checking face skewness" << endl;
631  }
632 
633  // Warn if the skew correction vector is more than skewWarning times
634  // larger than the face area vector
635 
636  const pointField& p = points();
637  const faceList& fcs = faces();
638 
639  const labelList& own = faceOwner();
640  const labelList& nei = faceNeighbour();
641  const vectorField& cellCtrs = cellCentres();
642  const vectorField& faceCtrs = faceCentres();
643  const vectorField& fAreas = faceAreas();
644 
645  scalar maxSkew = 0;
646  label nWarnSkew = 0;
647 
648  forAll(nei, faceI)
649  {
650  vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
651  vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
652 
653  // Skewness vector
654  vector sv =
655  Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
656  vector svHat = sv/(mag(sv) + VSMALL);
657 
658  // Normalisation distance calculated as the approximate distance
659  // from the face centre to the edge of the face in the direction of
660  // the skewness
661  scalar fd = 0.2*mag(d) + VSMALL;
662  const face& f = fcs[faceI];
663  forAll(f, pi)
664  {
665  fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
666  }
667 
668  // Normalised skewness
669  scalar skewness = mag(sv)/fd;
670 
671  // Check if the skewness vector is greater than the PN vector.
672  // This does not cause trouble but is a good indication of a poor mesh.
673  if (skewness > skewThreshold_)
674  {
675  if (setPtr)
676  {
677  setPtr->insert(faceI);
678  }
679 
680  nWarnSkew++;
681  }
682 
683  if(skewness > maxSkew)
684  {
685  maxSkew = skewness;
686  }
687  }
688 
689 
690  // Boundary faces: consider them to have only skewness error.
691  // (i.e. treat as if mirror cell on other side)
692 
693  for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
694  {
695  vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
696 
697  vector normal = fAreas[faceI];
698  normal /= mag(normal) + VSMALL;
699  vector d = normal*(normal & Cpf);
700 
701 
702  // Skewness vector
703  vector sv = Cpf - ((fAreas[faceI]&Cpf)/((fAreas[faceI]&d)+VSMALL))*d;
704  vector svHat = sv/(mag(sv) + VSMALL);
705 
706  // Normalisation distance calculated as the approximate distance
707  // from the face centre to the edge of the face in the direction of
708  // the skewness
709  scalar fd = 0.4*mag(d) + VSMALL;
710  const face& f = fcs[faceI];
711  forAll(f, pi)
712  {
713  fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
714  }
715 
716  // Normalised skewness
717  scalar skewness = mag(sv)/fd;
718 
719  // Check if the skewness vector is greater than the PN vector.
720  // This does not cause trouble but is a good indication of a poor mesh.
721  if (skewness > skewThreshold_)
722  {
723  if (setPtr)
724  {
725  setPtr->insert(faceI);
726  }
727 
728  nWarnSkew++;
729  }
730 
731  if(skewness > maxSkew)
732  {
733  maxSkew = skewness;
734  }
735  }
736 
737 
738  reduce(maxSkew, maxOp<scalar>());
739  reduce(nWarnSkew, sumOp<label>());
740 
741  if (nWarnSkew > 0)
742  {
743  if (debug || report)
744  {
745  Info<< " ***Max skewness = " << maxSkew
746  << ", " << nWarnSkew << " highly skew faces detected"
747  " which may impair the quality of the results"
748  << endl;
749  }
750 
751  return true;
752  }
753  else
754  {
755  if (debug || report)
756  {
757  Info<< " Max skewness = " << maxSkew << " OK." << endl;
758  }
759 
760  return false;
761  }
762 }
763 
764 
766 (
767  const bool report,
768  labelHashSet* setPtr
769 ) const
770 {
771  if (debug)
772  {
773  Info<< "bool primitiveMesh::checkPoints"
774  << "(const bool, labelHashSet*) const: "
775  << "checking points" << endl;
776  }
777 
778  label nFaceErrors = 0;
779  label nCellErrors = 0;
780 
781  const labelListList& pf = pointFaces();
782 
783  forAll (pf, pointI)
784  {
785  if (pf[pointI].empty())
786  {
787  if (setPtr)
788  {
789  setPtr->insert(pointI);
790  }
791 
792  nFaceErrors++;
793  }
794  }
795 
796 
797  forAll (pf, pointI)
798  {
799  const labelList& pc = pointCells(pointI);
800 
801  if (pc.empty())
802  {
803  if (setPtr)
804  {
805  setPtr->insert(pointI);
806  }
807 
808  nCellErrors++;
809  }
810  }
811 
812  reduce(nFaceErrors, sumOp<label>());
813  reduce(nCellErrors, sumOp<label>());
814 
815  if (nFaceErrors > 0 || nCellErrors > 0)
816  {
817  if (debug || report)
818  {
819  Info<< " ***Unused points found in the mesh, "
820  "number unused by faces: " << nFaceErrors
821  << " number unused by cells: " << nCellErrors
822  << endl;
823  }
824 
825  return true;
826  }
827  else
828  {
829  if (debug || report)
830  {
831  Info<< " Point usage OK." << endl;
832  }
833 
834  return false;
835  }
836 }
837 
838 
839 // Check convexity of angles in a face. Allow a slight non-convexity.
840 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
841 // (if truly concave and points not visible from face centre the face-pyramid
842 // check in checkMesh will fail)
844 (
845  const bool report,
846  const scalar maxDeg,
847  labelHashSet* setPtr
848 ) const
849 {
850  if (debug)
851  {
852  Info<< "bool primitiveMesh::checkFaceAngles"
853  << "(const bool, const scalar, labelHashSet*) const: "
854  << "checking face angles" << endl;
855  }
856 
857  if (maxDeg < -SMALL || maxDeg > 180+SMALL)
858  {
860  (
861  "primitiveMesh::checkFaceAngles"
862  "(const bool, const scalar, labelHashSet*)"
863  ) << "maxDeg should be [0..180] but is now " << maxDeg
864  << exit(FatalError);
865  }
866 
867  const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
868 
869  const pointField& p = points();
870  const faceList& fcs = faces();
871  vectorField faceNormals(faceAreas());
872  faceNormals /= mag(faceNormals) + VSMALL;
873 
874  scalar maxEdgeSin = 0.0;
875 
876  label nConcave = 0;
877 
878  label errorFaceI = -1;
879 
880  forAll(fcs, faceI)
881  {
882  const face& f = fcs[faceI];
883 
884  // Get edge from f[0] to f[size-1];
885  vector ePrev(p[f[0]] - p[f[f.size()-1]]);
886  scalar magEPrev = mag(ePrev);
887  ePrev /= magEPrev + VSMALL;
888 
889  forAll(f, fp0)
890  {
891  // Get vertex after fp
892  label fp1 = f.fcIndex(fp0);
893 
894  // Normalized vector between two consecutive points
895  vector e10(p[f[fp1]] - p[f[fp0]]);
896  scalar magE10 = mag(e10);
897  e10 /= magE10 + VSMALL;
898 
899  if (magEPrev > SMALL && magE10 > SMALL)
900  {
901  vector edgeNormal = ePrev ^ e10;
902  scalar magEdgeNormal = mag(edgeNormal);
903 
904  if (magEdgeNormal < maxSin)
905  {
906  // Edges (almost) aligned -> face is ok.
907  }
908  else
909  {
910  // Check normal
911  edgeNormal /= magEdgeNormal;
912 
913  if ((edgeNormal & faceNormals[faceI]) < SMALL)
914  {
915  if (faceI != errorFaceI)
916  {
917  // Count only one error per face.
918  errorFaceI = faceI;
919  nConcave++;
920  }
921 
922  if (setPtr)
923  {
924  setPtr->insert(faceI);
925  }
926 
927  maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
928  }
929  }
930  }
931 
932  ePrev = e10;
933  magEPrev = magE10;
934  }
935  }
936 
937  reduce(nConcave, sumOp<label>());
938  reduce(maxEdgeSin, maxOp<scalar>());
939 
940  if (nConcave > 0)
941  {
942  scalar maxConcaveDegr =
943  Foam::asin(Foam::min(1.0, maxEdgeSin))
945 
946  if (debug || report)
947  {
948  Info<< " *There are " << nConcave
949  << " faces with concave angles between consecutive"
950  << " edges. Max concave angle = " << maxConcaveDegr
951  << " degrees." << endl;
952  }
953 
954  return true;
955  }
956  else
957  {
958  if (debug || report)
959  {
960  Info<< " All angles in faces OK." << endl;
961  }
962 
963  return false;
964  }
965 }
966 
967 
968 // Check warpage of faces. Is calculated as the difference between areas of
969 // individual triangles and the overall area of the face (which ifself is
970 // is the average of the areas of the individual triangles).
972 (
973  const bool report,
974  const scalar warnFlatness,
975  labelHashSet* setPtr
976 ) const
977 {
978  if (debug)
979  {
980  Info<< "bool primitiveMesh::checkFaceFlatness"
981  << "(const bool, const scalar, labelHashSet*) const: "
982  << "checking face flatness" << endl;
983  }
984 
985  if (warnFlatness < 0 || warnFlatness > 1)
986  {
988  (
989  "primitiveMesh::checkFaceFlatness"
990  "(const bool, const scalar, labelHashSet*)"
991  ) << "warnFlatness should be [0..1] but is now " << warnFlatness
992  << exit(FatalError);
993  }
994 
995 
996  const pointField& p = points();
997  const faceList& fcs = faces();
998  const pointField& fctrs = faceCentres();
999 
1000  // Areas are calculated as the sum of areas. (see
1001  // primitiveMeshFaceCentresAndAreas.C)
1002  scalarField magAreas(mag(faceAreas()));
1003 
1004  label nWarped = 0;
1005 
1006  scalar minFlatness = GREAT;
1007  scalar sumFlatness = 0;
1008  label nSummed = 0;
1009 
1010  forAll(fcs, faceI)
1011  {
1012  const face& f = fcs[faceI];
1013 
1014  if (f.size() > 3 && magAreas[faceI] > VSMALL)
1015  {
1016  const point& fc = fctrs[faceI];
1017 
1018  // Calculate the sum of magnitude of areas and compare to magnitude
1019  // of sum of areas.
1020 
1021  scalar sumA = 0.0;
1022 
1023  forAll(f, fp)
1024  {
1025  const point& thisPoint = p[f[fp]];
1026  const point& nextPoint = p[f.nextLabel(fp)];
1027 
1028  // Triangle around fc.
1029  vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
1030  sumA += mag(n);
1031  }
1032 
1033  scalar flatness = magAreas[faceI] / (sumA+VSMALL);
1034 
1035  sumFlatness += flatness;
1036  nSummed++;
1037 
1038  minFlatness = min(minFlatness, flatness);
1039 
1040  if (flatness < warnFlatness)
1041  {
1042  nWarped++;
1043 
1044  if (setPtr)
1045  {
1046  setPtr->insert(faceI);
1047  }
1048  }
1049  }
1050  }
1051 
1052 
1053  reduce(nWarped, sumOp<label>());
1054  reduce(minFlatness, minOp<scalar>());
1055 
1056  reduce(nSummed, sumOp<label>());
1057  reduce(sumFlatness, sumOp<scalar>());
1058 
1059  if (debug || report)
1060  {
1061  if (nSummed > 0)
1062  {
1063  Info<< " Face flatness (1 = flat, 0 = butterfly) : average = "
1064  << sumFlatness / nSummed << " min = " << minFlatness << endl;
1065  }
1066  }
1067 
1068 
1069  if (nWarped> 0)
1070  {
1071  if (debug || report)
1072  {
1073  Info<< " *There are " << nWarped
1074  << " faces with ratio between projected and actual area < "
1075  << warnFlatness << endl;
1076 
1077  Info<< " Minimum ratio (minimum flatness, maximum warpage) = "
1078  << minFlatness << endl;
1079  }
1080 
1081  return true;
1082  }
1083  else
1084  {
1085  if (debug || report)
1086  {
1087  Info<< " All face flatness OK." << endl;
1088  }
1089 
1090  return false;
1091  }
1092 }
1093 
1094 
1095 // Check 1D/2Dness of edges. Gets passed the non-empty directions and
1096 // checks all edges in the mesh whether they:
1097 // - have no component in a non-empty direction or
1098 // - are only in a singe non-empty direction.
1099 // Empty direction info is passed in as a vector of labels (synchronised)
1100 // which are 1 if the direction is non-empty, 0 if it is.
1103  const bool report,
1104  const Vector<label>& directions,
1105  labelHashSet* setPtr
1106 ) const
1107 {
1108  if (debug)
1109  {
1110  Info<< "bool primitiveMesh::checkEdgeAlignment("
1111  << "const bool, const Vector<label>&, labelHashSet*) const: "
1112  << "checking edge alignment" << endl;
1113  }
1114 
1115  label nDirs = 0;
1116  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
1117  {
1118  if (directions[cmpt] == 1)
1119  {
1120  nDirs++;
1121  }
1122  else if (directions[cmpt] != 0)
1123  {
1124  FatalErrorIn
1125  (
1126  "primitiveMesh::checkEdgeAlignment"
1127  "(const bool, const Vector<label>&, labelHashSet*)"
1128  ) << "directions should contain 0 or 1 but is now " << directions
1129  << exit(FatalError);
1130  }
1131  }
1132 
1133  if (nDirs == vector::nComponents)
1134  {
1135  return false;
1136  }
1137 
1138 
1139 
1140  const pointField& p = points();
1141  const faceList& fcs = faces();
1142 
1143  EdgeMap<label> edgesInError;
1144 
1145  forAll(fcs, faceI)
1146  {
1147  const face& f = fcs[faceI];
1148 
1149  forAll(f, fp)
1150  {
1151  label p0 = f[fp];
1152  label p1 = f.nextLabel(fp);
1153  if (p0 < p1)
1154  {
1155  vector d(p[p1]-p[p0]);
1156  scalar magD = mag(d);
1157 
1158  if (magD > ROOTVSMALL)
1159  {
1160  d /= magD;
1161 
1162  // Check how many empty directions are used by the edge.
1163  label nEmptyDirs = 0;
1164  label nNonEmptyDirs = 0;
1165  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
1166  {
1167  if (mag(d[cmpt]) > 1e-6)
1168  {
1169  if (directions[cmpt] == 0)
1170  {
1171  nEmptyDirs++;
1172  }
1173  else
1174  {
1175  nNonEmptyDirs++;
1176  }
1177  }
1178  }
1179 
1180  if (nEmptyDirs == 0)
1181  {
1182  // Purely in ok directions.
1183  }
1184  else if (nEmptyDirs == 1)
1185  {
1186  // Ok if purely in empty directions.
1187  if (nNonEmptyDirs > 0)
1188  {
1189  edgesInError.insert(edge(p0, p1), faceI);
1190  }
1191  }
1192  else if (nEmptyDirs > 1)
1193  {
1194  // Always an error
1195  edgesInError.insert(edge(p0, p1), faceI);
1196  }
1197  }
1198  }
1199  }
1200  }
1201 
1202  label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
1203 
1204  if (nErrorEdges > 0)
1205  {
1206  if (debug || report)
1207  {
1208  Info<< " ***Number of edges not aligned with or perpendicular to "
1209  << "non-empty directions: " << nErrorEdges << endl;
1210  }
1211 
1212  if (setPtr)
1213  {
1214  setPtr->resize(2*edgesInError.size());
1215  forAllConstIter(EdgeMap<label>, edgesInError, iter)
1216  {
1217  setPtr->insert(iter.key()[0]);
1218  setPtr->insert(iter.key()[1]);
1219  }
1220  }
1221 
1222  return true;
1223  }
1224  else
1225  {
1226  if (debug || report)
1227  {
1228  Info<< " All edges aligned with or perpendicular to "
1229  << "non-empty directions." << endl;
1230  }
1231  return false;
1232  }
1233 }
1234 
1235 
1238  const bool report,
1239  labelHashSet* setPtr
1240 ) const
1241 {
1242  if (debug)
1243  {
1244  Info<< "bool primitiveMesh::checkUpperTriangular("
1245  << "const bool, labelHashSet*) const: "
1246  << "checking face ordering" << endl;
1247  }
1248 
1249  // Check whether internal faces are ordered in the upper triangular order
1250  const labelList& own = faceOwner();
1251  const labelList& nei = faceNeighbour();
1252 
1253  const cellList& c = cells();
1254 
1255  label internal = nInternalFaces();
1256 
1257  // Has error occurred?
1258  bool error = false;
1259  // Have multiple faces been detected?
1260  label nMultipleCells = false;
1261 
1262  // Loop through faceCells once more and make sure that for internal cell
1263  // the first label is smaller
1264  for (label faceI = 0; faceI < internal; faceI++)
1265  {
1266  if (own[faceI] >= nei[faceI])
1267  {
1268  error = true;
1269 
1270  if (setPtr)
1271  {
1272  setPtr->insert(faceI);
1273  }
1274  }
1275  }
1276 
1277  // Loop through all cells. For each cell, find the face that is internal
1278  // and add it to the check list (upper triangular order).
1279  // Once the list is completed, check it against the faceCell list
1280 
1281  forAll (c, cellI)
1282  {
1283  const labelList& curFaces = c[cellI];
1284 
1285  // Neighbouring cells
1286  SortableList<label> nbr(curFaces.size());
1287 
1288  forAll(curFaces, i)
1289  {
1290  label faceI = curFaces[i];
1291 
1292  if (faceI >= nInternalFaces())
1293  {
1294  // Sort last
1295  nbr[i] = labelMax;
1296  }
1297  else
1298  {
1299  label nbrCellI = nei[faceI];
1300 
1301  if (nbrCellI == cellI)
1302  {
1303  nbrCellI = own[faceI];
1304  }
1305 
1306  if (cellI < nbrCellI)
1307  {
1308  // cellI is master
1309  nbr[i] = nbrCellI;
1310  }
1311  else
1312  {
1313  // nbrCell is master. Let it handle this face.
1314  nbr[i] = labelMax;
1315  }
1316  }
1317  }
1318 
1319  nbr.sort();
1320 
1321  // Now nbr holds the cellCells in incremental order. Check:
1322  // - neighbouring cells appear only once. Since nbr is sorted this
1323  // is simple check on consecutive elements
1324  // - faces indexed in same order as nbr are incrementing as well.
1325 
1326  label prevCell = nbr[0];
1327  label prevFace = curFaces[nbr.indices()[0]];
1328 
1329  bool hasMultipleFaces = false;
1330 
1331  for (label i = 1; i < nbr.size(); i++)
1332  {
1333  label thisCell = nbr[i];
1334  label thisFace = curFaces[nbr.indices()[i]];
1335 
1336  if (thisCell == labelMax)
1337  {
1338  break;
1339  }
1340 
1341  if (thisCell == prevCell)
1342  {
1343  hasMultipleFaces = true;
1344 
1345  if (setPtr)
1346  {
1347  setPtr->insert(prevFace);
1348  setPtr->insert(thisFace);
1349  }
1350  }
1351  else if (thisFace < prevFace)
1352  {
1353  error = true;
1354 
1355  if (setPtr)
1356  {
1357  setPtr->insert(thisFace);
1358  }
1359  }
1360 
1361  prevCell = thisCell;
1362  prevFace = thisFace;
1363  }
1364 
1365  if (hasMultipleFaces)
1366  {
1367  nMultipleCells++;
1368  }
1369  }
1370 
1371  reduce(error, orOp<bool>());
1372  reduce(nMultipleCells, sumOp<label>());
1373 
1374  if ((debug || report) && nMultipleCells > 0)
1375  {
1376  Info<< " <<Found " << nMultipleCells
1377  << " neighbouring cells with multiple inbetween faces." << endl;
1378  }
1379 
1380  if (error)
1381  {
1382  if (debug || report)
1383  {
1384  Info<< " ***Faces not in upper triangular order." << endl;
1385  }
1386 
1387  return true;
1388  }
1389  else
1390  {
1391  if (debug || report)
1392  {
1393  Info<< " Upper triangular ordering OK." << endl;
1394  }
1395 
1396  return false;
1397  }
1398 }
1399 
1400 
1403  const bool report,
1404  labelHashSet* setPtr
1405 ) const
1406 {
1407  if (debug)
1408  {
1409  Info<< "bool primitiveMesh::checkCellsZipUp("
1410  << "const bool, labelHashSet*) const: "
1411  << "checking topological cell openness" << endl;
1412  }
1413 
1414  label nOpenCells = 0;
1415 
1416  const faceList& f = faces();
1417  const cellList& c = cells();
1418 
1419  forAll (c, cellI)
1420  {
1421  const labelList& curFaces = c[cellI];
1422 
1423  const edgeList cellEdges = c[cellI].edges(f);
1424 
1425  labelList edgeUsage(cellEdges.size(), 0);
1426 
1427  forAll (curFaces, faceI)
1428  {
1429  edgeList curFaceEdges = f[curFaces[faceI]].edges();
1430 
1431  forAll (curFaceEdges, faceEdgeI)
1432  {
1433  const edge& curEdge = curFaceEdges[faceEdgeI];
1434 
1435  forAll (cellEdges, cellEdgeI)
1436  {
1437  if (cellEdges[cellEdgeI] == curEdge)
1438  {
1439  edgeUsage[cellEdgeI]++;
1440  break;
1441  }
1442  }
1443  }
1444  }
1445 
1446  edgeList singleEdges(cellEdges.size());
1447  label nSingleEdges = 0;
1448 
1449  forAll (edgeUsage, edgeI)
1450  {
1451  if (edgeUsage[edgeI] == 1)
1452  {
1453  singleEdges[nSingleEdges] = cellEdges[edgeI];
1454  nSingleEdges++;
1455  }
1456  else if (edgeUsage[edgeI] != 2)
1457  {
1458  if (setPtr)
1459  {
1460  setPtr->insert(cellI);
1461  }
1462  }
1463  }
1464 
1465  if (nSingleEdges > 0)
1466  {
1467  if (setPtr)
1468  {
1469  setPtr->insert(cellI);
1470  }
1471 
1472  nOpenCells++;
1473  }
1474  }
1475 
1476  reduce(nOpenCells, sumOp<label>());
1477 
1478  if (nOpenCells > 0)
1479  {
1480  if (debug || report)
1481  {
1482  Info<< " ***Open cells found, number of cells: " << nOpenCells
1483  << ". This problem may be fixable using the zipUpMesh utility."
1484  << endl;
1485  }
1486 
1487  return true;
1488  }
1489  else
1490  {
1491  if (debug || report)
1492  {
1493  Info<< " Topological cell zip-up check OK." << endl;
1494  }
1495 
1496  return false;
1497  }
1498 }
1499 
1500 
1501 // Vertices of face within point range and unique.
1504  const bool report,
1505  labelHashSet* setPtr
1506 ) const
1507 {
1508  if (debug)
1509  {
1510  Info<< "bool primitiveMesh::checkFaceVertices("
1511  << "const bool, labelHashSet*) const: "
1512  << "checking face vertices" << endl;
1513  }
1514 
1515  // Check that all vertex labels are valid
1516  const faceList& f = faces();
1517 
1518  label nErrorFaces = 0;
1519 
1520  forAll (f, fI)
1521  {
1522  const face& curFace = f[fI];
1523 
1524  if (min(curFace) < 0 || max(curFace) > nPoints())
1525  {
1526  if (setPtr)
1527  {
1528  setPtr->insert(fI);
1529  }
1530 
1531  nErrorFaces++;
1532  }
1533 
1534  // Uniqueness of vertices
1535  labelHashSet facePoints(2*curFace.size());
1536 
1537  forAll(curFace, fp)
1538  {
1539  bool inserted = facePoints.insert(curFace[fp]);
1540 
1541  if (!inserted)
1542  {
1543  if (setPtr)
1544  {
1545  setPtr->insert(fI);
1546  }
1547 
1548  nErrorFaces++;
1549  }
1550  }
1551  }
1552 
1553  reduce(nErrorFaces, sumOp<label>());
1554 
1555  if (nErrorFaces > 0)
1556  {
1557  if (debug || report)
1558  {
1559  Info<< " Faces with invalid vertex labels found, "
1560  << " number of faces: " << nErrorFaces << endl;
1561  }
1562 
1563  return true;
1564  }
1565  else
1566  {
1567  if (debug || report)
1568  {
1569  Info<< " Face vertices OK." << endl;
1570  }
1571 
1572  return false;
1573  }
1574 }
1575 
1576 
1577 // Check if all points on face are shared between faces.
1578 bool Foam::primitiveMesh::checkDuplicateFaces
1579 (
1580  const label faceI,
1581  const Map<label>& nCommonPoints,
1582  label& nBaffleFaces,
1583  labelHashSet* setPtr
1584 ) const
1585 {
1586  bool error = false;
1587 
1588  forAllConstIter(Map<label>, nCommonPoints, iter)
1589  {
1590  label nbFaceI = iter.key();
1591  label nCommon = iter();
1592 
1593  const face& curFace = faces()[faceI];
1594  const face& nbFace = faces()[nbFaceI];
1595 
1596  if (nCommon == nbFace.size() || nCommon == curFace.size())
1597  {
1598  if (nbFace.size() != curFace.size())
1599  {
1600  error = true;
1601  }
1602  else
1603  {
1604  nBaffleFaces++;
1605  }
1606 
1607  if (setPtr)
1608  {
1609  setPtr->insert(faceI);
1610  setPtr->insert(nbFaceI);
1611  }
1612  }
1613  }
1614 
1615  return error;
1616 }
1617 
1618 
1619 // Check that shared points are in consecutive order.
1620 bool Foam::primitiveMesh::checkCommonOrder
1621 (
1622  const label faceI,
1623  const Map<label>& nCommonPoints,
1624  labelHashSet* setPtr
1625 ) const
1626 {
1627  bool error = false;
1628 
1629  forAllConstIter(Map<label>, nCommonPoints, iter)
1630  {
1631  label nbFaceI = iter.key();
1632  label nCommon = iter();
1633 
1634  const face& curFace = faces()[faceI];
1635  const face& nbFace = faces()[nbFaceI];
1636 
1637  if
1638  (
1639  nCommon >= 2
1640  && nCommon != nbFace.size()
1641  && nCommon != curFace.size()
1642  )
1643  {
1644  forAll(curFace, fp)
1645  {
1646  // Get the index in the neighbouring face shared with curFace
1647  label nb = findIndex(nbFace, curFace[fp]);
1648 
1649  if (nb != -1)
1650  {
1651 
1652  // Check the whole face from nb onwards for shared vertices
1653  // with neighbouring face. Rule is that any shared vertices
1654  // should be consecutive on both faces i.e. if they are
1655  // vertices fp,fp+1,fp+2 on one face they should be
1656  // vertices nb, nb+1, nb+2 (or nb+2, nb+1, nb) on the
1657  // other face.
1658 
1659 
1660  // Vertices before and after on curFace
1661  label fpPlus1 = curFace.fcIndex(fp);
1662  label fpMin1 = curFace.rcIndex(fp);
1663 
1664  // Vertices before and after on nbFace
1665  label nbPlus1 = nbFace.fcIndex(nb);
1666  label nbMin1 = nbFace.rcIndex(nb);
1667 
1668  // Find order of walking by comparing next points on both
1669  // faces.
1670  label curInc = labelMax;
1671  label nbInc = labelMax;
1672 
1673  if (nbFace[nbPlus1] == curFace[fpPlus1])
1674  {
1675  curInc = 1;
1676  nbInc = 1;
1677  }
1678  else if (nbFace[nbPlus1] == curFace[fpMin1])
1679  {
1680  curInc = -1;
1681  nbInc = 1;
1682  }
1683  else if (nbFace[nbMin1] == curFace[fpMin1])
1684  {
1685  curInc = -1;
1686  nbInc = -1;
1687  }
1688  else
1689  {
1690  curInc = 1;
1691  nbInc = -1;
1692  }
1693 
1694 
1695  // Pass1: loop until start of common vertices found.
1696  label curNb = nb;
1697  label curFp = fp;
1698 
1699  do
1700  {
1701  curFp += curInc;
1702 
1703  if (curFp >= curFace.size())
1704  {
1705  curFp = 0;
1706  }
1707  else if (curFp < 0)
1708  {
1709  curFp = curFace.size()-1;
1710  }
1711 
1712  curNb += nbInc;
1713 
1714  if (curNb >= nbFace.size())
1715  {
1716  curNb = 0;
1717  }
1718  else if (curNb < 0)
1719  {
1720  curNb = nbFace.size()-1;
1721  }
1722  } while (curFace[curFp] == nbFace[curNb]);
1723 
1724 
1725  // Pass2: check equality walking from curFp, curNb
1726  // in opposite order.
1727 
1728  curInc = -curInc;
1729  nbInc = -nbInc;
1730 
1731  for (label commonI = 0; commonI < nCommon; commonI++)
1732  {
1733  curFp += curInc;
1734 
1735  if (curFp >= curFace.size())
1736  {
1737  curFp = 0;
1738  }
1739  else if (curFp < 0)
1740  {
1741  curFp = curFace.size()-1;
1742  }
1743 
1744  curNb += nbInc;
1745 
1746  if (curNb >= nbFace.size())
1747  {
1748  curNb = 0;
1749  }
1750  else if (curNb < 0)
1751  {
1752  curNb = nbFace.size()-1;
1753  }
1754 
1755  if (curFace[curFp] != nbFace[curNb])
1756  {
1757  if (setPtr)
1758  {
1759  setPtr->insert(faceI);
1760  setPtr->insert(nbFaceI);
1761  }
1762 
1763  error = true;
1764 
1765  break;
1766  }
1767  }
1768 
1769 
1770  // Done the curFace - nbFace combination.
1771  break;
1772  }
1773  }
1774  }
1775  }
1776 
1777  return error;
1778 }
1779 
1780 
1781 // Checks common vertices between faces. If more than 2 they should be
1782 // consecutive on both faces.
1785  const bool report,
1786  labelHashSet* setPtr
1787 ) const
1788 {
1789  if (debug)
1790  {
1791  Info<< "bool primitiveMesh::checkFaceFaces(const bool, labelHashSet*)"
1792  << " const: " << "checking face-face connectivity" << endl;
1793  }
1794 
1795  const labelListList& pf = pointFaces();
1796 
1797  label nBaffleFaces = 0;
1798  label nErrorDuplicate = 0;
1799  label nErrorOrder = 0;
1800  Map<label> nCommonPoints(100);
1801 
1802  for (label faceI = 0; faceI < nFaces(); faceI++)
1803  {
1804  const face& curFace = faces()[faceI];
1805 
1806  // Calculate number of common points between current faceI and
1807  // neighbouring face. Store on map.
1808  nCommonPoints.clear();
1809 
1810  forAll(curFace, fp)
1811  {
1812  label pointI = curFace[fp];
1813 
1814  const labelList& nbs = pf[pointI];
1815 
1816  forAll(nbs, nbI)
1817  {
1818  label nbFaceI = nbs[nbI];
1819 
1820  if (faceI < nbFaceI)
1821  {
1822  // Only check once for each combination of two faces.
1823 
1824  Map<label>::iterator fnd = nCommonPoints.find(nbFaceI);
1825 
1826  if (fnd == nCommonPoints.end())
1827  {
1828  // First common vertex found.
1829  nCommonPoints.insert(nbFaceI, 1);
1830  }
1831  else
1832  {
1833  fnd()++;
1834  }
1835  }
1836  }
1837  }
1838 
1839  // Perform various checks on common points
1840 
1841  // Check all vertices shared (duplicate point)
1842  if (checkDuplicateFaces(faceI, nCommonPoints, nBaffleFaces, setPtr))
1843  {
1844  nErrorDuplicate++;
1845  }
1846 
1847  // Check common vertices are consecutive on both faces
1848  if (checkCommonOrder(faceI, nCommonPoints, setPtr))
1849  {
1850  nErrorOrder++;
1851  }
1852  }
1853 
1854  reduce(nBaffleFaces, sumOp<label>());
1855  reduce(nErrorDuplicate, sumOp<label>());
1856  reduce(nErrorOrder, sumOp<label>());
1857 
1858  if (nBaffleFaces)
1859  {
1860  Info<< " Number of identical duplicate faces (baffle faces): "
1861  << nBaffleFaces << endl;
1862  }
1863 
1864  if (nErrorDuplicate > 0 || nErrorOrder > 0)
1865  {
1866  if (nErrorDuplicate > 0)
1867  {
1868  Info<< " ***Number of duplicate (not baffle) faces found: "
1869  << nErrorDuplicate << endl;
1870  }
1871 
1872  if (nErrorOrder > 0)
1873  {
1874  Info<< " ***Number of faces with non-consecutive shared points: "
1875  << nErrorOrder << endl;
1876  }
1877 
1878  return true;
1879  }
1880  else
1881  {
1882  if (debug || report)
1883  {
1884  Info<< " Face-face connectivity OK." << endl;
1885  }
1886 
1887  return false;
1888  }
1889 }
1890 
1891 
1892 // Checks cells with 1 or less internal faces. Give numerical problems.
1895  const bool report, // report,
1896  labelHashSet* setPtr, // setPtr
1897  const Vector<label>& meshD
1898 ) const
1899 {
1900  if (debug)
1901  {
1902  Info<< "bool primitiveMesh::checkCellDeterminant(const bool"
1903  << ", labelHashSet*) const: "
1904  << "checking for under-determined cells" << endl;
1905  }
1906 
1907  // Determine number of dimensions and (for 2D) missing dimension
1908  label nDims = 0;
1909  label twoD = -1;
1910  for (direction dir = 0; dir < vector::nComponents; dir++)
1911  {
1912  if (meshD[dir] == 1)
1913  {
1914  nDims++;
1915  }
1916  else
1917  {
1918  twoD = dir;
1919  }
1920  }
1921 
1922 
1923  const cellList& c = cells();
1924 
1925  label nErrorCells = 0;
1926 
1927  scalar minDet = GREAT;
1928  scalar sumDet = 0;
1929  label nSummed = 0;
1930 
1931  if (nDims == 1)
1932  {
1933  minDet = 1;
1934  sumDet = c.size()*minDet;
1935  nSummed = c.size();
1936  }
1937  else
1938  {
1939  forAll (c, cellI)
1940  {
1941  const labelList& curFaces = c[cellI];
1942 
1943  // Calculate local normalization factor
1944  scalar avgArea = 0;
1945 
1946  label nInternalFaces = 0;
1947 
1948  forAll(curFaces, i)
1949  {
1950  if (isInternalFace(curFaces[i]))
1951  {
1952  avgArea += mag(faceAreas()[curFaces[i]]);
1953 
1954  nInternalFaces++;
1955  }
1956  }
1957 
1958  if (nInternalFaces == 0)
1959  {
1960  if (setPtr)
1961  {
1962  setPtr->insert(cellI);
1963  }
1964 
1965  nErrorCells++;
1966  }
1967  else
1968  {
1969  avgArea /= nInternalFaces;
1970 
1971  symmTensor areaTensor(symmTensor::zero);
1972 
1973  forAll(curFaces, i)
1974  {
1975  if (isInternalFace(curFaces[i]))
1976  {
1977  areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
1978  }
1979  }
1980 
1981  if (nDims == 2)
1982  {
1983  // Add the missing eigenvector (such that it does not
1984  // affect the determinant)
1985  if (twoD == 0)
1986  {
1987  areaTensor.xx() = 1;
1988  }
1989  else if (twoD == 1)
1990  {
1991  areaTensor.yy() = 1;
1992  }
1993  else
1994  {
1995  areaTensor.zz() = 1;
1996  }
1997  }
1998 
1999  scalar determinant = mag(det(areaTensor));
2000 
2001  minDet = min(determinant, minDet);
2002  sumDet += determinant;
2003  nSummed++;
2004 
2005  if (determinant < 1e-3)
2006  {
2007  if (setPtr)
2008  {
2009  setPtr->insert(cellI);
2010  }
2011 
2012  nErrorCells++;
2013  }
2014  }
2015  }
2016  }
2017 
2018  reduce(nErrorCells, sumOp<label>());
2019  reduce(minDet, minOp<scalar>());
2020  reduce(sumDet, sumOp<scalar>());
2021  reduce(nSummed, sumOp<label>());
2022 
2023  if (debug || report)
2024  {
2025  if (nSummed > 0)
2026  {
2027  Info<< " Cell determinant (wellposedness) : minimum: " << minDet
2028  << " average: " << sumDet/nSummed
2029  << endl;
2030  }
2031  }
2032 
2033  if (nErrorCells > 0)
2034  {
2035  if (debug || report)
2036  {
2037  Info<< " ***Cells with small determinant found, number of cells: "
2038  << nErrorCells << endl;
2039  }
2040 
2041  return true;
2042  }
2043  else
2044  {
2045  if (debug || report)
2046  {
2047  Info<< " Cell determinant check OK." << endl;
2048  }
2049 
2050  return false;
2051  }
2052 
2053  return false;
2054 }
2055 
2056 
2057 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2058 
2059 bool Foam::primitiveMesh::checkTopology(const bool report) const
2060 {
2061  label noFailedChecks = 0;
2062 
2063  if (checkPoints(report)) noFailedChecks++;
2064  if (checkUpperTriangular(report)) noFailedChecks++;
2065  if (checkCellsZipUp(report)) noFailedChecks++;
2066  if (checkFaceVertices(report)) noFailedChecks++;
2067  if (checkFaceFaces(report)) noFailedChecks++;
2068  //if (checkCellDeterminant(report)) noFailedChecks++;
2069 
2070  if (noFailedChecks == 0)
2071  {
2072  if (debug || report)
2073  {
2074  Info<< " Mesh topology OK." << endl;
2075  }
2076 
2077  return false;
2078  }
2079  else
2080  {
2081  if (debug || report)
2082  {
2083  Info<< " Failed " << noFailedChecks
2084  << " mesh topology checks." << endl;
2085  }
2086 
2087  return true;
2088  }
2089 }
2090 
2091 
2092 bool Foam::primitiveMesh::checkGeometry(const bool report) const
2093 {
2094  label noFailedChecks = 0;
2095 
2096  if (checkClosedBoundary(report)) noFailedChecks++;
2097  if (checkClosedCells(report)) noFailedChecks++;
2098  if (checkFaceAreas(report)) noFailedChecks++;
2099  if (checkCellVolumes(report)) noFailedChecks++;
2100  if (checkFaceOrthogonality(report)) noFailedChecks++;
2101  if (checkFacePyramids(report)) noFailedChecks++;
2102  if (checkFaceSkewness(report)) noFailedChecks++;
2103 
2104  if (noFailedChecks == 0)
2105  {
2106  if (debug || report)
2107  {
2108  Info<< " Mesh geometry OK." << endl;
2109  }
2110  return false;
2111  }
2112  else
2113  {
2114  if (debug || report)
2115  {
2116  Info<< " Failed " << noFailedChecks
2117  << " mesh geometry checks." << endl;
2118  }
2119 
2120  return true;
2121  }
2122 }
2123 
2124 
2125 bool Foam::primitiveMesh::checkMesh(const bool report) const
2126 {
2127  if (debug)
2128  {
2129  Info<< "bool primitiveMesh::checkMesh(const bool report) const: "
2130  << "checking primitiveMesh" << endl;
2131  }
2132 
2133  label noFailedChecks = checkTopology(report) + checkGeometry(report);
2134 
2135  if (noFailedChecks == 0)
2136  {
2137  if (debug || report)
2138  {
2139  Info<< "Mesh OK." << endl;
2140  }
2141 
2142  return false;
2143  }
2144  else
2145  {
2146  if (debug || report)
2147  {
2148  Info<< " Failed " << noFailedChecks
2149  << " mesh checks." << endl;
2150  }
2151 
2152  return true;
2153  }
2154 }
2155 
2156 
2157 Foam::scalar Foam::primitiveMesh::setClosedThreshold(const scalar val)
2158 {
2159  scalar prev = closedThreshold_;
2160  closedThreshold_ = val;
2161 
2162  return prev;
2163 }
2164 
2165 
2166 Foam::scalar Foam::primitiveMesh::setAspectThreshold(const scalar val)
2167 {
2168  scalar prev = aspectThreshold_;
2169  aspectThreshold_ = val;
2170 
2171  return prev;
2172 }
2173 
2174 
2175 Foam::scalar Foam::primitiveMesh::setNonOrthThreshold(const scalar val)
2176 {
2177  scalar prev = nonOrthThreshold_;
2178  nonOrthThreshold_ = val;
2179 
2180  return prev;
2181 }
2182 
2183 
2184 Foam::scalar Foam::primitiveMesh::setSkewThreshold(const scalar val)
2185 {
2186  scalar prev = skewThreshold_;
2187  skewThreshold_ = val;
2188 
2189  return prev;
2190 }
2191 
2192 
2193 // ************************ vim: set sw=4 sts=4 et: ************************ //