Merge commit 'OpenCFD/master' into olesenm

applications/solvers/combustion/dieselEngineFoam/createSpray.H

@@ -32,7 +32,5 @@ if (dieselSpray.twoD())
     gasMass0 *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
 }
 
-reduce(gasMass0, sumOp<scalar>());
-
 gasMass0 -=
     dieselSpray.injectedMass(runTime.value()) - dieselSpray.liquidMass();

applications/solvers/combustion/dieselEngineFoam/spraySummary.H

@@ -23,8 +23,6 @@
             gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
         }
 
-        reduce(gasMass, sumOp<scalar>());
-
         scalar addedMass = gasMass - gasMass0;
 
         Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -274,6 +274,14 @@ addLayersControls
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations
+    nLayerIter 50;
+
+    // Max number of iterations after which relaxed meshQuality controls
+    // get used.
+    nRelaxedIter 20;
 }
 
 
@@ -335,6 +343,16 @@ meshQualityControls
     nSmoothScale 4;
     //- amount to scale back displacement at error points
     errorReduction 0.75;
+
+
+
+    // Optional : some meshing phases allow usage of relaxed rules.
+    // See e.g. addLayersControls::nRelaxedIter.
+    relaxed
+    {
+        //- Maximum non-orthogonality allowed. Set to 180 to disable.
+        maxNonOrtho 75;
+    }
 }
 
 

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -171,11 +171,11 @@ int main(int argc, char *argv[])
 
     argList::validArgs.clear();
     argList::validArgs.append("surface file");
-    argList::validOptions.insert("noSelfIntersection", "");
+    argList::validOptions.insert("checkSelfIntersection", "");
     argList::validOptions.insert("verbose", "");
     argList args(argc, argv);
 
-    bool checkSelfIntersection = !args.options().found("noSelfIntersection");
+    bool checkSelfIntersection = args.options().found("checkSelfIntersection");
     bool verbose = args.options().found("verbose");
 
     fileName surfFileName(args.additionalArgs()[0]);

src/OpenFOAM/db/IOstreams/Sstreams/ISstream.C

@@ -224,7 +224,16 @@ Foam::Istream& Foam::ISstream::read(token& t)
                 }
                 else
                 {
-                    t = label(atol(numberBuffer));
+                    long lt = atol(numberBuffer);
+                    t = label(lt);
+
+                    // If the integer is too large to be represented as a label
+                    // return it as a scalar
+                    if (t.labelToken() != lt)
+                    {
+                        isScalar = true;
+                        t = scalar(atof(numberBuffer));
+                    }
                 }
             }
             else

src/OpenFOAM/meshes/boundBox/boundBox.H

@@ -174,7 +174,7 @@ public:
 
         // Query
 
-            //- Completely contains other boundingBox? (inside or on edge)
+            //- Overlaps/touches boundingBox?
             bool overlaps(const boundBox& bb) const
             {
                 return

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -222,16 +222,16 @@ public:
         ) const;
 
         //- Fast intersection with a ray.
-        //  For a hit, the pointHit.distance() is the line parameter t :
-        //  intersection=p+t*q. Only defined for FULL_RAY or
-        //  HALF_RAY.
+        //  Does face-center decomposition and returns triangle intersection
+        //  point closest to p. See triangle::intersection for details.
         pointHit intersection
         (
             const point& p,
             const vector& q,
             const point& ctr,
             const pointField& meshPoints,
-            const intersection::algorithm alg
+            const intersection::algorithm alg,
+            const scalar tol = 0.0
         ) const;
 
         //- Return nearest point to face

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -136,7 +136,8 @@ Foam::pointHit Foam::face::intersection
     const vector& q,
     const point& ctr,
     const pointField& meshPoints,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     scalar nearestHitDist = VGREAT;
@@ -154,7 +155,7 @@ Foam::pointHit Foam::face::intersection
             meshPoints[f[pI]],
             meshPoints[f[fcIndex(pI)]],
             ctr
-        ).intersection(p, q, alg);
+        ).intersection(p, q, alg, tol);
 
         if (curHit.hit())
         {

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -52,24 +52,13 @@ void Foam::mapDistribute::distribute
 
             if (domain != Pstream::myProcNo() && map.size())
             {
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, map);
             }
         }
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -126,13 +115,7 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize);
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -152,16 +135,8 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
-                const labelList& map = subMap[recvProc];
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, subMap[recvProc]);
             }
             else
             {
@@ -374,24 +349,13 @@ void Foam::mapDistribute::distribute
 
             if (domain != Pstream::myProcNo() && map.size())
             {
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, map);
             }
         }
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -449,13 +413,7 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize, nullValue);
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -475,16 +433,8 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
-                const labelList& map = subMap[recvProc];
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, subMap[recvProc]);
             }
             else
             {

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -153,7 +153,7 @@ public:
             inline scalar sweptVol(const triangle& t) const;
 
             //- Return point intersection with a ray.
-            //  For a hit, the distance is signed.  Positive number
+            //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.
             //  In case of miss pointHit is set to nearest point
             //  on triangle and its distance to the distance between
@@ -169,12 +169,14 @@ public:
             //- Fast intersection with a ray.
             //  For a hit, the pointHit.distance() is the line parameter t :
             //  intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or
-            //  HALF_RAY.
+            //  HALF_RAY. tol increases the virtual size of the triangle
+            // by a relative factor.
             inline pointHit intersection
             (
                 const point& p,
                 const vector& q,
-                const intersection::algorithm alg
+                const intersection::algorithm alg,
+                const scalar tol = 0.0
             ) const;
 
             //- Return nearest point to p on triangle

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -439,7 +439,8 @@ inline pointHit triangle<Point, PointRef>::intersection
 (
     const point& orig,
     const vector& dir,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     const vector edge1 = b_ - a_;
@@ -457,99 +458,61 @@ inline pointHit triangle<Point, PointRef>::intersection
     if (alg == intersection::VISIBLE)
     {
         // Culling branch
-        if (det < SMALL)
+        if (det < ROOTVSMALL)
         {
-            // return miss
+            // Ray on wrong side of triangle. Return miss
             return intersection;
         }
-        /* calculate distance from a_ to ray origin */
-        const vector tVec = orig-a_;
-
-        /* calculate U parameter and test bounds */
-        scalar u = tVec & pVec;
-
-        if (u < 0.0 || u > det)
-        {
-            // return miss
-            return intersection;
-        }
-
-        /* prepare to test V parameter */
-        const vector qVec = tVec ^ edge1;
-
-        /* calculate V parameter and test bounds */
-        scalar v = dir & qVec;
-
-        if (v < 0.0 || u + v > det)
-        {
-            // return miss
-            return intersection;
-        }
-
-        /* calculate t, scale parameters, ray intersects triangle */
-        scalar t = edge2 & qVec;
-        scalar inv_det = 1.0 / det;
-        t *= inv_det;
-        u *= inv_det;
-        v *= inv_det;
-
-        intersection.setHit();
-        intersection.setPoint(a_ + u*edge1 + v*edge2);
-        intersection.setDistance(t);
     }
     else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
     {
         // Non-culling branch
-        if (det > -SMALL && det < SMALL)
+        if (det > -ROOTVSMALL && det < ROOTVSMALL)
         {
-            // return miss
+            // Ray parallel to triangle. Return miss
             return intersection;
         }
-        const scalar inv_det = 1.0 / det;
+    }
 
-        /* calculate distance from a_ to ray origin */
-        const vector tVec = orig - a_;
-        /* calculate U parameter and test bounds */
-        const scalar u = (tVec & pVec)*inv_det;
+    const scalar inv_det = 1.0 / det;
 
-        if (u < 0.0 || u > 1.0)
-        {
-            // return miss
-            return intersection;
-        }
-        /* prepare to test V parameter */
-        const vector qVec = tVec ^ edge1;
-        /* calculate V parameter and test bounds */
-        const scalar v = (dir & qVec) * inv_det;
+    /* calculate distance from a_ to ray origin */
+    const vector tVec = orig-a_;
 
-        if (v < 0.0 || u + v > 1.0)
-        {
-            // return miss
-            return intersection;
-        }
-        /* calculate t, ray intersects triangle */
-        const scalar t = (edge2 & qVec) * inv_det;
+    /* calculate U parameter and test bounds */
+    const scalar u = (tVec & pVec)*inv_det;
 
-        if (alg == intersection::HALF_RAY && t < 0)
-        {
-            // return miss
-            return intersection;
-        }
+    if (u < -tol || u > 1.0+tol)
+    {
+        // return miss
+        return intersection;
+    }
 
-        intersection.setHit();
-        intersection.setPoint(a_ + u*edge1 + v*edge2);
-        intersection.setDistance(t);
+    /* prepare to test V parameter */
+    const vector qVec = tVec ^ edge1;
+
+    /* calculate V parameter and test bounds */
+    const scalar v = (dir & qVec) * inv_det;
+
+    if (v < -tol || u + v > 1.0+tol)
+    {
+        // return miss
+        return intersection;
     }
-    else
+
+    /* calculate t, scale parameters, ray intersects triangle */
+    const scalar t = (edge2 & qVec) * inv_det;
+
+    if (alg == intersection::HALF_RAY && t < -tol)
     {
-        FatalErrorIn
-        (
-            "triangle<Point, PointRef>::intersection(const point&"
-            ", const vector&, const intersection::algorithm)"
-        )   << "intersection only defined for VISIBLE, FULL_RAY or HALF_RAY"
-            << abort(FatalError);
+        // Wrong side of orig. Return miss
+        return intersection;
     }
 
+    intersection.setHit();
+    intersection.setPoint(a_ + u*edge1 + v*edge2);
+    intersection.setDistance(t);
+
     return intersection;
 }
 
@@ -622,7 +585,7 @@ inline bool triangle<Point, PointRef>::classify
 
     bool hit = false;
 
-    if (Foam::mag(u1) < SMALL)
+    if (Foam::mag(u1) < ROOTVSMALL)
     {
         beta = u0/u2;
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H

@@ -38,7 +38,6 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "dictionary.H"
-#include "boolList.H"
 #include "wallPoint.H"
 #include "searchableSurfaces.H"
 #include "refinementSurfaces.H"

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -45,6 +45,7 @@ Description
 #include "OFstream.H"
 #include "layerParameters.H"
 #include "combineFaces.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1414,8 +1415,11 @@ void Foam::autoLayerDriver::calculateLayerThickness
     const indirectPrimitivePatch& pp,
     const labelList& patchIDs,
     const scalarField& patchExpansionRatio,
-    const scalarField& patchFinalLayerRatio,
-    const scalarField& patchRelMinThickness,
+
+    const bool relativeSizes,
+    const scalarField& patchFinalLayerThickness,
+    const scalarField& patchMinThickness,
+
     const labelList& cellLevel,
     const labelList& patchNLayers,
     const scalar edge0Len,
@@ -1428,22 +1432,100 @@ void Foam::autoLayerDriver::calculateLayerThickness
     const fvMesh& mesh = meshRefiner_.mesh();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
+
+    // Rework patch-wise layer parameters into minimum per point
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Reuse input fields
+    expansionRatio.setSize(pp.nPoints());
+    expansionRatio = GREAT;
+    thickness.setSize(pp.nPoints());
+    thickness = GREAT;
+    minThickness.setSize(pp.nPoints());
+    minThickness = GREAT;
+
+    forAll(patchIDs, i)
     {
-        FatalErrorIn("calculateLayerThickness(..)")
-            << "Thickness should be factor of local undistorted cell size."
-            << " Valid values are [0..2]." << nl
-            << " minThickness:" << patchRelMinThickness
-            << exit(FatalError);
+        label patchI = patchIDs[i];
+
+        const labelList& meshPoints = patches[patchI].meshPoints();
+
+        forAll(meshPoints, patchPointI)
+        {
+            label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+
+            expansionRatio[ppPointI] = min
+            (
+                expansionRatio[ppPointI],
+                patchExpansionRatio[patchI]
+            );
+            thickness[ppPointI] = min
+            (
+                thickness[ppPointI],
+                patchFinalLayerThickness[patchI]
+            );
+            minThickness[ppPointI] = min
+            (
+                minThickness[ppPointI],
+                patchMinThickness[patchI]
+            );
+        }
     }
 
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        expansionRatio,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        thickness,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        minThickness,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+
 
-    // Per point the max cell level of connected cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Now the thicknesses are set according to the minimum of connected
+    // patches.
 
-    labelList maxPointLevel(pp.nPoints(), labelMin);
 
+    // Rework relative thickness into absolute
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // by multiplying with the internal cell size.
+
+    if (relativeSizes)
     {
+        if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
+        {
+            FatalErrorIn("calculateLayerThickness(..)")
+                << "Thickness should be factor of local undistorted cell size."
+                << " Valid values are [0..2]." << nl
+                << " minThickness:" << patchMinThickness
+                << exit(FatalError);
+        }
+
+
+        // Determine per point the max cell level of connected cells
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        labelList maxPointLevel(pp.nPoints(), labelMin);
+
         forAll(pp, i)
         {
             label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
@@ -1465,113 +1547,44 @@ void Foam::autoLayerDriver::calculateLayerThickness
             labelMin,           // null value
             false               // no separation
         );
-    }
 
 
-    // Rework patch-wise layer parameters into minimum per point
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    expansionRatio.setSize(pp.nPoints());
-    expansionRatio = GREAT;
-    scalarField finalLayerRatio(pp.nPoints(), GREAT);
-    scalarField relMinThickness(pp.nPoints(), GREAT);
-
-    {
-        forAll(patchIDs, i)