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

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H

@@ -53,7 +53,6 @@ mesh.setFluxRequired(p.name());
         suppressDict.add("cellMask", true);
         suppressDict.add("cellDisplacement", true);
         suppressDict.add("interpolatedCells", true);
-        suppressDict.add("cellInterpolationWeight", true);
     }
 
     const_cast<dictionary&>

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H

@@ -239,11 +239,9 @@ forAll(isNeiInterpolatedFace, faceI)
                 n1 = vector(-n.z(), 0, n.x());
             }
         }
+        n1.normalise();
 
-        n1 /= mag(n1);
-
-        vector n2 = n ^ n1;
-        n2 /= mag(n2);
+        const vector n2 = normalised(n ^ n1);
 
         tensor rot =
             tensor

applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C

@@ -26,7 +26,7 @@ License
 #include "twoPhaseMixtureThermo.H"
 #include "gradientEnergyFvPatchScalarField.H"
 #include "mixedEnergyFvPatchScalarField.H"
-
+#include "collatedFileOperation.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -60,6 +60,10 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo
         T2.write();
     }
 
+    // Note: we're writing files to be read in immediately afterwards.
+    //       Avoid any thread-writing problems.
+    fileHandler().flush();
+
     thermo1_ = rhoThermo::New(U.mesh(), phase1Name());
     thermo2_ = rhoThermo::New(U.mesh(), phase2Name());
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C

@@ -165,8 +165,7 @@ void Foam::radiation::laserDTRM::initialise()
 
     const scalar t = mesh_.time().value();
     const vector lPosition = focalLaserPosition_->value(t);
-    vector lDir = laserDirection_->value(t);
-    lDir /= mag(lDir);
+    const vector lDir = normalised(laserDirection_->value(t));
 
     if (debug)
     {
@@ -175,7 +174,7 @@ void Foam::radiation::laserDTRM::initialise()
     }
 
     // Find a vector on the area plane. Normal to laser direction
-    vector rArea = vector::zero;
+    vector rArea = Zero;
     scalar magr = 0.0;
 
     {
@@ -188,7 +187,7 @@ void Foam::radiation::laserDTRM::initialise()
             magr = mag(rArea);
         }
     }
-    rArea /= mag(rArea);
+    rArea.normalise();
 
     scalar dr =  focalLaserRadius_/ndr_;
     scalar dTheta =  mathematical::twoPi/ndTheta_;

applications/test/cubicEqn/Test-cubicEqn.C

@@ -51,7 +51,7 @@ void test(const Type& polynomialEqn, const scalar tol)
             case roots::real:
                 v[i] = polynomialEqn.value(r[i]);
                 e[i] = polynomialEqn.error(r[i]);
-                ok = ok && mag(v[i]) < tol*mag(e[i]);
+                ok = ok && mag(v[i]) <= tol*mag(e[i]);
                 break;
             case roots::posInf:
                 v[i] = + inf;
@@ -79,8 +79,10 @@ void test(const Type& polynomialEqn, const scalar tol)
 
 int main()
 {
-    const int nTests = 1000000;
-    for (int t = 0; t < nTests; ++ t)
+    const scalar tol = 5;
+
+    const label nTests = 1000000;
+    for (label t = 0; t < nTests; ++ t)
     {
         test
         (
@@ -91,11 +93,26 @@ int main()
                 randomScalar(1e-50, 1e+50),
                 randomScalar(1e-50, 1e+50)
             ),
-            100
+            tol
         );
     }
+    Info << nTests << " random cubics tested" << endl;
 
-    Info << nTests << " cubics tested" << endl;
+    const label coeffMin = -9, coeffMax = 10, nCoeff = coeffMax - coeffMin;
+    for (label a = coeffMin; a < coeffMax; ++ a)
+    {
+        for (label b = coeffMin; b < coeffMax; ++ b)
+        {
+            for (label c = coeffMin; c < coeffMax; ++ c)
+            {
+                for (label d = coeffMin; d < coeffMax; ++ d)
+                {
+                    test(cubicEqn(a, b, c, d), tol);
+                }
+            }
+        }
+    }
+    Info<< nCoeff*nCoeff*nCoeff*nCoeff << " integer cubics tested" << endl;
 
     return 0;
 }

applications/test/faces/Test-faces.C

@@ -31,19 +31,81 @@ Description
 
 #include "argList.H"
 #include "labelledTri.H"
+#include "pointList.H"
+#include "ListOps.H"
 
 using namespace Foam;
 
+
+template<class Face>
+void faceInfo(const Face& f, const UList<point>& points)
+{
+    Info<< f
+        << " points:" << f.points(points)
+        << " normal:" << f.unitNormal(points);
+}
+
+
+template<class Face>
+void testSign
+(
+    const Face& f,
+    const UList<point>& points,
+    const UList<point>& testPoints
+)
+{
+    for (const point& p : testPoints)
+    {
+        Info<< "  point:" << p << " sign=" << f.sign(p, points) << nl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 //  Main program:
 
 int main(int argc, char *argv[])
 {
+    pointList points1
+    ({
+        { 0, 0, 0},
+        { -1, -1, 1},
+        { 1, -1, -1},
+        { 1, 1, -1},
+        { -1, 1, 1}
+    });
+
+    pointList points2 = ListOps::create<point>
+    (
+        points1,
+        [](const point& p){ return point(p.x(), p.y(), -p.z()); }
+    );
+
+    pointList testPoints
+    ({
+        { -2, -2, -2},
+        { -2, -2, 2},
+        { 0, 0, 0},
+        { 2, 2, -2},
+        { 2, 2, 2}
+    });
+
     face f1{ 1, 2, 3, 4 };
-    Info<< "face:" << f1 << nl;
+    Info<< "face:"; faceInfo(f1, points1); Info << nl;
+    testSign(f1, points1, testPoints);
+
+    Info<< "face:"; faceInfo(f1, points2); Info << nl;
+    testSign(f1, points2, testPoints);
+    Info<< nl;
 
     triFace t1{ 1, 2, 3 };
-    Info<< "triFace:" << t1 << nl;
+    Info<< "triFace:"; faceInfo(t1, points1); Info << nl;
+    testSign(t1, points1, testPoints);
+
+    Info<< "triFace:"; faceInfo(t1, points2); Info << nl;
+    testSign(t1, points2, testPoints);
+    Info<< nl;
+
 
     f1 = t1;
     Info<< "face:" << f1 << nl;

applications/utilities/mesh/advanced/splitCells/splitCells.C

@@ -135,8 +135,11 @@ bool largerAngle
     // Get cos between faceCentre and normal vector to determine in
     // which quadrant angle is. (Is correct for unwarped faces only!)
     // Correct for non-outwards pointing normal.
-    vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
-    c1c0 /= mag(c1c0) + VSMALL;
+    const vector c1c0 =
+        normalised
+        (
+            mesh.faceCentres()[f1] - mesh.faceCentres()[f0]
+        );
 
     scalar fcCosAngle = n0 & c1c0;
 
@@ -319,8 +322,7 @@ bool splitCell
 {
     const edge& e = mesh.edges()[edgeI];
 
-    vector eVec = e.vec(mesh.points());
-    eVec /= mag(eVec);
+    const vector eVec = e.unitVec(mesh.points());
 
     vector planeN = eVec ^ halfNorm;
 
@@ -328,8 +330,7 @@ bool splitCell
     // halfway on fully regular meshes (since we want cuts
     // to be snapped to vertices)
     planeN += 0.01*halfNorm;
-
-    planeN /= mag(planeN);
+    planeN.normalise();
 
     // Define plane through edge
     plane cutPlane(mesh.points()[e.start()], planeN);
@@ -410,11 +411,8 @@ void collectCuts
             label f0, f1;
             meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1);
 
-            vector n0 = faceAreas[f0];
-            n0 /= mag(n0);
-
-            vector n1 = faceAreas[f1];
-            n1 /= mag(n1);
+            const vector n0 = normalised(faceAreas[f0]);
+            const vector n1 = normalised(faceAreas[f1]);
 
             if
             (

applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C

@@ -188,19 +188,17 @@ label findInternalFace(const primitiveMesh& mesh, const labelList& meshF)
 bool correctOrientation(const pointField& points, const cellShape& shape)
 {
     // Get centre of shape.
-    point cc(shape.centre(points));
+    const point cc(shape.centre(points));
 
     // Get outwards pointing faces.
     faceList faces(shape.faces());
 
-    forAll(faces, i)
+    for (const face& f : faces)
     {
-        const face& f = faces[i];
-
-        vector n(f.normal(points));
+        const vector areaNorm(f.areaNormal(points));
 
         // Check if vector from any point on face to cc points outwards
-        if (((points[f[0]] - cc) & n) < 0)
+        if (((points[f[0]] - cc) & areaNorm) < 0)
         {
             // Incorrectly oriented
             return false;

applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C

@@ -225,9 +225,10 @@ int main(int argc, char *argv[])
                 // Determine orientation of tri v.s. cell centre.
                 point cc(cll.centre(points));
                 point fc(tri.centre(points));
-                vector fn(tri.normal(points));
 
-                if (((fc - cc) & fn) < 0)
+                const vector areaNorm(tri.areaNormal(points));
+
+                if (((fc - cc) & areaNorm) < 0)
                 {
                     // Boundary face points inwards. Flip.
                     boundaryFaces[facei].flip();

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C

@@ -45,7 +45,7 @@ Foam::cellAspectRatioControl::cellAspectRatioControl
     )
 {
     // Normalise the direction
-    aspectRatioDirection_ /= mag(aspectRatioDirection_) + SMALL;
+    aspectRatioDirection_.normalise();
 
     Info<< nl
         << "Cell Aspect Ratio Control" << nl

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C

@@ -428,7 +428,7 @@ void Foam::searchableSurfaceControl::initialVertices
 
                 if (mag(normals[0]) < SMALL)
                 {
-                    normals[0] = vector(1, 1, 1);
+                    normals[0] = vector::one;
                 }
 
                 pointAlignment.reset(new triad(normals[0]));

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -1137,7 +1137,7 @@ void Foam::conformalVoronoiMesh::move()
 
                         alignmentDirs[aA] = a + sign(dotProduct)*b;
 
-                        alignmentDirs[aA] /= mag(alignmentDirs[aA]);
+                        alignmentDirs[aA].normalise();
                     }
                 }
             }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -1832,11 +1832,11 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                         }
 
                         vector correctNormal = calcSharedPatchNormal(vc1, vc2);
-                        correctNormal /= mag(correctNormal);
+                        correctNormal.normalise();
 
                         Info<< "    cN " << correctNormal << endl;
 
-                        vector fN = f.normal(pts);
+                        vector fN = f.areaNormal(pts);
 
                         if (mag(fN) < SMALL)
                         {
@@ -1844,7 +1844,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                             continue;
                         }
 
-                        fN /= mag(fN);
+                        fN.normalise();
                         Info<< "    fN " << fN << endl;
 
                         if ((fN & correctNormal) > 0)

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C

@@ -220,8 +220,7 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
 //            << endl;
 
         // Calculate master point
-        vector masterPtVec(normalDir + nextNormalDir);
-        masterPtVec /= mag(masterPtVec) + SMALL;
+        const vector masterPtVec = normalised(normalDir + nextNormalDir);
 
         if
         (

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C

@@ -481,10 +481,9 @@ void Foam::conformalVoronoiMesh::calcFaceZones
                     norm
                 );
 
-                vector fN = faces[facei].normal(mesh.points());
-                fN /= mag(fN) + SMALL;
+                const vector areaNorm = faces[facei].areaNormal(mesh.points());
 
-                if ((norm[0] & fN) < 0)
+                if ((norm[0] & areaNorm) < 0)
                 {
                     flipMap[facei] = true;
                 }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -89,11 +89,9 @@ void Foam::conformationSurfaces::hasBoundedVolume
             Info<< "        Index = " << surfaces_[s] << endl;
             Info<< "        Offset = " << regionOffset_[s] << endl;
 
-            forAll(triSurf, sI)
+            for (const labelledTri& f : triSurf)
             {
-                const label patchID =
-                    triSurf[sI].region()
-                  + regionOffset_[s];
+                const label patchID = f.region() + regionOffset_[s];
 
                 // Don't include baffle surfaces in the calculation
                 if
@@ -102,15 +100,15 @@ void Foam::conformationSurfaces::hasBoundedVolume
                  != extendedFeatureEdgeMesh::BOTH
                 )
                 {
-                    sum += triSurf[sI].normal(surfPts);
+                    sum += f.areaNormal(surfPts);
                 }
                 else
                 {
-                    nBaffles++;
+                    ++nBaffles;
                 }
             }
             Info<< "        has " << nBaffles << " baffles out of "
-                << triSurf.size() << " triangles" << endl;
+                << triSurf.size() << " triangles" << nl;
 
             totalTriangles += triSurf.size();
         }
@@ -740,8 +738,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
 
                 surface.findNearest(sample, nearestDistSqr, info);
 
-                vector hitDir = info[0].rawPoint() - samplePts[i];
-                hitDir /= mag(hitDir) + SMALL;
+                const vector hitDir =
+                    normalised
+                    (
+                        info[0].rawPoint() - samplePts[i]
+                    );
 
                 pointIndexHit surfHit;
                 label hitSurface;

applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C

@@ -558,7 +558,7 @@ void Foam::CV2D::newPoints()
 
                     alignmentDirs[aA] = a + sign(dotProduct)*b;
 
-                    alignmentDirs[aA] /= mag(alignmentDirs[aA]);
+                    alignmentDirs[aA].normalise();
                 }
             }
         }
@@ -845,7 +845,7 @@ void Foam::CV2D::newPoints()
             cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x());
 
             // Normalise the primary coordinate direction
-            cd0 /= mag(cd0);
+            cd0.normalise();
 
             // Calculate the orthogonal coordinate direction
             vector2D cd1(-cd0.y(), cd0.x());

applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/insertFeaturePoints.C

@@ -123,7 +123,7 @@ void Foam::CV2D::insertFeaturePoints()
                 vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI));
 
                 vector2D cornerNormal = sum(fpn);
-                cornerNormal /= mag(cornerNormal);
+                cornerNormal.normalise();
 
                 if (debug)
                 {

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -665,6 +665,13 @@ addLayersControls
         // Default is 0.5*featureAngle. Set to -180 always attempt extrusion
         //layerTerminationAngle 25;
 
+        // Optional: disable shrinking of edges that have one (or two) points
+        // on an extruded patch.
+        // Default is false to enable single/two cell thick channels to still
+        // have layers. In <=1806 this was true by default. On larger gaps it
+        // should have no effect.
+        //disableWallEdges true;
+
         // Optional: at non-patched sides allow mesh to slip if extrusion
         // direction makes angle larger than slipFeatureAngle. Default is
         // 0.5*featureAngle.

applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C

@@ -112,7 +112,7 @@ int main(int argc, char *argv[])
     const word addRegion =
         args.lookupOrDefault<word>
         (
-            "masterRegion",
+            "addRegion",
             polyMesh::defaultRegion
         );