ENH: finiteArea: improve robustness in code sections vulnerable to math errors

It has been observed that the finite-area framework is prone to numerical
issues when zero-valued edge lenghts, edge/face normals and face areas exist.

To improve exception handling at identified code sections to gracefully
overcome math errors, the problematic entities are lower-bounded by SMALL.

applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H

@@ -4,7 +4,7 @@
         sqrt
         (
             2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
-            *aMesh.edgeInterpolation::deltaCoeffs()
+            *mag(aMesh.edgeInterpolation::deltaCoeffs())
             /rhol
         )
     ).value()*runTime.deltaT().value();

src/finiteArea/faMesh/faMeshDemandDrivenData.C

@@ -594,6 +594,12 @@ Foam::tmp<Foam::vectorField> Foam::faMesh::calcRawEdgeNormals(int order) const
 
         edgeNormals[edgei].removeCollinear(edgeLine.unitVec());
         edgeNormals[edgei].normalise();
+
+        // Do not allow any mag(val) < SMALL
+        if (mag(edgeNormals[edgei]) < SMALL)
+        {
+            edgeNormals[edgei] = vector::uniform(SMALL);
+        }
     }
 
     return tedgeNormals;
@@ -658,6 +664,12 @@ void Foam::faMesh::calcLe() const
                 edges_[edgei].line(localPoints),
                 edgeNormals[edgei]
             );
+
+            // Do not allow any mag(val) < SMALL
+            if (mag(leVectors[edgei]) < SMALL)
+            {
+                leVectors[edgei] = vector::uniform(SMALL);
+            }
         }
 
         // Copy internal field
@@ -672,6 +684,15 @@ void Foam::faMesh::calcLe() const
         {
             const faPatch& fap = boundary()[patchi];
             bfld[patchi] = fap.patchRawSlice(leVectors);
+
+            for (auto& patchEdge : bfld[patchi])
+            {
+                // Do not allow any mag(val) < SMALL
+                if (mag(patchEdge) < SMALL)
+                {
+                    patchEdge = vector::uniform(SMALL);
+                }
+            }
         }
     }
     else
@@ -692,6 +713,12 @@ void Foam::faMesh::calcLe() const
                     edges_[edgei].line(localPoints),
                     edgeNormals[edgei]
                 );
+
+                // Do not allow any mag(val) < SMALL
+                if (mag(fld[edgei]) < SMALL)
+                {
+                    fld[edgei] = vector::uniform(SMALL);
+                }
             }
         }
 
@@ -715,6 +742,12 @@ void Foam::faMesh::calcLe() const
                     bndEdgeNormals[patchEdgei]
                 );
 
+                // Do not allow any mag(val) < SMALL
+                if (mag(pfld[patchEdgei]) < SMALL)
+                {
+                    pfld[patchEdgei] = vector::uniform(SMALL);
+                }
+
                 ++edgei;
             }
         }
@@ -759,6 +792,13 @@ void Foam::faMesh::calcMagLe() const
         for (const edge& e : internalEdges())
         {
             *iter = e.mag(localPoints);
+
+            // Do not allow any mag(val) < SMALL
+            if (mag(*iter) < SMALL)
+            {
+                *iter = SMALL;
+            }
+
             ++iter;
         }
     }
@@ -774,6 +814,13 @@ void Foam::faMesh::calcMagLe() const
             for (const edge& e : boundary()[patchi].patchSlice(edges_))
             {
                 *iter = e.mag(localPoints);
+
+                // Do not allow any mag(val) < SMALL
+                if (mag(*iter) < SMALL)
+                {
+                    *iter = SMALL;
+                }
+
                 ++iter;
             }
         }
@@ -955,6 +1002,12 @@ void Foam::faMesh::calcS() const
         forAll(fld, facei)
         {
             fld[facei] = Foam::mag(meshFaceAreas[facei]);
+
+            // Do not allow any mag(val) < SMALL
+            if (mag(fld[facei]) < SMALL)
+            {
+                fld[facei] = SMALL;
+            }
         }
     }
     else
@@ -968,6 +1021,13 @@ void Foam::faMesh::calcS() const
         for (const face& f : faces())
         {
             *iter = f.mag(localPoints);
+
+            // Do not allow any mag(val) < SMALL
+            if (mag(*iter) < SMALL)
+            {
+                *iter = SMALL;
+            }
+
             ++iter;
         }
     }
@@ -1028,6 +1088,15 @@ void Foam::faMesh::calcFaceAreaNormals() const
 
         // Make unit normals
         fld.normalise();
+
+        for (auto& f : fld)
+        {
+            // Do not allow any mag(val) < SMALL
+            if (mag(f) < SMALL)
+            {
+                f = vector::uniform(SMALL);
+            }
+        }
     }
 
 
@@ -1122,6 +1191,12 @@ void Foam::faMesh::calcEdgeAreaNormals() const
 
             fld[edgei].removeCollinear(edgeLine.unitVec());
             fld[edgei].normalise();
+
+            // Do not allow any mag(val) < SMALL
+            if (mag(fld[edgei]) < SMALL)
+            {
+                fld[edgei] = vector::uniform(SMALL);
+            }
         }
     }
 
@@ -1148,6 +1223,12 @@ void Foam::faMesh::calcEdgeAreaNormals() const
 
                 pfld[patchEdgei].removeCollinear(edgeLine.unitVec());
                 pfld[patchEdgei].normalise();
+
+                // Do not allow any mag(val) < SMALL
+                if (mag(pfld[patchEdgei]) < SMALL)
+                {
+                    pfld[patchEdgei] = vector::uniform(SMALL);
+                }
             }
         }
     }

src/finiteArea/faMesh/faPatches/faPatch/faPatch.C

@@ -478,7 +478,19 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
 {
     // Use patch-normal delta for all non-coupled BCs
     const vectorField nHat(edgeNormals());
-    return nHat*(nHat & (edgeCentres() - edgeFaceCentres()));
+
+    vectorField edgePN(edgeCentres() - edgeFaceCentres());
+
+    // Do not allow any mag(val) < SMALL
+    for (vector& edgei : edgePN)
+    {
+        if (mag(edgei) < SMALL)
+        {
+            edgei = vector::uniform(SMALL);
+        }
+    }
+
+    return nHat*(nHat & edgePN);
 }
 
 

src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C

@@ -115,6 +115,13 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
         label nei = neighbour[facei];
 
         vector d = C[nei] - C[own];
+
+        // Do not allow any mag(val) < SMALL
+        if (mag(d) < SMALL)
+        {
+            d = vector::uniform(SMALL);
+        }
+
         symmTensor wdd = (1.0/magSqr(d))*sqr(d);
 
         dd[own] += wdd;
@@ -169,6 +176,13 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
         label nei = neighbour[facei];
 
         vector d = C[nei] - C[own];
+
+        // Do not allow any mag(val) < SMALL
+        if (mag(d) < SMALL)
+        {
+            d = vector::uniform(SMALL);
+        }
+
         scalar magSfByMagSqrd = 1.0/magSqr(d);
 
         lsP[facei] = magSfByMagSqrd*(invDd[own] & d);

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

@@ -248,6 +248,12 @@ void Foam::edgeInterpolation::makeLPN() const
             );
 
         lPNIn[edgeI] = (lPE + lEN);
+
+        // Do not allow any mag(val) < SMALL
+        if (mag(lPNIn[edgeI]) < SMALL)
+        {
+            lPNIn[edgeI] = SMALL;
+        }
     }
 
 
@@ -287,7 +293,7 @@ void Foam::edgeInterpolation::makeWeights() const
             false
         ),
         mesh(),
-        dimless
+        dimensionedScalar(dimless, 1)
     );
     edgeScalarField& weightingFactors = *weightingFactors_;
 
@@ -323,15 +329,12 @@ void Foam::edgeInterpolation::makeWeights() const
               + skewCorrEdge
             );
 
-        weightingFactorsIn[edgeI] =
-            lEN
-            /(
-                lPE
-#               ifdef BAD_MESH_STABILISATION
-              + VSMALL
-#               endif
-              + lEN
-            );
+        // weight = (0,1]
+        const scalar lPN = lPE + lEN;
+        if (mag(lPN) > SMALL)
+        {
+            weightingFactorsIn[edgeI] = lEN/lPN;
+        }
     }
 
     forAll(mesh().boundary(), patchI)
@@ -370,7 +373,7 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
             false
         ),
         mesh(),
-        dimless/dimLength
+        dimensionedScalar(dimless/dimLength, SMALL)
     );
     edgeScalarField& DeltaCoeffs = *differenceFactors_;
     scalarField& dc = DeltaCoeffs.primitiveFieldRef();
@@ -427,11 +430,12 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
         // Edge normal - area tangent
         edgeNormal = normalised(lengths[edgeI]);
 
-        // Delta coefficient as per definition
-//         dc[edgeI] = 1.0/(lPN*(unitDelta & edgeNormal));
-
-        // Stabilised form for bad meshes.  HJ, 23/Jul/2009
-        dc[edgeI] = 1.0/max((lPN*(unitDelta & edgeNormal)), 0.05*lPN);
+        // Do not allow any mag(val) < SMALL
+        const scalar alpha = lPN*(unitDelta & edgeNormal);
+        if (mag(alpha) > SMALL)
+        {
+            dc[edgeI] = scalar(1)/max(alpha, 0.05*lPN);
+        }
     }
 
 
@@ -478,7 +482,7 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
     const edgeList& edges = mesh().edges();
     const pointField& points = mesh().points();
 
-    scalarField deltaCoeffs(owner.size());
+    scalarField deltaCoeffs(owner.size(), SMALL);
 
     vectorField& CorrVecsIn = CorrVecs.primitiveFieldRef();
 
@@ -499,8 +503,12 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
         // Edge normal - area tangent
         edgeNormal = normalised(lengths[edgeI]);
 
-        // Delta coeffs
-        deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal);
+        // Do not allow any mag(val) < SMALL
+        const scalar alpha = unitDelta & edgeNormal;
+        if (mag(alpha) > SMALL)
+        {
+            deltaCoeffs[edgeI] = scalar(1)/alpha;
+        }
 
         // Edge correction vector
         CorrVecsIn[edgeI] =
@@ -570,7 +578,7 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
             false
         ),
         mesh(),
-        dimless
+        dimensionedVector(dimless, Zero)
     );
     edgeVectorField& SkewCorrVecs = *skewCorrectionVectors_;
 
@@ -590,19 +598,22 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
         const vector& P = C[owner[edgeI]];
         const vector& N = C[neighbour[edgeI]];
         const vector& S = points[edges[edgeI].start()];
-        vector e = edges[edgeI].vec(points);
 
-        const scalar beta = magSqr((N - P)^e);
+        // (T:Eq. 5.4)
+        const vector d(N - P);
+        const vector e(edges[edgeI].vec(points));
+        const vector de(d^e);
+        const scalar alpha = magSqr(de);
 
-        if (beta < ROOTVSMALL)
+        if (alpha < SMALL)
         {
             // Too small - skew correction remains zero
             continue;
         }
+        const scalar beta = -((d^(S - P)) & de)/alpha;
 
-        const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
-
-        vector E(S + alpha*e);
+        // (T:Eq. 5.3)
+        const vector E(S + beta*e);
 
         SkewCorrVecs[edgeI] = Ce[edgeI] - E;
     }
@@ -630,28 +641,26 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
                 const vector& P = C[edgeFaces[edgeI]];
                 const vector& N = ngbC[edgeI];
                 const vector& S = points[patchEdges[edgeI].start()];
-                vector e = patchEdges[edgeI].vec(points);
 
-                const scalar beta = magSqr((N - P)^e);
+                // (T:Eq. 5.4)
+                const vector d(N - P);
+                const vector e(patchEdges[edgeI].vec(points));
+                const vector de(d^e);
+                const scalar alpha = magSqr(de);
 
-                if (beta < ROOTVSMALL)
+                if (alpha < SMALL)
                 {
                     // Too small - skew correction remains zero
                     continue;
                 }
+                const scalar beta = -((d^(S - P)) & de)/alpha;
 
-                const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
-
-                vector E(S + alpha*e);
+                const vector E(S + beta*e);
 
                 patchSkewCorrVecs[edgeI] =
                     Ce.boundaryField()[patchI][edgeI] - E;
             }
         }
-        else
-        {
-            patchSkewCorrVecs = vector::zero;
-        }
     }
 
     #ifdef FA_SKEW_CORRECTION

src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C

@@ -108,6 +108,12 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
         d.normalise();
         d *= mesh.edgeInterpolation::lPN().internalField()[edge];
 
+        // Do not allow any mag(val) < SMALL
+        if (mag(d) < SMALL)
+        {
+            d = vector::uniform(SMALL);
+        }
+
         weights[edge] =
             this->weight
             (
@@ -189,6 +195,12 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
                 d.normalise();
                 d *= pLPN[edgeI];
 
+                // Do not allow any mag(val) < SMALL
+                if (mag(d) < SMALL)
+                {
+                    d = vector::uniform(SMALL);
+                }
+
                 pWeights[edgeI] =
                     this->weight
                     (