1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "primitiveMesh.H"
28 #include "pyramidPointFaceRef.H"
30 #include "mathematicalConstants.H"
31 #include "SortableList.H"
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 scalar primitiveMesh::closedThreshold_ = 1.0e-6;
41 scalar primitiveMesh::aspectThreshold_ = 1000;
42 scalar primitiveMesh::nonOrthThreshold_ = 70; // deg
43 scalar primitiveMesh::skewThreshold_ = 4;
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 bool primitiveMesh::checkClosedBoundary(const bool report) const
52 Info<< "bool primitiveMesh::checkClosedBoundary("
53 << "const bool) const: "
54 << "checking whether the boundary is closed" << endl;
57 // Loop through all boundary faces and sum up the face area vectors.
58 // For a closed boundary, this should be zero in all vector components
60 vector sumClosed(vector::zero);
61 scalar sumMagClosedBoundary = 0;
63 const vectorField& areas = faceAreas();
65 for (label faceI = nInternalFaces(); faceI < areas.size(); faceI++)
67 sumClosed += areas[faceI];
68 sumMagClosedBoundary += mag(areas[faceI]);
71 reduce(sumClosed, sumOp<vector>());
72 reduce(sumMagClosedBoundary, sumOp<scalar>());
74 vector openness = sumClosed/(sumMagClosedBoundary + VSMALL);
76 if (cmptMax(cmptMag(openness)) > closedThreshold_)
80 Info<< " ***Boundary openness " << openness
81 << " possible hole in boundary description."
91 Info<< " Boundary openness " << openness << " OK."
100 bool primitiveMesh::checkClosedCells
103 labelHashSet* setPtr,
104 labelHashSet* aspectSetPtr
109 Info<< "bool primitiveMesh::checkClosedCells("
110 << "const bool, labelHashSet*, labelHashSet*) const: "
111 << "checking whether cells are closed" << endl;
114 // Check that all cells labels are valid
115 const cellList& c = cells();
117 label nErrorClosed = 0;
121 const cell& curCell = c[cI];
123 if (min(curCell) < 0 || max(curCell) > nFaces())
134 if (nErrorClosed > 0)
138 Info<< " ***Cells with invalid face labels found, number of cells "
139 << nErrorClosed << endl;
145 // Loop through cell faces and sum up the face area vectors for each cell.
146 // This should be zero in all vector components
148 vectorField sumClosed(nCells(), vector::zero);
149 vectorField sumMagClosed(nCells(), vector::zero);
151 const labelList& own = faceOwner();
152 const labelList& nei = faceNeighbour();
153 const vectorField& areas = faceAreas();
158 sumClosed[own[faceI]] += areas[faceI];
159 sumMagClosed[own[faceI]] += cmptMag(areas[faceI]);
164 // Subtract from neighbour
165 sumClosed[nei[faceI]] -= areas[faceI];
166 sumMagClosed[nei[faceI]] += cmptMag(areas[faceI]);
170 scalar maxOpennessCell = 0;
173 scalar maxAspectRatio = 0;
175 const scalarField& vols = cellVolumes();
178 forAll (sumClosed, cellI)
180 scalar maxOpenness = 0;
182 for(direction cmpt=0; cmpt<vector::nComponents; cmpt++)
187 mag(sumClosed[cellI][cmpt])
188 /(sumMagClosed[cellI][cmpt] + VSMALL)
192 maxOpennessCell = max(maxOpennessCell, maxOpenness);
194 if (maxOpenness > closedThreshold_)
198 setPtr->insert(cellI);
204 // Calculate the aspect ration as the maximum of Cartesian component
205 // aspect ratio to the total area hydraulic area aspect ratio
206 scalar aspectRatio = max
208 cmptMax(sumMagClosed[cellI])
209 /(cmptMin(sumMagClosed[cellI]) + VSMALL),
210 1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
213 maxAspectRatio = max(maxAspectRatio, aspectRatio);
215 if (aspectRatio > aspectThreshold_)
219 aspectSetPtr->insert(cellI);
226 reduce(nOpen, sumOp<label>());
227 reduce(maxOpennessCell, maxOp<scalar>());
229 reduce(nAspect, sumOp<label>());
230 reduce(maxAspectRatio, maxOp<scalar>());
237 Info<< " ***Open cells found, max cell openness: "
238 << maxOpennessCell << ", number of open cells " << nOpen
249 Info<< " ***High aspect ratio cells found, Max aspect ratio: "
251 << ", number of cells " << nAspect
260 Info<< " Max cell openness = " << maxOpennessCell << " OK." << nl
261 << " Max aspect ratio = " << maxAspectRatio << " OK."
269 bool primitiveMesh::checkFaceAreas
277 Info<< "bool primitiveMesh::checkFaceAreas("
278 << "const bool, labelHashSet*) const: "
279 << "checking face area magnitudes" << endl;
282 const scalarField magFaceAreas = mag(faceAreas());
284 scalar minArea = GREAT;
285 scalar maxArea = -GREAT;
287 forAll (magFaceAreas, faceI)
289 if (magFaceAreas[faceI] < VSMALL)
293 setPtr->insert(faceI);
297 minArea = min(minArea, magFaceAreas[faceI]);
298 maxArea = max(maxArea, magFaceAreas[faceI]);
301 reduce(minArea, minOp<scalar>());
302 reduce(maxArea, maxOp<scalar>());
304 if (minArea < VSMALL)
308 Info<< " ***Zero or negative face area detected. "
309 "Minimum area: " << minArea << endl;
318 Info<< " Minumum face area = " << minArea
319 << ". Maximum face area = " << maxArea
320 << ". Face area magnitudes OK." << endl;
328 bool primitiveMesh::checkCellVolumes
336 Info<< "bool primitiveMesh::checkCellVolumes("
337 << "const bool, labelHashSet*) const: "
338 << "checking cell volumes" << endl;
341 const scalarField& vols = cellVolumes();
343 scalar minVolume = GREAT;
344 scalar maxVolume = -GREAT;
346 label nNegVolCells = 0;
350 if (vols[cellI] < VSMALL)
354 setPtr->insert(cellI);
360 minVolume = min(minVolume, vols[cellI]);
361 maxVolume = max(maxVolume, vols[cellI]);
364 reduce(minVolume, minOp<scalar>());
365 reduce(maxVolume, maxOp<scalar>());
366 reduce(nNegVolCells, sumOp<label>());
368 if (minVolume < VSMALL)
372 Info<< " ***Zero or negative cell volume detected. "
373 << "Minimum negative volume: " << minVolume
374 << ", Number of negative volume cells: " << nNegVolCells
384 Info<< " Min volume = " << minVolume
385 << ". Max volume = " << maxVolume
386 << ". Total volume = " << gSum(vols)
387 << ". Cell volumes OK." << endl;
395 bool primitiveMesh::checkFaceOrthogonality
403 Info<< "bool primitiveMesh::checkFaceOrthogonality("
404 << "const bool, labelHashSet*) const: "
405 << "checking mesh non-orthogonality" << endl;
408 // for all internal faces check theat the d dot S product is positive
409 const vectorField& centres = cellCentres();
410 const vectorField& areas = faceAreas();
412 const labelList& own = faceOwner();
413 const labelList& nei = faceNeighbour();
415 // Severe nonorthogonality threshold
416 const scalar severeNonorthogonalityThreshold =
417 ::cos(nonOrthThreshold_/180.0*mathematicalConstant::pi);
419 scalar minDDotS = GREAT;
423 label severeNonOrth = 0;
425 label errorNonOrth = 0;
429 vector d = centres[nei[faceI]] - centres[own[faceI]];
430 const vector& s = areas[faceI];
432 scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
434 if (dDotS < severeNonorthogonalityThreshold)
440 setPtr->insert(faceI);
449 setPtr->insert(faceI);
456 if (dDotS < minDDotS)
464 reduce(minDDotS, minOp<scalar>());
465 reduce(sumDDotS, sumOp<scalar>());
466 reduce(severeNonOrth, sumOp<label>());
467 reduce(errorNonOrth, sumOp<label>());
471 label neiSize = nei.size();
472 reduce(neiSize, sumOp<label>());
478 Info<< " Mesh non-orthogonality Max: "
479 << ::acos(minDDotS)/mathematicalConstant::pi*180.0
481 ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0
486 if (severeNonOrth > 0)
488 Info<< " *Number of severely non-orthogonal faces: "
489 << severeNonOrth << "." << endl;
493 if (errorNonOrth > 0)
497 Info<< " ***Number of non-orthogonality errors: "
498 << errorNonOrth << "." << endl;
507 Info<< " Non-orthogonality check OK." << endl;
515 bool primitiveMesh::checkFacePyramids
518 const scalar minPyrVol,
524 Info<< "bool primitiveMesh::checkFacePyramids("
525 << "const bool, const scalar, labelHashSet*) const: "
526 << "checking face orientation" << endl;
529 // check whether face area vector points to the cell with higher label
530 const vectorField& ctrs = cellCentres();
532 const labelList& own = faceOwner();
533 const labelList& nei = faceNeighbour();
535 const faceList& f = faces();
537 const pointField& p = points();
539 label nErrorPyrs = 0;
543 // Create the owner pyramid - it will have negative volume
544 scalar pyrVol = pyramidPointFaceRef(f[faceI], ctrs[own[faceI]]).mag(p);
546 if (pyrVol > -minPyrVol)
550 setPtr->insert(faceI);
556 if (isInternalFace(faceI))
558 // Create the neighbour pyramid - it will have positive volume
560 pyramidPointFaceRef(f[faceI], ctrs[nei[faceI]]).mag(p);
562 if (pyrVol < minPyrVol)
566 setPtr->insert(faceI);
574 reduce(nErrorPyrs, sumOp<label>());
580 Info<< " ***Error in face pyramids: "
581 << nErrorPyrs << " faces are incorrectly oriented."
591 Info<< " Face pyramids OK." << endl;
599 bool primitiveMesh::checkFaceSkewness
607 Info<< "bool primitiveMesh::checkFaceSkewnesss("
608 << "const bool, labelHashSet*) const: "
609 << "checking face skewness" << endl;
612 // Warn if the skew correction vector is more than skewWarning times
613 // larger than the face area vector
615 const pointField& p = points();
616 const faceList& fcs = faces();
618 const labelList& own = faceOwner();
619 const labelList& nei = faceNeighbour();
620 const vectorField& cellCtrs = cellCentres();
621 const vectorField& faceCtrs = faceCentres();
622 const vectorField& fAreas = faceAreas();
629 vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
630 vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
633 vector sv = Cpf - ((fAreas[faceI] & Cpf)/(fAreas[faceI] & d))*d;
634 vector svHat = sv/(mag(sv) + VSMALL);
636 // Normalisation distance calculated as the approximate distance
637 // from the face centre to the edge of the face in the direction of
639 scalar fd = 0.2*mag(d) + VSMALL;
640 const face& f = fcs[faceI];
643 fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
646 // Normalised skewness
647 scalar skewness = mag(sv)/fd;
649 // Check if the skewness vector is greater than the PN vector.
650 // This does not cause trouble but is a good indication of a poor mesh.
651 if (skewness > skewThreshold_)
655 setPtr->insert(faceI);
661 if(skewness > maxSkew)
668 // Boundary faces: consider them to have only skewness error.
669 // (i.e. treat as if mirror cell on other side)
671 for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
673 vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
675 vector normal = fAreas[faceI];
676 normal /= mag(normal) + VSMALL;
677 vector d = normal*(normal & Cpf);
681 vector sv = Cpf - ((fAreas[faceI]&Cpf)/((fAreas[faceI]&d)+VSMALL))*d;
682 vector svHat = sv/(mag(sv) + VSMALL);
684 // Normalisation distance calculated as the approximate distance
685 // from the face centre to the edge of the face in the direction of
687 scalar fd = 0.4*mag(d) + VSMALL;
688 const face& f = fcs[faceI];
691 fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
694 // Normalised skewness
695 scalar skewness = mag(sv)/fd;
697 // Check if the skewness vector is greater than the PN vector.
698 // This does not cause trouble but is a good indication of a poor mesh.
699 if (skewness > skewThreshold_)
703 setPtr->insert(faceI);
709 if(skewness > maxSkew)
716 reduce(maxSkew, maxOp<scalar>());
717 reduce(nWarnSkew, sumOp<label>());
723 Info<< " ***Max skewness = " << maxSkew
724 << ", " << nWarnSkew << " highly skew faces detected"
725 " which may impair the quality of the results"
735 Info<< " Max skewness = " << maxSkew << " OK." << endl;
743 bool primitiveMesh::checkPoints
751 Info<< "bool primitiveMesh::checkPoints"
752 << "(const bool, labelHashSet*) const: "
753 << "checking points" << endl;
756 label nFaceErrors = 0;
757 label nCellErrors = 0;
759 const labelListList& pf = pointFaces();
763 if (pf[pointI].size() == 0)
767 setPtr->insert(pointI);
774 const labelListList& pc = pointCells();
778 if (pc[pointI].size() == 0)
782 setPtr->insert(pointI);
789 reduce(nFaceErrors, sumOp<label>());
790 reduce(nCellErrors, sumOp<label>());
792 if (nFaceErrors > 0 || nCellErrors > 0)
796 Info<< " ***Unsed points found in the mesh, "
797 "number unused by faces: " << nFaceErrors
798 << " number unused by cells: " << nCellErrors
808 Info<< " Point usage OK." << endl;
816 // Check convexity of angles in a face. Allow a slight non-convexity.
817 // E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
818 // (if truly concave and points not visible from face centre the face-pyramid
819 // check in checkMesh will fail)
820 bool primitiveMesh::checkFaceAngles
829 Info<< "bool primitiveMesh::checkFaceAngles"
830 << "(const bool, const scalar, labelHashSet*) const: "
831 << "checking face angles" << endl;
834 if (maxDeg < -SMALL || maxDeg > 180+SMALL)
838 "primitiveMesh::checkFaceAngles"
839 "(const bool, const scalar, labelHashSet*)"
840 ) << "maxDeg should be [0..180] but is now " << maxDeg
844 const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
846 const pointField& p = points();
847 const faceList& fcs = faces();
848 vectorField faceNormals(faceAreas());
849 faceNormals /= mag(faceNormals) + VSMALL;
851 scalar maxEdgeSin = 0.0;
855 label errorFaceI = -1;
859 const face& f = fcs[faceI];
861 // Get edge from f[0] to f[size-1];
862 vector ePrev(p[f[0]] - p[f[f.size()-1]]);
863 scalar magEPrev = mag(ePrev);
864 ePrev /= magEPrev + VSMALL;
868 // Get vertex after fp
869 label fp1 = (fp0 + 1) % f.size();
871 // Normalized vector between two consecutive points
872 vector e10(p[f[fp1]] - p[f[fp0]]);
873 scalar magE10 = mag(e10);
874 e10 /= magE10 + VSMALL;
876 if (magEPrev > SMALL && magE10 > SMALL)
878 vector edgeNormal = ePrev ^ e10;
879 scalar magEdgeNormal = mag(edgeNormal);
881 if (magEdgeNormal < maxSin)
883 // Edges (almost) aligned -> face is ok.
888 edgeNormal /= magEdgeNormal;
890 if ((edgeNormal & faceNormals[faceI]) < SMALL)
892 if (faceI != errorFaceI)
894 // Count only one error per face.
901 setPtr->insert(faceI);
904 maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
914 reduce(nConcave, sumOp<label>());
915 reduce(maxEdgeSin, maxOp<scalar>());
919 scalar maxConcaveDegr =
920 Foam::asin(Foam::min(1.0, maxEdgeSin))
921 *180.0/mathematicalConstant::pi;
925 Info<< " *There are " << nConcave
926 << " faces with concave angles between consecutive"
927 << " edges. Max concave angle = " << maxConcaveDegr
928 << " degrees." << endl;
937 Info<< " All angles in faces OK." << endl;
945 // Check warpage of faces. Is calculated as the difference between areas of
946 // individual triangles and the overall area of the face (which ifself is
947 // is the average of the areas of the individual triangles).
948 bool primitiveMesh::checkFaceFlatness
951 const scalar warnFlatness,
957 Info<< "bool primitiveMesh::checkFaceFlatness"
958 << "(const bool, const scalar, labelHashSet*) const: "
959 << "checking face flatness" << endl;
962 if (warnFlatness < 0 || warnFlatness > 1)
966 "primitiveMesh::checkFaceFlatness"
967 "(const bool, const scalar, labelHashSet*)"
968 ) << "warnFlatness should be [0..1] but is now " << warnFlatness
973 const pointField& p = points();
974 const faceList& fcs = faces();
975 const pointField& fctrs = faceCentres();
977 // Areas are calculated as the sum of areas. (see
978 // primitiveMeshFaceCentresAndAreas.C)
979 scalarField magAreas(mag(faceAreas()));
983 scalar minFlatness = GREAT;
984 scalar sumFlatness = 0;
989 const face& f = fcs[faceI];
991 if (f.size() > 3 && magAreas[faceI] > VSMALL)
993 const point& fc = fctrs[faceI];
995 // Calculate the sum of magnitude of areas and compare to magnitude
1002 const point& thisPoint = p[f[fp]];
1003 const point& nextPoint = p[f.nextLabel(fp)];
1005 // Triangle around fc.
1006 vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
1010 scalar flatness = magAreas[faceI] / (sumA+VSMALL);
1012 sumFlatness += flatness;
1015 minFlatness = min(minFlatness, flatness);
1017 if (flatness < warnFlatness)
1023 setPtr->insert(faceI);
1030 reduce(nWarped, sumOp<label>());
1031 reduce(minFlatness, minOp<scalar>());
1033 reduce(nSummed, sumOp<label>());
1034 reduce(sumFlatness, sumOp<scalar>());
1036 if (debug || report)
1040 Info<< " Face flatness (1 = flat, 0 = butterfly) : average = "
1041 << sumFlatness / nSummed << " min = " << minFlatness << endl;
1048 if (debug || report)
1050 Info<< " *There are " << nWarped
1051 << " faces with ratio between projected and actual area < "
1052 << warnFlatness << endl;
1054 Info<< " Minimum ratio (minimum flatness, maximum warpage) = "
1055 << minFlatness << endl;
1062 if (debug || report)
1064 Info<< " All face flatness OK." << endl;
1072 // Check 1D/2Dness of edges. Gets passed the non-empty directions and
1073 // checks all edges in the mesh whether they:
1074 // - have no component in a non-empty direction or
1075 // - are only in a singe non-empty direction.
1076 // Empty direction info is passed in as a vector of labels (synchronised)
1077 // which are 1 if the direction is non-empty, 0 if it is.
1078 bool primitiveMesh::checkEdgeAlignment
1081 const Vector<label>& directions,
1082 labelHashSet* setPtr
1087 Info<< "bool primitiveMesh::checkEdgeAlignment("
1088 << "const bool, const Vector<label>&, labelHashSet*) const: "
1089 << "checking edge alignment" << endl;
1093 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
1095 if (directions[cmpt] == 1)
1099 else if (directions[cmpt] != 0)
1103 "primitiveMesh::checkEdgeAlignment"
1104 "(const bool, const Vector<label>&, labelHashSet*)"
1105 ) << "directions should contain 0 or 1 but is now " << directions
1106 << exit(FatalError);
1110 if (nDirs == vector::nComponents)
1117 const pointField& p = points();
1118 const faceList& fcs = faces();
1120 EdgeMap<label> edgesInError;
1124 const face& f = fcs[faceI];
1129 label p1 = f.nextLabel(fp);
1132 vector d(p[p1]-p[p0]);
1133 scalar magD = mag(d);
1135 if (magD > ROOTVSMALL)
1139 // Check how many empty directions are used by the edge.
1140 label nEmptyDirs = 0;
1141 label nNonEmptyDirs = 0;
1142 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
1144 if (mag(d[cmpt]) > 1e-6)
1146 if (directions[cmpt] == 0)
1157 if (nEmptyDirs == 0)
1159 // Purely in ok directions.
1161 else if (nEmptyDirs == 1)
1163 // Ok if purely in empty directions.
1164 if (nNonEmptyDirs > 0)
1166 edgesInError.insert(edge(p0, p1), faceI);
1169 else if (nEmptyDirs > 1)
1172 edgesInError.insert(edge(p0, p1), faceI);
1179 label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
1181 if (nErrorEdges > 0)
1183 if (debug || report)
1185 Info<< " ***Number of edges not aligned with or perpendicular to "
1186 << "non-empty directions: " << nErrorEdges << endl;
1191 setPtr->resize(2*edgesInError.size());
1192 forAllConstIter(EdgeMap<label>, edgesInError, iter)
1194 setPtr->insert(iter.key()[0]);
1195 setPtr->insert(iter.key()[1]);
1203 if (debug || report)
1205 Info<< " All edges aligned with or perpendicular to "
1206 << "non-empty directions." << endl;
1213 bool primitiveMesh::checkUpperTriangular
1216 labelHashSet* setPtr
1221 Info<< "bool primitiveMesh::checkUpperTriangular("
1222 << "const bool, labelHashSet*) const: "
1223 << "checking face ordering" << endl;
1226 // Check whether internal faces are ordered in the upper triangular order
1227 const labelList& own = faceOwner();
1228 const labelList& nei = faceNeighbour();
1230 const cellList& c = cells();
1232 label internal = nInternalFaces();
1234 // Has error occurred?
1236 // Have multiple faces been detected?
1237 label nMultipleCells = false;
1239 // Loop through faceCells once more and make sure that for internal cell
1240 // the first label is smaller
1241 for (label faceI = 0; faceI < internal; faceI++)
1243 if (own[faceI] >= nei[faceI])
1249 setPtr->insert(faceI);
1254 // Loop through all cells. For each cell, find the face that is internal
1255 // and add it to the check list (upper triangular order).
1256 // Once the list is completed, check it against the faceCell list
1260 const labelList& curFaces = c[cellI];
1262 // Neighbouring cells
1263 SortableList<label> nbr(curFaces.size());
1267 label faceI = curFaces[i];
1269 if (faceI >= nInternalFaces())
1276 label nbrCellI = nei[faceI];
1278 if (nbrCellI == cellI)
1280 nbrCellI = own[faceI];
1283 if (cellI < nbrCellI)
1290 // nbrCell is master. Let it handle this face.
1298 // Now nbr holds the cellCells in incremental order. Check:
1299 // - neighbouring cells appear only once. Since nbr is sorted this
1300 // is simple check on consecutive elements
1301 // - faces indexed in same order as nbr are incrementing as well.
1303 label prevCell = nbr[0];
1304 label prevFace = curFaces[nbr.indices()[0]];
1306 bool hasMultipleFaces = false;
1308 for (label i = 1; i < nbr.size(); i++)
1310 label thisCell = nbr[i];
1311 label thisFace = curFaces[nbr.indices()[i]];
1313 if (thisCell == labelMax)
1318 if (thisCell == prevCell)
1320 hasMultipleFaces = true;
1324 setPtr->insert(prevFace);
1325 setPtr->insert(thisFace);
1328 else if (thisFace < prevFace)
1334 setPtr->insert(thisFace);
1338 prevCell = thisCell;
1339 prevFace = thisFace;
1342 if (hasMultipleFaces)
1348 reduce(error, orOp<bool>());
1349 reduce(nMultipleCells, sumOp<label>());
1351 if ((debug || report) && nMultipleCells > 0)
1353 Info<< " <<Found " << nMultipleCells
1354 << " neighbouring cells with multiple inbetween faces." << endl;
1359 if (debug || report)
1361 Info<< " ***Faces not in upper triangular order." << endl;
1368 if (debug || report)
1370 Info<< " Upper triangular ordering OK." << endl;
1378 bool primitiveMesh::checkCellsZipUp
1381 labelHashSet* setPtr
1386 Info<< "bool primitiveMesh::checkCellsZipUp("
1387 << "const bool, labelHashSet*) const: "
1388 << "checking topological cell openness" << endl;
1391 label nOpenCells = 0;
1393 const faceList& f = faces();
1394 const cellList& c = cells();
1398 const labelList& curFaces = c[cellI];
1400 const edgeList cellEdges = c[cellI].edges(f);
1402 labelList edgeUsage(cellEdges.size(), 0);
1404 forAll (curFaces, faceI)
1406 edgeList curFaceEdges = f[curFaces[faceI]].edges();
1408 forAll (curFaceEdges, faceEdgeI)
1410 const edge& curEdge = curFaceEdges[faceEdgeI];
1412 forAll (cellEdges, cellEdgeI)
1414 if (cellEdges[cellEdgeI] == curEdge)
1416 edgeUsage[cellEdgeI]++;
1423 edgeList singleEdges(cellEdges.size());
1424 label nSingleEdges = 0;
1426 forAll (edgeUsage, edgeI)
1428 if (edgeUsage[edgeI] == 1)
1430 singleEdges[nSingleEdges] = cellEdges[edgeI];
1433 else if (edgeUsage[edgeI] != 2)
1437 setPtr->insert(cellI);
1442 if (nSingleEdges > 0)
1446 setPtr->insert(cellI);
1453 reduce(nOpenCells, sumOp<label>());
1457 if (debug || report)
1459 Info<< " ***Open cells found, number of cells: " << nOpenCells
1460 << ". This problem may be fixable using the zipUpMesh utility."
1468 if (debug || report)
1470 Info<< " Topological cell zip-up check OK." << endl;
1478 // Vertices of face within point range and unique.
1479 bool primitiveMesh::checkFaceVertices
1482 labelHashSet* setPtr
1487 Info<< "bool primitiveMesh::checkFaceVertices("
1488 << "const bool, labelHashSet*) const: "
1489 << "checking face vertices" << endl;
1492 // Check that all vertex labels are valid
1493 const faceList& f = faces();
1495 label nErrorFaces = 0;
1499 const face& curFace = f[fI];
1501 if (min(curFace) < 0 || max(curFace) > nPoints())
1511 // Uniqueness of vertices
1512 labelHashSet facePoints(2*curFace.size());
1516 bool inserted = facePoints.insert(curFace[fp]);
1530 reduce(nErrorFaces, sumOp<label>());
1532 if (nErrorFaces > 0)
1534 if (debug || report)
1536 Info<< " Faces with invalid vertex labels found, "
1537 << " number of faces: " << nErrorFaces << endl;
1544 if (debug || report)
1546 Info<< " Face vertices OK." << endl;
1554 // Check if all points on face are shared between faces.
1555 bool primitiveMesh::checkDuplicateFaces
1558 const Map<label>& nCommonPoints,
1559 label& nBaffleFaces,
1560 labelHashSet* setPtr
1565 forAllConstIter(Map<label>, nCommonPoints, iter)
1567 label nbFaceI = iter.key();
1568 label nCommon = iter();
1570 const face& curFace = faces()[faceI];
1571 const face& nbFace = faces()[nbFaceI];
1573 if (nCommon == nbFace.size() || nCommon == curFace.size())
1575 if (nbFace.size() != curFace.size())
1586 setPtr->insert(faceI);
1587 setPtr->insert(nbFaceI);
1596 // Check that shared points are in consecutive order.
1597 bool primitiveMesh::checkCommonOrder
1600 const Map<label>& nCommonPoints,
1601 labelHashSet* setPtr
1606 forAllConstIter(Map<label>, nCommonPoints, iter)
1608 label nbFaceI = iter.key();
1609 label nCommon = iter();
1611 const face& curFace = faces()[faceI];
1612 const face& nbFace = faces()[nbFaceI];
1617 && nCommon != nbFace.size()
1618 && nCommon != curFace.size()
1623 // Get the index in the neighbouring face shared with curFace
1624 label nb = findIndex(nbFace, curFace[fp]);
1629 // Check the whole face from nb onwards for shared vertices
1630 // with neighbouring face. Rule is that any shared vertices
1631 // should be consecutive on both faces i.e. if they are
1632 // vertices fp,fp+1,fp+2 on one face they should be
1633 // vertices nb, nb+1, nb+2 (or nb+2, nb+1, nb) on the
1637 // Vertices before and after on curFace
1638 label fpPlus1 = (fp+1) % curFace.size();
1639 label fpMin1 = (fp == 0 ? curFace.size()-1 : fp-1);
1641 // Vertices before and after on nbFace
1642 label nbPlus1 = (nb+1) % nbFace.size();
1643 label nbMin1 = (nb == 0 ? nbFace.size()-1 : nb-1);
1645 // Find order of walking by comparing next points on both
1647 label curInc = labelMax;
1648 label nbInc = labelMax;
1650 if (nbFace[nbPlus1] == curFace[fpPlus1])
1655 else if (nbFace[nbPlus1] == curFace[fpMin1])
1660 else if (nbFace[nbMin1] == curFace[fpMin1])
1672 // Pass1: loop until start of common vertices found.
1680 if (curFp >= curFace.size())
1686 curFp = curFace.size()-1;
1691 if (curNb >= nbFace.size())
1697 curNb = nbFace.size()-1;
1699 } while (curFace[curFp] == nbFace[curNb]);
1702 // Pass2: check equality walking from curFp, curNb
1703 // in opposite order.
1708 for (label commonI = 0; commonI < nCommon; commonI++)
1712 if (curFp >= curFace.size())
1718 curFp = curFace.size()-1;
1723 if (curNb >= nbFace.size())
1729 curNb = nbFace.size()-1;
1732 if (curFace[curFp] != nbFace[curNb])
1736 setPtr->insert(faceI);
1737 setPtr->insert(nbFaceI);
1747 // Done the curFace - nbFace combination.
1758 // Checks common vertices between faces. If more than 2 they should be
1759 // consecutive on both faces.
1760 bool primitiveMesh::checkFaceFaces
1763 labelHashSet* setPtr
1768 Info<< "bool primitiveMesh::checkFaceFaces(const bool, labelHashSet*)"
1769 << " const: " << "checking face-face connectivity" << endl;
1772 const labelListList& pf = pointFaces();
1774 label nBaffleFaces = 0;
1775 label nErrorDuplicate = 0;
1776 label nErrorOrder = 0;
1777 Map<label> nCommonPoints(100);
1779 for (label faceI = 0; faceI < nFaces(); faceI++)
1781 const face& curFace = faces()[faceI];
1783 // Calculate number of common points between current faceI and
1784 // neighbouring face. Store on map.
1785 nCommonPoints.clear();
1789 label pointI = curFace[fp];
1791 const labelList& nbs = pf[pointI];
1795 label nbFaceI = nbs[nbI];
1797 if (faceI < nbFaceI)
1799 // Only check once for each combination of two faces.
1801 Map<label>::iterator fnd = nCommonPoints.find(nbFaceI);
1803 if (fnd == nCommonPoints.end())
1805 // First common vertex found.
1806 nCommonPoints.insert(nbFaceI, 1);
1816 // Perform various checks on common points
1818 // Check all vertices shared (duplicate point)
1819 if (checkDuplicateFaces(faceI, nCommonPoints, nBaffleFaces, setPtr))
1824 // Check common vertices are consecutive on both faces
1825 if (checkCommonOrder(faceI, nCommonPoints, setPtr))
1831 reduce(nBaffleFaces, sumOp<label>());
1832 reduce(nErrorDuplicate, sumOp<label>());
1833 reduce(nErrorOrder, sumOp<label>());
1837 Info<< " Number of identical duplicate faces (baffle faces): "
1838 << nBaffleFaces << endl;
1841 if (nErrorDuplicate > 0 || nErrorOrder > 0)
1843 if (nErrorDuplicate > 0)
1845 Info<< " ***Number of duplicate (not baffle) faces found: "
1846 << nErrorDuplicate << endl;
1849 if (nErrorOrder > 0)
1851 Info<< " ***Number of faces with non-consecutive shared points: "
1852 << nErrorOrder << endl;
1859 if (debug || report)
1861 Info<< " Face-face connectivity OK." << endl;
1869 // Checks cells with 1 or less internal faces. Give numerical problems.
1870 bool primitiveMesh::checkCellDeterminant
1872 const bool report, // report,
1873 labelHashSet* setPtr // setPtr
1878 Info<< "bool primitiveMesh::checkCellDeterminant(const bool"
1879 << ", labelHashSet*) const: "
1880 << "checking for under-determined cells" << endl;
1883 const cellList& c = cells();
1885 label nErrorCells = 0;
1887 scalar minDet = GREAT;
1893 const labelList& curFaces = c[cellI];
1895 // Calculate local normalization factor
1898 label nInternalFaces = 0;
1902 if (isInternalFace(curFaces[i]))
1904 avgArea += mag(faceAreas()[curFaces[i]]);
1910 if (nInternalFaces == 0)
1914 setPtr->insert(cellI);
1921 avgArea /= nInternalFaces;
1923 symmTensor areaTensor(symmTensor::zero);
1927 if (isInternalFace(curFaces[i]))
1929 areaTensor += sqr(faceAreas()[curFaces[i]]/avgArea);
1933 scalar determinant = mag(det(areaTensor));
1935 minDet = min(determinant, minDet);
1936 sumDet += determinant;
1939 if (determinant < 1e-3)
1943 setPtr->insert(cellI);
1951 reduce(nErrorCells, sumOp<label>());
1952 reduce(minDet, minOp<scalar>());
1953 reduce(sumDet, sumOp<scalar>());
1954 reduce(nSummed, sumOp<label>());
1956 if (debug || report)
1960 Info<< " Cell determinant (wellposedness) : minimum: " << minDet
1961 << " average: " << sumDet/nSummed
1966 if (nErrorCells > 0)
1968 if (debug || report)
1970 Info<< " ***Cells with small determinant found, number of cells: "
1971 << nErrorCells << endl;
1978 if (debug || report)
1980 Info<< " Cell determinant check OK." << endl;
1990 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1992 bool primitiveMesh::checkTopology(const bool report) const
1994 label noFailedChecks = 0;
1996 if (checkPoints(report)) noFailedChecks++;
1997 if (checkUpperTriangular(report)) noFailedChecks++;
1998 if (checkCellsZipUp(report)) noFailedChecks++;
1999 if (checkFaceVertices(report)) noFailedChecks++;
2000 if (checkFaceFaces(report)) noFailedChecks++;
2001 //if (checkCellDeterminant(report)) noFailedChecks++;
2003 if (noFailedChecks == 0)
2005 if (debug || report)
2007 Info<< " Mesh topology OK." << endl;
2014 if (debug || report)
2016 Info<< " Failed " << noFailedChecks
2017 << " mesh topology checks." << endl;
2025 bool primitiveMesh::checkGeometry(const bool report) const
2027 label noFailedChecks = 0;
2029 if (checkClosedBoundary(report)) noFailedChecks++;
2030 if (checkClosedCells(report)) noFailedChecks++;
2031 if (checkFaceAreas(report)) noFailedChecks++;
2032 if (checkCellVolumes(report)) noFailedChecks++;
2033 if (checkFaceOrthogonality(report)) noFailedChecks++;
2034 if (checkFacePyramids(report)) noFailedChecks++;
2035 if (checkFaceSkewness(report)) noFailedChecks++;
2037 if (noFailedChecks == 0)
2039 if (debug || report)
2041 Info<< " Mesh geometry OK." << endl;
2047 if (debug || report)
2049 Info<< " Failed " << noFailedChecks
2050 << " mesh geometry checks." << endl;
2058 bool primitiveMesh::checkMesh(const bool report) const
2062 Info<< "bool primitiveMesh::checkMesh(const bool report) const: "
2063 << "checking primitiveMesh" << endl;
2066 label noFailedChecks = checkTopology(report) + checkGeometry(report);
2068 if (noFailedChecks == 0)
2070 if (debug || report)
2072 Info<< "Mesh OK." << endl;
2079 if (debug || report)
2081 Info<< " Failed " << noFailedChecks
2082 << " mesh checks." << endl;
2090 scalar primitiveMesh::setClosedThreshold(const scalar ct)
2092 scalar oldClosedThreshold = closedThreshold_;
2093 closedThreshold_ = ct;
2095 return oldClosedThreshold;
2099 scalar primitiveMesh::setAspectThreshold(const scalar at)
2101 scalar oldAspectThreshold = aspectThreshold_;
2102 aspectThreshold_ = at;
2104 return oldAspectThreshold;
2108 scalar primitiveMesh::setNonOrthThreshold(const scalar nt)
2110 scalar oldNonOrthThreshold = nonOrthThreshold_;
2111 nonOrthThreshold_ = nt;
2113 return oldNonOrthThreshold;
2117 scalar primitiveMesh::setSkewThreshold(const scalar st)
2119 scalar oldSkewThreshold = skewThreshold_;
2120 skewThreshold_ = st;
2122 return oldSkewThreshold;
2126 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
2128 } // End namespace Foam
2130 // ************************************************************************* //