Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -35,9 +35,78 @@ License
 #include "processorPolyPatch.H"
 #include "surfaceWriter.H"
 #include "checkTools.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+template<class PatchType>
+void Foam::checkPatch
+(
+    const bool allGeometry,
+    const word& name,
+    const PatchType& pp,
+    pointSet& points
+)
+{
+    Info<< "    "
+        << setw(20) << name
+        << setw(9) << returnReduce(pp.size(), sumOp<label>())
+        << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
+
+    if (!Pstream::parRun())
+    {
+        typedef typename PatchType::surfaceTopo TopoType;
+        TopoType pTyp = pp.surfaceType();
+
+        if (pp.empty())
+        {
+            Info<< setw(34) << "ok (empty)";
+        }
+        else if (pTyp == TopoType::MANIFOLD)
+        {
+            if (pp.checkPointManifold(true, &points))
+            {
+                Info<< setw(34)
+                    << "multiply connected (shared point)";
+            }
+            else
+            {
+                Info<< setw(34) << "ok (closed singly connected)";
+            }
+
+            // Add points on non-manifold edges to make set complete
+            pp.checkTopology(false, &points);
+        }
+        else
+        {
+            pp.checkTopology(false, &points);
+
+            if (pTyp == TopoType::OPEN)
+            {
+                Info<< setw(34)
+                    << "ok (non-closed singly connected)";
+            }
+            else
+            {
+                Info<< setw(34)
+                    << "multiply connected (shared edge)";
+            }
+        }
+    }
+
+    if (allGeometry)
+    {
+        const labelList& mp = pp.meshPoints();
+
+        if (returnReduce(mp.size(), sumOp<label>()) > 0)
+        {
+            boundBox bb(pp.points(), mp, true); // reduce
+            Info<< ' ' << bb;
+        }
+    }
+}
+
+
 Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
@@ -466,6 +535,13 @@ Foam::label Foam::checkTopology
         }
     }
 
+    // Non-manifold points
+    pointSet points
+    (
+        mesh,
+        "nonManifoldPoints",
+        mesh.nPoints()/1000
+    );
 
     {
         if (!Pstream::parRun())
@@ -478,17 +554,8 @@ Foam::label Foam::checkTopology
             Info<< "\nChecking basic patch addressing..." << endl;
         }
 
-
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-        // Non-manifold points
-        pointSet points
-        (
-            mesh,
-            "nonManifoldPoints",
-            mesh.nPoints()/1000
-        );
-
         Pout.setf(ios_base::left);
 
         Info<< "    "
@@ -511,76 +578,127 @@ Foam::label Foam::checkTopology
 
             if (!isA<processorPolyPatch>(pp))
             {
-                Info<< "    "
+                checkPatch(allGeometry, pp.name(), pp, points);
-                    << setw(20) << pp.name()
-                    << setw(9) << returnReduce(pp.size(), sumOp<label>())
-                    << setw(9) << returnReduce(pp.nPoints(), sumOp<label>());
-
-                if (!Pstream::parRun())
-                {
-                    primitivePatch::surfaceTopo pTyp = pp.surfaceType();
-
-                    if (pp.empty())
-                    {
-                        Info<< setw(34) << "ok (empty)";
-                    }
-                    else if (pTyp == primitivePatch::MANIFOLD)
-                    {
-                        if (pp.checkPointManifold(true, &points))
-                        {
-                            Info<< setw(34)
-                                << "multiply connected (shared point)";
-                        }
-                        else
-                        {
-                            Info<< setw(34) << "ok (closed singly connected)";
-                        }
-
-                        // Add points on non-manifold edges to make set complete
-                        pp.checkTopology(false, &points);
-                    }
-                    else
-                    {
-                        pp.checkTopology(false, &points);
-
-                        if (pTyp == primitivePatch::OPEN)
-                        {
-                            Info<< setw(34)
-                                << "ok (non-closed singly connected)";
-                        }
-                        else
-                        {
-                            Info<< setw(34)
-                                << "multiply connected (shared edge)";
-                        }
-                    }
-                }
-
-                if (allGeometry)
-                {
-                    const labelList& mp = pp.meshPoints();
-
-                    if (returnReduce(mp.size(), sumOp<label>()) > 0)
-                    {
-                        boundBox bb(pp.points(), mp, true); // reduce
-                        Info<< ' ' << bb;
-                    }
-                }
                 Info<< endl;
             }
         }
 
-        if (points.size())
+        //Info.setf(ios_base::right);
-        {
-            Info<< "  <<Writing " << returnReduce(points.size(), sumOp<label>())
-                << " conflicting points to set "
-                << points.name() << endl;
-
-            points.instance() = mesh.pointsInstance();
-            points.write();
     }
 
-        //Info.setf(ios_base::right);
+    {
+        if (!Pstream::parRun())
+        {
+            Info<< "\nChecking faceZone topology for multiply connected"
+                << " surfaces..." << endl;
+        }
+        else
+        {
+            Info<< "\nChecking basic faceZone addressing..." << endl;
+        }
+
+        Pout.setf(ios_base::left);
+
+        const faceZoneMesh& faceZones = mesh.faceZones();
+
+        if (faceZones.size())
+        {
+            Info<< "    "
+                << setw(20) << "FaceZone"
+                << setw(9) << "Faces"
+                << setw(9) << "Points";
+
+            if (!Pstream::parRun())
+            {
+                Info<< setw(34) << "Surface topology";
+            }
+            if (allGeometry)
+            {
+                Info<< " Bounding box";
+            }
+            Info<< endl;
+
+            forAll(faceZones, zoneI)
+            {
+                const faceZone& fz = faceZones[zoneI];
+                checkPatch(allGeometry, fz.name(), fz(), points);
+                Info<< endl;
+            }
+        }
+        else
+        {
+            Info<< "    No faceZones found."<<endl;
+        }
+    }
+
+    label nPoints = returnReduce(points.size(), sumOp<label>());
+
+    if (nPoints)
+    {
+        Info<< "  <<Writing " << nPoints
+            << " conflicting points to set " << points.name() << endl;
+        points.instance() = mesh.pointsInstance();
+        points.write();
+        if (setWriter.valid())
+        {
+            mergeAndWrite(setWriter, points);
+        }
+    }
+
+    {
+        Info<< "\nChecking basic cellZone addressing..." << endl;
+
+        Pout.setf(ios_base::left);
+
+        const cellZoneMesh& cellZones = mesh.cellZones();
+
+        if (cellZones.size())
+        {
+
+            Info<< "    "
+                << setw(20) << "CellZone"
+                << setw(9) << "Cells"
+                << setw(9) << "Points"
+                << setw(13) << "BoundingBox" <<endl;
+
+            const cellList& cells = mesh.cells();
+            const faceList& faces = mesh.faces();
+            treeBoundBox bb(boundBox::invertedBox);
+            PackedBoolList isZonePoint(mesh.nPoints());
+
+            forAll(cellZones, zoneI)
+            {
+                const cellZone& cZone = cellZones[zoneI];
+
+                forAll(cZone, i)
+                {
+                    const label cellI = cZone[i];
+                    const cell& cFaces = cells[cellI];
+                    forAll(cFaces, cFacei)
+                    {
+                        const face& f = faces[cFaces[cFacei]];
+                        forAll(f, fp)
+                        {
+                            if (isZonePoint.set(f[fp]))
+                            {
+                                bb.add(mesh.points()[f[fp]]);
+                            }
+                         }
+                    }
+                }
+
+                Info<< "    "
+                    << setw(20) << cZone.name()
+                    << setw(9) << returnReduce(cZone.size(), sumOp<label>())
+                    << setw(9) 
+                    << returnReduce(isZonePoint.count(), sumOp<label>())
+                    << setw(3) << bb << endl;
+            }
+        }
+        else
+        {
+            Info<< "    No cellZones found."<<endl;
+        }
     }
 
     // Force creation of all addressing if requested.

src/OpenFOAM/db/error/error.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -182,14 +182,6 @@ void Foam::error::exit(const int errNo)
         abort();
     }
 
-    if (Pstream::parRun())
-    {
-        Perr<< endl << *this << endl
-            << "\nFOAM parallel run exiting\n" << endl;
-        Pstream::exit(errNo);
-    }
-    else
-    {
     if (throwExceptions_)
     {
         // Make a copy of the error to throw
@@ -200,13 +192,18 @@ void Foam::error::exit(const int errNo)
 
         throw errorException;
     }
+    else if (Pstream::parRun())
+    {
+        Perr<< endl << *this << endl
+            << "\nFOAM parallel run exiting\n" << endl;
+        Pstream::exit(errNo);
+    }
     else
     {
         Perr<< endl << *this << endl
             << "\nFOAM exiting\n" << endl;
         ::exit(1);
     }
-    }
 }
 
 
@@ -226,15 +223,6 @@ void Foam::error::abort()
         ::abort();
     }
 
-    if (Pstream::parRun())
-    {
-        Perr<< endl << *this << endl
-            << "\nFOAM parallel run aborting\n" << endl;
-        printStack(Perr);
-        Pstream::abort();
-    }
-    else
-    {
     if (throwExceptions_)
     {
         // Make a copy of the error to throw
@@ -245,6 +233,13 @@ void Foam::error::abort()
 
         throw errorException;
     }
+    else if (Pstream::parRun())
+    {
+        Perr<< endl << *this << endl
+            << "\nFOAM parallel run aborting\n" << endl;
+        printStack(Perr);
+        Pstream::abort();
+    }
     else
     {
         Perr<< endl << *this << endl
@@ -252,7 +247,6 @@ void Foam::error::abort()
         printStack(Perr);
         ::abort();
     }
-    }
 }
 
 

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.C

@@ -791,7 +791,9 @@ bool Foam::functionObjectList::read()
                 FatalError.throwExceptions(throwingError);
                 FatalIOError.throwExceptions(throwingIOerr);
 
-                if (foPtr.valid())
+                // If one processor only has thrown an exception (so exited the
+                // constructor) invalidate the whole functionObject
+                if (returnReduce(foPtr.valid(), andOp<bool>()))
                 {
                     objPtr = foPtr.ptr();
                 }

src/OpenFOAM/meshes/boundBox/boundBoxI.H

@@ -59,7 +59,16 @@ inline Foam::boundBox::boundBox(Istream& is)
 
 inline bool Foam::boundBox::empty() const
 {
-    return (min_ > max_);
+    // Note: cannot use min_ > max_ here since that tests for -all- components
+    // of min_ being larger than max_.
+    for (direction dir = 0; dir < vector::nComponents; ++dir)
+    {
+        if (min_[dir] > max_[dir])
+        {
+            return true;
+        }
+    }
+    return false;
 }
 
 

src/functionObjects/field/fluxSummary/fluxSummary.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -258,7 +258,7 @@ void Foam::functionObjects::fluxSummary::initialiseFaceZone
         if (faceID >= 0)
         {
             // Orientation set by faceZone flip map
-            if (fZone.flipMap()[facei])
+            if (fZone.flipMap()[i])
             {
                 flips.append(true);
             }

src/functionObjects/field/nearWallFields/nearWallFields.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,6 +77,7 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
 
         vectorField nf(patch.nf());
         vectorField faceCellCentres(patch.patch().faceCellCentres());
+        const labelUList& faceCells = patch.patch().faceCells();
 
         forAll(patch, patchFacei)
         {
@@ -95,25 +96,47 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
             );
 
 
+            // Starting point and tet
             point start;
+            label tetFacei = -1;
+            label tetPti = -1;
+            const label celli = faceCells[patchFacei];
+
             if (startInfo.hit())
             {
+                // Move start point slightly in so it is inside the tet
+                const face& f = mesh_.faces()[meshFacei];
+
+                tetFacei = meshFacei;
+                tetPti = (startInfo.index()+1) % f.size();
+
                 start = startInfo.hitPoint();
+                //// Uncomment below to shift slightly in:
+                //tetIndices tet(celli, meshFacei, tetPti, mesh_);
+                //start =
+                //    (1.0-1e-6)*startInfo.hitPoint()
+                //  + 1e-6*tet.tet(mesh_).centre();
             }
             else
             {
                 // Fallback: start tracking from neighbouring cell centre
                 start = faceCellCentres[patchFacei];
+                mesh_.findTetFacePt(celli, start, tetFacei, tetPti);
             }
 
             const point end = start-distance_*nf[patchFacei];
 
-            // Find tet for starting location
+            if (tetFacei == -1)
-            label celli = -1;
+            {
-            label tetFacei = -1;
+                WarningInFunction << "Did not find point " << start
-            label tetPti = -1;
+                    << " inside cell " << celli
-            mesh_.findCellFacePt(start, celli, tetFacei, tetPti);
+                    << ". Not seeding particle originating from face centre "
-
+                    << mesh_.faceCentres()[meshFacei]
+                    << " with cell centre " << mesh_.cellCentres()[celli]
+                    << endl;
+            }
+            else
+            {
                 // Add to cloud. Add originating face as passive data
                 cloud.addParticle
                 (
@@ -128,7 +151,7 @@ void Foam::functionObjects::nearWallFields::calcAddressing()
                         globalWalls.toGlobal(nPatchFaces)    // passive data
                     )
                 );
-
+            }
             nPatchFaces++;
         }
     }

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -683,7 +683,8 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
             IOobject::MUST_READ,
             IOobject::NO_WRITE,
             false
-        )
+        ),
+        dict
     ),
     sampleSource_(samplingSourceNames_.lookup("source", dict)),
     needsUpdate_(true),
@@ -779,15 +780,28 @@ bool Foam::sampledTriSurfaceMesh::update()
         surface_.triSurface::meshPoints()
     );
 
-    bb.intersect(mesh().bounds());
+    // Check for overlap with (global!) mesh bb
+    const bool intersect = bb.intersect(mesh().bounds());
 
-    // Extend a bit
+    if (!intersect)
+    {
+        // Surface and mesh do not overlap at all. Guarantee a valid
+        // bounding box so we don't get any 'invalid bounding box' errors.
+        bb = treeBoundBox(mesh().bounds());
         const vector span(bb.span());
 
+        bb.min() += (0.5-1e-6)*span;
+        bb.max() -= (0.5-1e-6)*span;
+    }
+    else
+    {
+        // Extend a bit
+        const vector span(bb.span());
         bb.min() -= 0.5*span;
         bb.max() += 0.5*span;
 
         bb.inflate(1e-6);
+    }
 
     // Mesh search engine, no triangulation of faces.
     meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);