ENH: overset: GAMG normalisation of hole cells

src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C

@@ -357,6 +357,71 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
 }
 
 
+Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
+(
+    const labelList& types,
+    const labelList& nbrTypes,
+    const scalarField& norm,
+    const scalarField& nbrNorm,
+    const label celli,
+    bitSet& isFront
+) const
+{
+    const labelList& own = faceOwner();
+    const labelList& nei = faceNeighbour();
+    const cell& cFaces = cells()[celli];
+
+    scalar avg = 0.0;
+    label n = 0;
+    label nFront = 0;
+    for (const label facei : cFaces)
+    {
+        if (isInternalFace(facei))
+        {
+            label nbrCelli = (own[facei] == celli ? nei[facei] : own[facei]);
+            if (norm[nbrCelli] == -GREAT)
+            {
+                // Invalid neighbour. Add to front
+                if (isFront.set(facei))
+                {
+                    nFront++;
+                }
+            }
+            else
+            {
+                // Valid neighbour. Add to average
+                avg += norm[nbrCelli];
+                n++;
+            }
+        }
+        else
+        {
+            if (nbrNorm[facei-nInternalFaces()] == -GREAT)
+            {
+                if (isFront.set(facei))
+                {
+                    nFront++;
+                }
+            }
+            else
+            {
+                avg += nbrNorm[facei-nInternalFaces()];
+                n++;
+            }
+        }
+    }
+
+    if (n > 0)
+    {
+        return avg/n;
+    }
+    else
+    {
+        return norm[celli];
+    }
+}
+
+
 void Foam::dynamicOversetFvMesh::writeAgglomeration
 (
     const GAMGAgglomeration& agglom
@@ -375,7 +440,8 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
                 this->time().timeName(),
                 *this,
                 IOobject::NO_READ,
-                IOobject::AUTO_WRITE
+                IOobject::NO_WRITE,
+                false
             ),
             *this,
             dimensionedScalar(dimless, Zero)
@@ -430,7 +496,8 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
                     this->time().timeName(),
                     *this,
                     IOobject::NO_READ,
-                    IOobject::AUTO_WRITE
+                    IOobject::NO_WRITE,
+                    false
                 ),
                 *this,
                 dimensionedScalar(dimless, Zero)

src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H

@@ -146,6 +146,17 @@ protected:
         template<class GeoField>
         static void correctCoupledBoundaryConditions(GeoField& fld);
 
+        //- Average norm of valid neighbours
+        scalar cellAverage
+        (
+            const labelList& types,
+            const labelList& nbrTypes,
+            const scalarField& norm,
+            const scalarField& nbrNorm,
+            const label celli,
+            bitSet& isFront
+        ) const;
+
         //- Debug: dump agglomeration
         void writeAgglomeration
         (

src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C

@@ -30,6 +30,7 @@ License
 #include "calculatedProcessorFvPatchField.H"
 #include "lduInterfaceFieldPtrsList.H"
 #include "processorFvPatch.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
@@ -171,18 +172,151 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
         }
     }
 
-    forAll(norm, celli)
+    // Count number of problematic cells
+    label nZeroDiag = 0;
+    for (const scalar n : norm)
     {
-        scalar& n = norm[celli];
-        if (mag(n) < SMALL)
+        if (magSqr(n) < sqr(SMALL))
         {
-            n = 1.0; //?
+            nZeroDiag++;
         }
-        else
+    }
+
+    reduce(nZeroDiag, sumOp<label>());
+
+    if (debug)
+    {
+        Pout<< "For field " << m.psi().name() << " have zero diagonals for "
+            << nZeroDiag << " cells" << endl;
+    }
+
+    if (nZeroDiag > 0)
+    {
+        // Walk out the norm across hole cells
+
+        const labelList& own = faceOwner();
+        const labelList& nei = faceNeighbour();
+        const cellCellStencilObject& overlap = Stencil::New(*this);
+        const labelUList& types = overlap.cellTypes();
+
+        label nHoles = 0;
+        scalarField extrapolatedNorm(norm);
+        forAll(types, celli)
         {
-            // Restore original diagonal
-            n = m.diag()[celli];
+            if (types[celli] == cellCellStencil::HOLE)
+            {
+                extrapolatedNorm[celli] = -GREAT;
+                nHoles++;
+            }
         }
+
+        bitSet isFront(nFaces());
+        for (label facei = 0; facei < nInternalFaces(); facei++)
+        {
+            label ownType = types[own[facei]];
+            label neiType = types[nei[facei]];
+            if
+            (
+                (ownType == cellCellStencil::HOLE)
+             != (neiType == cellCellStencil::HOLE)
+            )
+            {
+                isFront.set(facei);
+            }
+        }
+        labelList nbrTypes;
+        syncTools::swapBoundaryCellList(*this, types, nbrTypes);
+        for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+        {
+            label ownType = types[own[facei]];
+            label neiType = nbrTypes[facei-nInternalFaces()];
+            if
+            (
+                (ownType == cellCellStencil::HOLE)
+             != (neiType == cellCellStencil::HOLE)
+            )
+            {
+                isFront.set(facei);
+            }
+        }
+
+
+        while (true)
+        {
+            scalarField nbrNorm;
+            syncTools::swapBoundaryCellList(*this, extrapolatedNorm, nbrNorm);
+
+            bitSet newIsFront(nFaces());
+            scalarField newNorm(extrapolatedNorm);
+
+            label nChanged = 0;
+            for (const label facei : isFront)
+            {
+                if (extrapolatedNorm[own[facei]] == -GREAT)
+                {
+                    // Average owner cell, add faces to newFront
+                    newNorm[own[facei]] = cellAverage
+                    (
+                        types,
+                        nbrTypes,
+                        extrapolatedNorm,
+                        nbrNorm,
+                        own[facei],
+                        newIsFront
+                    );
+                    nChanged++;
+                }
+                if
+                (
+                    isInternalFace(facei)
+                 && extrapolatedNorm[nei[facei]] == -GREAT
+                )
+                {
+                    // Average nei cell, add faces to newFront
+                    newNorm[nei[facei]] = cellAverage
+                    (
+                        types,
+                        nbrTypes,
+                        extrapolatedNorm,
+                        nbrNorm,
+                        nei[facei],
+                        newIsFront
+                    );
+                    nChanged++;
+                }
+            }
+
+            reduce(nChanged, sumOp<label>());
+            if (nChanged == 0)
+            {
+                break;
+            }
+
+            // Transfer new front
+            extrapolatedNorm.transfer(newNorm);
+            isFront.transfer(newIsFront);
+            syncTools::syncFaceList(*this, isFront, maxEqOp<unsigned int>());
+        }
+
+
+        forAll(norm, celli)
+        {
+            scalar& n = norm[celli];
+            if (mag(n) < SMALL)
+            {
+                n = extrapolatedNorm[celli];
+            }
+            else
+            {
+                // Use original diagonal
+                n = m.diag()[celli];
+            }
+        }
+    }
+    else
+    {
+        // Use original diagonal
+        norm = m.diag();
     }
     return tnorm;
 }
@@ -545,6 +679,38 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     // Calculate stabilised diagonal as normalisation for interpolation
     const scalarField norm(normalisation(m));
 
+    if (debug)
+    {
+        volScalarField scale
+        (
+            IOobject
+            (
+                m.psi().name() + "_normalisation",
+                this->time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            *this,
+            dimensionedScalar(dimless, Zero)
+        );
+        scale.ref().field() = norm;
+        correctBoundaryConditions
+        <
+            volScalarField,
+            oversetFvPatchField<scalar>
+        >(scale.boundaryFieldRef(), false);
+        scale.write();
+
+        if (debug)
+        {
+            Pout<< "dynamicOversetFvMesh::solve() :"
+                << " writing matrix normalisation for field " << m.psi().name()
+                << " to " << scale.name() << endl;
+        }
+    }
+
 
     // Switch to extended addressing (requires mesh::update() having been
     // called)