Corrected delta coefficients and weighting factors on coupled patches. Henrik Rusche

src/finiteVolume/fvMesh/fvPatches/constraint/cyclic/cyclicFvPatch.C

@@ -46,7 +46,12 @@ void cyclicFvPatch::makeWeights(scalarField& w) const
 {
     const scalarField& magFa = magSf();
 
-    scalarField deltas = nf() & fvPatch::delta();
+    // Note: mag in the dot-product.
+    // For all valid meshes, the non-orthogonality will be less that
+    // 90 deg and the dot-product will be positive.  For invalid
+    // meshes (d & s <= 0), this will stabilise the calculation
+    // but the result will be poor.  HJ, 24/Aug/2011
+    scalarField deltas = mag(nf() & fvPatch::delta());
     label sizeby2 = deltas.size()/2;
 
     scalar maxMatchError = 0;
@@ -100,15 +105,14 @@ void cyclicFvPatch::makeWeights(scalarField& w) const
 // Make patch face - neighbour cell distances
 void cyclicFvPatch::makeDeltaCoeffs(scalarField& dc) const
 {
-    scalarField deltas = nf() & fvPatch::delta();
+    vectorField d = delta();
-    label sizeby2 = deltas.size()/2;
+    vectorField n = nf();
+    label sizeby2 = d.size()/2;
 
     for (label facei = 0; facei < sizeby2; facei++)
     {
-        scalar di = deltas[facei];
+        // Stabilised form for bad meshes.  HJ, 24/Aug/2011
-        scalar dni = deltas[facei + sizeby2];
+        dc[facei] = 1.0/mag(n[facei] & d[facei]), 0.05*mag(d[facei]);
-
-        dc[facei] = 1.0/(di + dni);
         dc[facei + sizeby2] = dc[facei];
     }
 }

src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.C

@@ -54,10 +54,19 @@ void Foam::cyclicGgiFvPatch::makeWeights(scalarField& w) const
     if (cyclicGgiPolyPatch_.master())
     {
         vectorField n = nf();
-        scalarField nfc =
-            n & (cyclicGgiPolyPatch_.reconFaceCellCentres() - Cf());
 
-        w = nfc/((n & (Cf() - Cn())) + nfc);
+        // Note: mag in the dot-product.
+        // For all valid meshes, the non-orthogonality will be less that
+        // 90 deg and the dot-product will be positive.  For invalid
+        // meshes (d & s <= 0), this will stabilise the calculation
+        // but the result will be poor.  HJ, 24/Aug/2011
+        scalarField nfc =
+            mag
+            (
+                n & (cyclicGgiPolyPatch_.reconFaceCellCentres() - Cf())
+            );
+
+        w = nfc/(mag(n & (Cf() - Cn())) + nfc);
 
         if (bridgeOverlap())
         {
@@ -89,7 +98,10 @@ void Foam::cyclicGgiFvPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     if (cyclicGgiPolyPatch_.master())
     {
-        dc = 1.0/max(nf() & delta(), 0.05*mag(delta()));
+        // Stabilised form for bad meshes.  HJ, 24/Aug/2011
+        vectorField d = delta();
+
+        dc = 1.0/max(nf() & d, 0.05*mag(d));
 
         if (bridgeOverlap())
         {

src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C

@@ -64,9 +64,16 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
     if (ggiPolyPatch_.master())
     {
         vectorField n = nf();
-        scalarField nfc = n & (ggiPolyPatch_.reconFaceCellCentres() - Cf());
 
-        w = nfc/((n & (Cf() - Cn())) + nfc);
+        // Note: mag in the dot-product.
+        // For all valid meshes, the non-orthogonality will be less that
+        // 90 deg and the dot-product will be positive.  For invalid
+        // meshes (d & s <= 0), this will stabilise the calculation
+        // but the result will be poor.  HJ, 24/Aug/2011
+        scalarField nfc =
+            mag(n & (ggiPolyPatch_.reconFaceCellCentres() - Cf()));
+
+        w = nfc/(mag(n & (Cf() - Cn())) + nfc);
 
         if (bridgeOverlap())
         {
@@ -100,7 +107,10 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     if (ggiPolyPatch_.master())
     {
-        dc = 1.0/max(nf() & delta(), 0.05*mag(delta()));
+        // Stabilised form for bad meshes.  HJ, 24/Aug/2011
+        vectorField d = delta();
+
+        dc = 1.0/max(nf() & d, 0.05*mag(d));
 
         if (bridgeOverlap())
         {

src/finiteVolume/fvMesh/fvPatches/constraint/overlapGgi/overlapGgiFvPatch.C

@@ -51,9 +51,16 @@ void Foam::overlapGgiFvPatch::makeWeights(scalarField& w) const
     if (overlapGgiPolyPatch_.master())
     {
         vectorField n = nf();
+
+        // Note: mag in the dot-product.
+        // For all valid meshes, the non-orthogonality will be less that
+        // 90 deg and the dot-product will be positive.  For invalid
+        // meshes (d & s <= 0), this will stabilise the calculation
+        // but the result will be poor.  HJ, 24/Aug/2011
         scalarField nfc =
-            n & (overlapGgiPolyPatch_.reconFaceCellCentres() - Cf());
+            mag(n & (overlapGgiPolyPatch_.reconFaceCellCentres() - Cf()));
-        w = nfc/((n & (Cf() - Cn())) + nfc);
+
+        w = nfc/(mag(n & (Cf() - Cn())) + nfc);
     }
     else
     {
@@ -71,7 +78,10 @@ void Foam::overlapGgiFvPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     if (overlapGgiPolyPatch_.master())
     {
-        dc = 1.0/max(nf() & delta(), 0.05*mag(delta()));
+        // Stabilised form for bad meshes.  HJ, 24/Aug/2011
+        vectorField d = delta();
+
+        dc = 1.0/max(nf() & d, 0.05*mag(d));
     }
     else
     {

src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.C

@@ -46,18 +46,26 @@ void processorFvPatch::makeWeights(scalarField& w) const
     if (Pstream::parRun())
     {
         // The face normals point in the opposite direction on the other side
-        scalarField neighbFaceCentresCn
-        (
-            (
-                procPolyPatch_.neighbFaceAreas()
-               /(mag(procPolyPatch_.neighbFaceAreas()) + VSMALL)
-            )
-          & (
-              procPolyPatch_.neighbFaceCentres()
-            - procPolyPatch_.neighbFaceCellCentres())
-        );
 
-        w = neighbFaceCentresCn/((nf() & fvPatch::delta())
+        // Note: mag in the dot-product.
+        // For all valid meshes, the non-orthogonality will be less that
+        // 90 deg and the dot-product will be positive.  For invalid
+        // meshes (d & s <= 0), this will stabilise the calculation
+        // but the result will be poor.  HJ, 24/Aug/2011
+        scalarField neighbFaceCentresCn =
+            mag
+            (
+                (
+                    procPolyPatch_.neighbFaceAreas()/
+                    (mag(procPolyPatch_.neighbFaceAreas()) + VSMALL)
+                )
+              & (
+                    procPolyPatch_.neighbFaceCentres()
+                  - procPolyPatch_.neighbFaceCellCentres()
+                )
+            );
+
+        w = neighbFaceCentresCn/(mag(nf() & fvPatch::delta())
             + neighbFaceCentresCn);
     }
     else
@@ -71,7 +79,10 @@ void processorFvPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     if (Pstream::parRun())
     {
-        dc = (1.0 - weights())/(nf() & fvPatch::delta());
+        vectorField d = delta();
+
+        // Stabilised form for bad meshes
+        dc = 1.0/max((nf() & d), 0.05*mag(d));
     }
     else
     {

src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C

@@ -53,9 +53,16 @@ void Foam::regionCoupleFvPatch::makeWeights(scalarField& w) const
     if (rcPolyPatch_.attached())
     {
         vectorField n = nf();
-        scalarField nfc = n & (rcPolyPatch_.reconFaceCellCentres() - Cf());
 
-        w = nfc/((n & (Cf() - Cn())) + nfc);
+        // Note: mag in the dot-product.
+        // For all valid meshes, the non-orthogonality will be less that
+        // 90 deg and the dot-product will be positive.  For invalid
+        // meshes (d & s <= 0), this will stabilise the calculation
+        // but the result will be poor.  HJ, 24/Aug/2011
+        scalarField nfc =
+            mag(n & (rcPolyPatch_.reconFaceCellCentres() - Cf()));
+
+        w = nfc/(mag(n & (Cf() - Cn())) + nfc);
     }
     else
     {
@@ -69,7 +76,10 @@ void Foam::regionCoupleFvPatch::makeDeltaCoeffs(scalarField& dc) const
 {
     if (rcPolyPatch_.attached())
     {
-        dc = (1.0 - weights())/(nf() & fvPatch::delta());
+        // Stabilised form for bad meshes.  HJ, 24/Aug/2011
+        vectorField d = delta();
+
+        dc = 1.0/max(nf() & d, 0.05*mag(d));
     }
     else
     {