Updates to fvcReconstruct and Overset include files. Vuko Vukcevic

applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C

@@ -25,7 +25,8 @@ Description
     Extrude mesh from existing patch (by default outwards facing normals;
     optional flips faces) or from patch read from file.
 
-    Note: Merges close points so be careful.
+    Note: Merges close points so be careful. This can be controlled with the
+    optional mergeTolerance parameter (1e-4) by default.
 
     Type of extrusion prescribed by run-time selectable model.
 
@@ -181,7 +182,8 @@ int main(int argc, char *argv[])
     const vector span = bb.span();
 
     // Read merge tolerance
-    const scalar mergeTolerance = readScalar(dict.lookup("mergeTolerance"));
+    const scalar mergeTolerance =
+        dict.lookupOrDefault<scalar>("mergeTolerance", 1e-4);
     const scalar mergeDim = mergeTolerance*bb.minDim();
 
     Info<< "Mesh bounding box : " << bb << nl

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C

@@ -459,7 +459,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 // cell with multiple faces, we need to collect neighbour cell
                 // indices on the other side and possibly swap the order of
                 // adding faces. The same "swapping" of insertion order needs to
-                // happend on the slave side, but now we sort on the remote
+                // happen on the slave side, but now we sort on the remote
                 // (master) data and not on the local (slave) data as we do for
                 // master processor. VV, 16/Feb/2019.
 
@@ -588,7 +588,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                         // above, owner and neighbour proc are the same. In
                         // order to avoid using cellToProc list for slave
                         // processor (where ownCellI is actually found on the
-                        // other side), use curNbrPtc with patchFace indes
+                        // other side), use curNbrPtc with patchFace index
                         const label& ownerProc = curNbrPtc[patchFaceI];
 
                         // Add the patch face into the list in the correct

src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C

@@ -76,34 +76,66 @@ reconstruct
     GeometricField<GradType, fvPatchField, volMesh>& reconField =
         treconField();
 
-    // Note:
-    // 1) Reconstruction is only available in cell centres: there is no need
+    // Notes regarding boundary:
+    // 1. Reconstruction is only available in cell centres: there is no need
     //    to invert the tensor on the boundary
-    // 2) For boundaries, the only reconstructed data is the flux times
+    // 2. For boundaries, the only reconstructed data is the flux times
     //    normal. Based on this guess, boundary conditions can adjust
-    //    patch values
-    // HJ, 12/Aug/2011
+    //    patch values. HJ, 12/Aug/2011
 
-    GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
-        surfaceSum((mesh.Sf()/mesh.magSf())*ssf);
+    // Notes regarding procedure:
+    // 1. hinv inverse must be used to stabilise the inverse on bad meshes
+    //    but it gives strange failures because it unnecessarily avoids
+    //    performing ordinary inverse for meshes with reasonably sized
+    //    determinant (e.g. if SfSf/magSf is small). HJ, 19/Aug/2015
+    // 2. hinv has been stabilised now. HJ, 22/Mar/2019
+    // 3. But we still need to make sure that the determinant is not extremely
+    //    small, which may happen for extremely small meshes. We avoid this
+    //    issue by dividing the reconstruction equation with magSf^2 (instead of
+    //    magSf), which basically makes the dyadic tensor that we need to invert
+    //    dimensionless. VV, 13/Jun/2019
 
-    // Note: hinv inverse must be used to stabilise the inverse on bad meshes
-    // but it gives strange failures.  Fixed hinv.  HJ, 22/Mar/2019
-    // HJ, 19/Aug/2015
-    // Temporarily reverting hinv: Vuko Vukcevic, 6/Jun/2019
+    // Here's a short derivation in a Latex--like notation, where:
+    // - Sf is the surface area vector
+    // - uf is the face velocity factor (or field to be reconstructed)
+    // - Ff is the face flux
+    // - magSf is the magnitude of the surface area vector
+    // - uP is the velocity field in the cell centre
+    // - G is the surface area dyadic tensor
+    // - Fn is the vector representing surface sum of directional fluxes
+    // - \dprod is a dot product
+
+    // 1. Sf \dprod uf = Ff
+    //    Multiply Eq (1) with Sf/magSf^2
+    // 2. \frac{Sf Sf}{magSf^2} \dprod uf = \frac{Sf Ff}{magSf^2}
+    //    Sum Eq (2) over all the faces
+    // 3. \sum_f(\frac{Sf Sf}{magSf^2} \dprod uf) = \sum_f(\frac{Sf Ff}{magSf^2})
+    //    Assume first order extrapolation of uf, e.g. uP = uf
+    // 4. \sum_f(\frac{Sf Sf}{magSf^2}) \dprod uP) = \sum_f(\frac{Sf Ff}{magSf^2})
+    //    Use shorthand notation
+    // 5. G \dprod uP = Fn
+    // 6. uP = G^-1 \dprod Fn
+
+    // Fn -> fluxTimesNormal
+    // G  -> G
+
+    // Calculate sum of the directional fluxes
+    const surfaceScalarField magSfSqr = sqr(mesh.magSf());
+    const GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
+        surfaceSum((mesh.Sf()/magSfSqr)*ssf);
+
+    // Calculate the G tensor
+    const volSymmTensorField G = surfaceSum(sqr(mesh.Sf())/magSfSqr);
+
+    // Finally calculate the reconstructed field using hinv for stabilisation on
+    // really bad fvMesh bits (uses ordinary inverse most of the time, see
+    // tensor.C)
     reconField.internalField() =
-    (
-        // hinv
-        inv
-        (
-            surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField()
-        )
-      & fluxTimesNormal.internalField()
-    );
+        hinv(G.internalField()) & fluxTimesNormal.internalField();
 
+    // Boundary value update
     reconField.boundaryField() = fluxTimesNormal.boundaryField();
-
-    treconField().correctBoundaryConditions();
+    reconField.correctBoundaryConditions();
 
     return treconField;
 }

src/overset/oversetMesh/include/oversetAlphaCourantNo.H

@@ -38,21 +38,28 @@ scalar maxAlphaCo
 scalar alphaCoNum = 0.0;
 scalar meanAlphaCoNum = 0.0;
 
-surfaceScalarField alpha1f =
+const surfaceScalarField alpha1f =
     fvc::interpolate(min(max(alpha1, scalar(0)), scalar(1)));
 
-const dimensionedScalar alphaOffset("alphaOffset", dimless, dAlpha);
+const dimensionedScalar alphaOffset
+(
+    "alphaOffset",
+    dimless,
+    runTime.controlDict().lookupOrDefault("dAlpha", 0.01)
+);
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField magAlphaPhi
+    const oversetMesh& om = oversetMesh::New(mesh);
+
+    const surfaceScalarField magAlphaPhi
     (
         pos(alpha1f - alphaOffset)*
         pos(scalar(1) - alphaOffset - alpha1f)*
-        mag(faceOversetMask*phi)
+        mag(om.sGamma()*phi)
     );
 
-    surfaceScalarField SfUfbyDelta =
+    const surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*magAlphaPhi;
 
     const scalar deltaT = runTime.deltaT().value();

src/overset/oversetMesh/include/oversetCourantNo.H

@@ -36,9 +36,11 @@ scalar velMag = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField magPhi = mag(faceOversetMask*phi);
+    const oversetMesh& om = oversetMesh::New(mesh);
 
-    surfaceScalarField SfUfbyDelta =
+    const surfaceScalarField magPhi = mag(om.sGamma()*phi);
+
+    const surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*magPhi;
 
     CoNum = max(SfUfbyDelta/mesh.magSf())