Initial version of surface checks

meshLibrary/utilities/triSurfaceTools/triSurfaceChecks/triSurfaceChecks.C

@@ -281,13 +281,13 @@ label checkForHoles(const triSurf& surf, labelLongList& badTriangles)
     return badTriangles.size();
 }
 
-label checkForHoles(triSurf& surf, const word subsetPrefix)
+label checkForHoles(triSurf& surf, const word subsetName)
 {
     labelLongList trianglesNearHoles;
 
     if( checkForHoles(surf, trianglesNearHoles) )
     {
-        const label setId = surf.addFacetSubset(subsetPrefix);
+        const label setId = surf.addFacetSubset(subsetName);
 
         forAll(trianglesNearHoles, i)
             surf.addFacetToSubset(setId, trianglesNearHoles[i]);
@@ -623,6 +623,13 @@ label checkSelfIntersections
     {
         const labelledTri& tri = surf[tI];
 
+        const triangle<point, point> currTri
+        (
+            pts[tri[0]],
+            pts[tri[1]],
+            pts[tri[2]]
+        );
+
         boundBox bb(pts[tri[0]], pts[tri[0]]);
         for(label i=1;i<3;++i)
         {
@@ -649,21 +656,23 @@ label checkSelfIntersections
                 if( tI >= triJ )
                     continue;
 
-                point int0, int1;
-
-                const bool intersect =
-                    triangleFuncs::intersect
+                const triangle<point, point> neiTri
                 (
-                        pts[tri[0]],
-                        pts[tri[1]],
-                        pts[tri[2]],
-
                     pts[nt[0]],
                     pts[nt[1]],
-                        pts[nt[2]],
+                    pts[nt[2]]
+                );
+
+                FixedList<point, 2> intersectionEdge;
 
-                        int0,
-                        int1
+                const bool intersect =
+                    help::doTrianglesIntersect
+                    (
+                        currTri,
+                        neiTri,
+                        intersectionEdge,
+                        tol,
+                        Foam::cos(5.0 * M_PI / 180.0)
                     );
 
                 if( intersect )
@@ -728,6 +737,13 @@ label checkOverlaps
     {
         const labelledTri& tri = surf[tI];
 
+        const triangle<point, point> currTri
+        (
+            pts[tri[0]],
+            pts[tri[1]],
+            pts[tri[2]]
+        );
+
         boundBox bb(pts[tri[0]], pts[tri[0]]);
         for(label i=1;i<3;++i)
         {
@@ -754,21 +770,23 @@ label checkOverlaps
                 if( tI >= triJ )
                     continue;
 
-                point int0, int1;
-
-                const bool intersect =
-                    triangleFuncs::intersect
+                const triangle<point, point> neiTri
                 (
-                        pts[tri[0]],
-                        pts[tri[1]],
-                        pts[tri[2]],
-
                     pts[nt[0]],
                     pts[nt[1]],
-                        pts[nt[2]],
+                    pts[nt[2]]
+                );
 
-                        int0,
-                        int1
+                DynList<point> intersectionPolygon;
+
+                const bool intersect =
+                    help::doTrianglesOverlap
+                    (
+                        currTri,
+                        neiTri,
+                        intersectionPolygon,
+                        tol,
+                        Foam::cos(angleTol * M_PI / 180.0)
                     );
 
                 if( intersect )

meshLibrary/utilities/triSurfaceTools/triSurfaceChecks/triSurfaceChecks.H

@@ -74,7 +74,7 @@ label checkSurfaceManifolds(triSurf&, const word subsetPrefix="manifold_");
 
 //- check for existence of holes in the surface mesh
 label checkForHoles(const triSurf&, labelLongList&);
-label checkForHoles(triSurf&, const word subsetPrefix="hole_");
+label checkForHoles(triSurf&, const word subsetName="holes");
 
 //- check for existence of non-manifold edges
 label checkForNonManifoldEdges(const triSurf&, labelLongList&);

utilities/checkSurfaceMesh/checkSurfaceMesh.C

@@ -30,12 +30,8 @@ Description
 #include "IFstream.H"
 #include "fileName.H"
 #include "triSurf.H"
-#include "triSurfModifier.H"
-#include "helperFunctions.H"
-#include "demandDrivenData.H"
-#include "coordinateModifier.H"
-#include "checkMeshDict.H"
-#include "surfaceMeshGeometryModification.H"
+#include "triSurfaceChecks.H"
+#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
@@ -43,79 +39,106 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
-#   include "createTime.H"
-
-    IOdictionary meshDict
-    (
-        IOobject
-        (
-            "meshDict",
-            runTime.system(),
-            runTime,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    checkMeshDict cmd(meshDict);
-
-    fileName surfaceFile = meshDict.lookup("surfaceFile");
-    if( Pstream::parRun() )
-        surfaceFile = ".."/surfaceFile;
-
-    triSurf surface(runTime.path()/surfaceFile);
-
-    surfaceMeshGeometryModification gMod(surface, meshDict);
-
-    //- modify points
-    const triSurf* modSurfPtr = gMod.modifyGeometry();
-
-    Info << "Writting modified surface" << endl;
-    modSurfPtr->writeSurface("modifiedSurf.stl");
-
-    # ifdef DEBUGScaling
-    //- apply backward modification
-    Info << "Here" << endl;
-    coordinateModifier cMod(meshDict.subDict("anisotropicSources"));
-    Info << "Starting modifications" << endl;
-    forAll(surface.points(), i)
+    argList::noParallel();
+    argList::validArgs.clear();
+    argList::validArgs.append("input surface file");
+    argList args(argc, argv);
+
+    fileName inFileName(args.args()[1]);
+
+    triSurf surf(inFileName);
+
+    label nFailed(0);
+
+    boundBox bb;
+    triSurfaceChecks::calculateBoundingBox(surf, bb);
+    Info << "Surface bounding box is " << bb << endl;
+
+    //- calculate manifolds
+    const label nManifolds = triSurfaceChecks::checkSurfaceManifolds(surf);
+    if( nManifolds > 1 )
+    {
+        ++nFailed;
+
+        Info << "Surface mesh consists of " << nManifolds
+             << " manifolds." << endl;
+        Warning << "You cannot mesh geometries consisting of more than"
+                << " one domain, and it must not contain baffles." << endl;
+    }
+    else
+    {
+        Info << "Surface mesh consists of a single manifold." << endl;
+    }
+
+    //- find open boundary edges
+    if( triSurfaceChecks::checkForHoles(surf) )
+    {
+        ++nFailed;
+
+        Info << "Surface mesh has open boundaries!!" << endl;
+        Warning << "This indicates that there may be some holes in the surface"
+                << " mesh. Holes in the mesh must be smaller than the specified"
+                << " cell size at this location. In addition, please avoid"
+                << " using automatic refinement (minCellSize)." << endl;
+    }
+    else
     {
-        Info << "\nOrig point " << i << " coordinates " << surface.points()[i]
-             << " modified point " << modSurfPtr->points()[i] << endl;
-        const point p = cMod.backwardModifiedPoint(modSurfPtr->points()[i]);
+        Info << "No holes found in the surface mesh." << endl;
+    }
 
-        if( mag(p - surface.points()[i]) > 1e-14 )
+    //- find non-manifold edges
+    if( triSurfaceChecks::checkForNonManifoldEdges(surf) )
     {
-            Warning << "Point " << i << " is different "
-                    << p
-                    << " from original " << surface.points()[i]
-                    << " modified point "
-                    << cMod.modifiedPoint(surface.points()[i]) << endl;
-            ::exit(0);
+        ++nFailed;
+
+        Info << "Surface mesh has non-manifold edges!!" << endl;
+        Warning << "This indicates that the surface mesh consists of multiple"
+                << " domains and/or baffles. Please make sure that they are not"
+                << " in the domain which shall be meshed." << endl;
     }
+    else
+    {
+        Info << "Surface does not have any non-manifold edges." << endl;
     }
 
-    Info << "Backscaling Ok" << endl;
-    ::exit(0);
-    # endif
+    //- check the number of disconnected parts
+    if( triSurfaceChecks::checkDisconnectedParts(surf) > 1 )
+    {
+        ++nFailed;
+
+        Info << "Surface mesh consists of disconnected parts." << endl;
+        Warning << "This is not a problem if there exists a region surrounding"
+                << " the other ones! In other case, the mesher will generate"
+                << " the mesh in the domains with most cells." << endl;
+    }
+    else
+    {
+        Info << "Surface mesh consists of a single region." << endl;
+    }
 
-    surfaceMeshGeometryModification bgMod(*modSurfPtr, meshDict);
-    const triSurf* backModSurfPtr = bgMod.revertGeometryModification();
+    //- find triangles with small angles
+    if( triSurfaceChecks::checkAngles(surf, "smallAngles", 1.0) )
+    {
+        ++nFailed;
 
-    Info << "Writting backward transformed surface" << endl;
-    backModSurfPtr->writeSurface("backwardModifiedSurf.stl");
+        Info << "Surface mesh has some bad-quality triangles." << endl;
+        Warning << "This may cause problems to the automatic refinement"
+                << " procedure (minCellSize). " << endl;
+    }
+    else
+    {
+        Info << "No sliver triangles found." << endl;
+    }
 
-    # ifdef DEBUGScaling
-    forAll(backModSurfPtr->points(), pI)
-        if( mag(backModSurfPtr->points()[pI] - surface.points()[pI]) > 1e-14 )
-            Warning << "Point " << pI << " is different "
-                    << backModSurfPtr->points()[pI]
-                    << " from original " << surface.points()[pI] << endl;
-    # endif
+    if( nFailed )
+    {
+        Warning << "Found " << nFailed
+                << " checks indicating potential problems." << endl;
+        Warning << "This does not mean that you cannot generate"
+                << " a valid mesh. " << endl;
 
-    deleteDemandDrivenData(modSurfPtr);
-    deleteDemandDrivenData(backModSurfPtr);
+        surf.writeSurface(inFileName);
+    }
 
     Info << "End\n" << endl;