ENH: processorFaPatch: use internal-edge algos for processor edges to ensure parallel consistency

- The edgeInterpolation::makeCorrectionVectors() disables the non-orthogonality
correction if the calculated non-orthogonality coefficient is below 0.1.
However, this activation routine only considers internal edges, and excludes
any processor edges, resulting in inconsistent parallel calculations. This
routine is removed.
- Fatal errors are replaced with zero-valued fields for non-orthogonality-
and skewness-correction routines.

src/finiteArea/faMesh/faPatches/basic/coupled/coupledFaPatch.H

@@ -79,9 +79,15 @@ protected:
         //- Make patch weighting factors
         virtual void makeWeights(scalarField&) const = 0;
 
+        //- Make patch geodesic distance between P and N
+        virtual void makeLPN(scalarField&) const = 0;
+
         //- Make patch face - neighbour cell distances
         virtual void makeDeltaCoeffs(scalarField&) const = 0;
 
+        //- Make non-orthogonality correction vectors
+        virtual void makeCorrectionVectors(vectorField&) const = 0;
+
         //- Calculate the uniform transformation tensors
         void calcTransformTensors
         (

src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C

@@ -218,6 +218,13 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
 }
 
 
+void Foam::cyclicFaPatch::makeLPN(scalarField& lPN) const
+{
+    makeDeltaCoeffs(lPN);
+    lPN = scalar(1)/lPN;
+}
+
+
 void Foam::cyclicFaPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     const scalarField deltas(edgeNormals() & coupledFaPatch::delta());

src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.H

@@ -78,9 +78,15 @@ protected:
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
+        //- Make patch geodesic distance between P and N
+        void makeLPN(scalarField&) const;
+
         //- Make patch face - neighbour cell distances
         void makeDeltaCoeffs(scalarField&) const;
 
+        //- Make non-orthogonality correction vectors
+        void makeCorrectionVectors(vectorField& cv) const { cv = Zero; }
+
 
 public:
 

src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C

@@ -398,18 +398,21 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const
 {
     if (Pstream::parRun())
     {
-        // The face normals point in the opposite direction on the other side
+        const vectorField& skew = skewCorrectionVectors();
-        scalarField neighbEdgeCentresCn
+        const scalarField& PN = lPN();
+
+        const scalarField lEN
         (
-            neighbEdgeNormals()
+            mag(neighbEdgeFaceCentres() - edgeCentres() + skew)
-          & (neighbEdgeCentres() - neighbEdgeFaceCentres())
         );
 
-        w = neighbEdgeCentresCn/
+        forAll(w, i)
-            (
+        {
-                (edgeNormals() & coupledFaPatch::delta())
+            if (mag(PN[i]) > SMALL)
-              + neighbEdgeCentresCn + VSMALL
+            {
-            );
+                w[i] = lEN[i]/PN[i];
+            }
+        }
     }
     else
     {
@@ -418,12 +421,61 @@ void Foam::processorFaPatch::makeWeights(scalarField& w) const
 }
 
 
+void Foam::processorFaPatch::makeLPN(scalarField& lPN) const
+{
+    const vectorField& skew = skewCorrectionVectors();
+
+    if (Pstream::parRun())
+    {
+        lPN =
+            mag(edgeCentres() - skew - edgeFaceCentres())         // lPE
+          + mag(neighbEdgeFaceCentres() - edgeCentres() + skew);  // lEN
+
+        for (scalar& lpn: lPN)
+        {
+            if (mag(lpn) < SMALL)
+            {
+                lpn = SMALL;
+            }
+        }
+    }
+    else
+    {
+        lPN =
+            mag(edgeCentres() - skew - edgeFaceCentres())         // lPE
+          + mag(edgeFaceCentres() - edgeCentres() + skew);        // lEN
+    }
+}
+
+
 void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     if (Pstream::parRun())
     {
-        dc = (1.0 - weights())
+        const edgeList& edges = boundaryMesh().mesh().edges();
-            /((edgeNormals() & coupledFaPatch::delta()) + VSMALL);
+        const pointField& points = boundaryMesh().mesh().points();
+        const vectorField& lengths = edgeLengths();
+        const scalarField& PN = lPN();
+
+        tmp<vectorField> tdelta = processorFaPatch::delta();
+        vectorField& unitDelta = tdelta.ref();
+
+        forAll(dc, i)
+        {
+            vector edgeNormal =
+                normalised(lengths[i] ^ edges[i + start()].vec(points));
+
+            unitDelta[i].removeCollinear(edgeNormal);
+            unitDelta[i].normalise();
+
+            edgeNormal = normalised(lengths[i]);
+
+            const scalar alpha = PN[i]*(edgeNormal & unitDelta[i]);
+            if (mag(alpha) > SMALL)
+            {
+                dc[i] = scalar(1)/max(alpha, 0.05*PN[i]);
+            }
+        }
     }
     else
     {
@@ -433,6 +485,44 @@ void Foam::processorFaPatch::makeDeltaCoeffs(scalarField& dc) const
 }
 
 
+void Foam::processorFaPatch::makeCorrectionVectors(vectorField& cv) const
+{
+    if (Pstream::parRun())
+    {
+        const edgeList& edges = boundaryMesh().mesh().edges();
+        const pointField& points = boundaryMesh().mesh().points();
+        const vectorField& lengths = edgeLengths();
+
+        tmp<vectorField> tdelta = processorFaPatch::delta();
+        vectorField& unitDelta = tdelta.ref();
+
+        forAll(cv, i)
+        {
+            vector edgeNormal =
+                normalised(lengths[i] ^ edges[i + start()].vec(points));
+
+            unitDelta[i].removeCollinear(edgeNormal);
+            unitDelta.normalise();
+
+            edgeNormal = normalised(lengths[i]);
+
+            const scalar alpha = unitDelta[i] & edgeNormal;
+            scalar dc = SMALL;
+            if (mag(alpha) > SMALL)
+            {
+                dc = scalar(1)/alpha;
+            }
+
+            cv[i] = edgeNormal - dc*unitDelta[i];
+        }
+    }
+    else
+    {
+        cv = Zero;
+    }
+}
+
+
 Foam::tmp<Foam::vectorField> Foam::processorFaPatch::delta() const
 {
     if (Pstream::parRun())

src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H

@@ -105,9 +105,15 @@ protected:
         //- Make patch weighting factors
         void makeWeights(scalarField&) const;
 
+        //- Make patch geodesic distance between P and N
+        void makeLPN(scalarField&) const;
+
         //- Make patch face - neighbour cell distances
         void makeDeltaCoeffs(scalarField&) const;
 
+        //- Make non-orthogonality correction vectors
+        void makeCorrectionVectors(vectorField&) const;
+
         //- Find non-globa patch points
         void makeNonGlobalPatchPoints() const;
 

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

@@ -500,12 +500,30 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
 }
 
 
+void Foam::faPatch::makeLPN(scalarField& lPN) const
+{
+    lPN = (edgeNormals() & delta());
+}
+
+
+const Foam::scalarField& Foam::faPatch::lPN() const
+{
+    return boundaryMesh().mesh().lPN().boundaryField()[index()];
+}
+
+
 void Foam::faPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     dc = scalar(1)/(edgeNormals() & delta());
 }
 
 
+const Foam::scalarField& Foam::faPatch::deltaCoeffs() const
+{
+    return boundaryMesh().mesh().deltaCoeffs().boundaryField()[index()];
+}
+
+
 void Foam::faPatch::makeCorrectionVectors(vectorField& k) const
 {
     const vectorField unitDelta(delta()/mag(delta()));
@@ -514,9 +532,10 @@ void Foam::faPatch::makeCorrectionVectors(vectorField& k) const
 }
 
 
-const Foam::scalarField& Foam::faPatch::deltaCoeffs() const
+const Foam::vectorField& Foam::faPatch::skewCorrectionVectors() const
 {
-    return boundaryMesh().mesh().deltaCoeffs().boundaryField()[index()];
+    return
+        boundaryMesh().mesh().skewCorrectionVectors().boundaryField()[index()];
 }
 
 

src/finiteArea/faMesh/faPatches/faPatch/faPatch.H

@@ -413,14 +413,24 @@ public:
             //- Return patch weighting factors
             const scalarField& weights() const;
 
+            //- Make patch geodesic distance between P and N
+            virtual void makeLPN(scalarField&) const;
+
+            //- Return patch geodesic distance between P and N
+            const scalarField& lPN() const;
+
             //- Make patch edge - neighbour face distances
             virtual void makeDeltaCoeffs(scalarField&) const;
 
-            void makeCorrectionVectors(vectorField&) const;
-
             //- Return patch edge - neighbour face distances
             const scalarField& deltaCoeffs() const;
 
+            //- Make non-orthogonality correction vectors
+            virtual void makeCorrectionVectors(vectorField&) const;
+
+            //- Return patch skew-correction vectors
+            const vectorField& skewCorrectionVectors() const;
+
 
         // Topological changes
 

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

@@ -98,9 +98,17 @@ const Foam::edgeVectorField& Foam::edgeInterpolation::correctionVectors() const
 {
     if (orthogonal())
     {
-        FatalErrorInFunction
+        return tmp<edgeVectorField>::New
-            << "cannot return correctionVectors; mesh is orthogonal"
+        (
-            << abort(FatalError);
+            IOobject
+            (
+                "correctionVectors",
+                mesh().pointsInstance(),
+                mesh().thisDb()
+            ),
+            mesh(),
+            dimensionedVector(dimless, Zero)
+        );
     }
 
     return (*correctionVectorsPtr_);
@@ -123,9 +131,17 @@ Foam::edgeInterpolation::skewCorrectionVectors() const
 {
     if (!skew())
     {
-        FatalErrorInFunction
+        return tmp<edgeVectorField>::New
-            << "cannot return skewCorrectionVectors; mesh is now skewed"
+        (
-            << abort(FatalError);
+            IOobject
+            (
+                "skewCorrectionVectors",
+                mesh().pointsInstance(),
+                mesh().thisDb()
+            ),
+            mesh(),
+            dimensionedVector(dimless, Zero)
+        );
     }
 
     return (*skewCorrectionVectorsPtr_);
@@ -233,12 +249,10 @@ void Foam::edgeInterpolation::makeLPN() const
 
     forAll(lPN.boundaryField(), patchI)
     {
-        mesh().boundary()[patchI].makeDeltaCoeffs
+        mesh().boundary()[patchI].makeLPN
         (
             lPN.boundaryFieldRef()[patchI]
         );
-
-        lPN.boundaryFieldRef()[patchI] = 1.0/lPN.boundaryField()[patchI];
     }
 
 
@@ -305,6 +319,7 @@ void Foam::edgeInterpolation::makeWeights() const
 
         // weight = (0,1]
         const scalar lPN = lPE + lEN;
+
         if (mag(lPN) > SMALL)
         {
             weightingFactorsIn[edgeI] = lEN/lPN;
@@ -498,31 +513,6 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
         mesh().boundary()[patchI].makeCorrectionVectors(CorrVecsbf[patchI]);
     }
 
-    scalar NonOrthogCoeff = 0.0;
-
-    if (owner.size() > 0)
-    {
-        scalarField sinAlpha(deltaCoeffs*mag(CorrVecs.internalField()));
-        sinAlpha.clamp_range(-1, 1);
-
-        NonOrthogCoeff = max(Foam::asin(sinAlpha)*180.0/M_PI);
-    }
-
-    reduce(NonOrthogCoeff, maxOp<scalar>());
-
-    DebugInFunction
-        << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
-        << endl;
-
-    if (NonOrthogCoeff < 0.1)
-    {
-        orthogonal_ = true;
-        correctionVectorsPtr_.reset(nullptr);
-    }
-    else
-    {
-        orthogonal_ = false;
-    }
 
     DebugInFunction
         << "Finished constructing non-orthogonal correction vectors"

tutorials/heatTransfer/buoyantPimpleFoam/hotRoomWithThermalShell/system/finite-area/faSchemes

@@ -21,7 +21,7 @@ ddtSchemes
 
 gradSchemes
 {
-    default         none;
+    default         Gauss linear;
     grad(jouleHeatingSource:V_ceilingShell) Gauss linear;
 }