ENH: assemblyFaceAreaPair: revert to faceAreaPair. Fixes #2

src/regionCoupled/assemblyFaceAreaPairGAMGAgglomeration/assemblyFaceAreaPairGAMGAgglomeration.C

@@ -27,7 +27,7 @@ License
 #include "lduMatrix.H"
 #include "addToRunTimeSelectionTable.H"
 #include "lduPrimitiveMeshAssembly.H"
-#include "fvMatrix.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -55,10 +55,7 @@ namespace Foam
 
 Foam::assemblyFaceAreaPairGAMGAgglomeration::
 ~assemblyFaceAreaPairGAMGAgglomeration()
-{
-    deleteDemandDrivenData(faceAreasPtr_);
-    deleteDemandDrivenData(cellVolumesPtr_);
-}
+{}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -69,123 +66,134 @@ assemblyFaceAreaPairGAMGAgglomeration
     const dictionary& controlDict
 )
 :
-    pairGAMGAgglomeration(matrix.mesh(), controlDict),
-    faceAreasPtr_(nullptr),
-    cellVolumesPtr_(nullptr)
+    pairGAMGAgglomeration(matrix.mesh(), controlDict)
 {
 
     const lduMesh& ldumesh = matrix.mesh();
 
-    if (!isA<lduPrimitiveMeshAssembly>(ldumesh))
+    if (isA<lduPrimitiveMeshAssembly>(ldumesh))
     {
-        FatalErrorInFunction
-            << "assemblyFaceAreaPairGAMGAgglomeration requested but "
-            << " lduMesh is not  lduPrimitiveMeshAssembly." << nl
-            << " Change agglomerator type to faceAreaPair "
-            << abort(FatalError);
-    }
-    const lduPrimitiveMeshAssembly& mesh =
-        refCast<const lduPrimitiveMeshAssembly>(ldumesh);
-
-    faceAreasPtr_ = new vectorField(mesh.lduAddr().upperAddr().size(), Zero);
-    vectorField& faceAreas = *faceAreasPtr_;
+        const lduPrimitiveMeshAssembly& mesh =
+            refCast<const lduPrimitiveMeshAssembly>(ldumesh);
 
-    cellVolumesPtr_ = new scalarField(mesh.lduAddr().size(), Zero);
-    scalarField& cellVolumes = *cellVolumesPtr_;
+        vectorField faceAreas(mesh.lduAddr().upperAddr().size(), Zero);
+        //scalarField cellVolumes(mesh.lduAddr().size(), Zero);
 
-    const labelListList& faceMap = mesh.faceMap();
+        const labelListList& faceMap = mesh.faceMap();
 
-    for (label i=0; i < mesh.meshes().size(); ++i)
-    {
-        const fvMesh& m = refCast<const fvMesh>(mesh.meshes()[i]);
-        const labelList& subFaceMap = faceMap[i];
-        const vectorField& areas = m.Sf();
-
-        forAll(subFaceMap, facei)
+        for (label i=0; i < mesh.meshes().size(); ++i)
         {
-            faceAreas[subFaceMap[facei]] = areas[facei];
-        }
+            const fvMesh& m = refCast<const fvMesh>(mesh.meshes()[i]);
+            const labelList& subFaceMap = faceMap[i];
+            const vectorField& areas = m.Sf().internalField();
 
-        const polyBoundaryMesh& patches = m.boundaryMesh();
+            forAll(subFaceMap, facei)
+            {
+                faceAreas[subFaceMap[facei]] = areas[facei];
+            }
 
-        // Fill faceAreas for new faces
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-            label globalPatchID = mesh.patchMap()[i][patchI];
+            const polyBoundaryMesh& patches = m.boundaryMesh();
 
-            if (globalPatchID == -1)
+            // Fill faceAreas for new faces
+            forAll(patches, patchI)
             {
-                //const lduInterface& inter = m.interfaces()[patchI];
+                const polyPatch& pp = patches[patchI];
+                label globalPatchID = mesh.patchMap()[i][patchI];
 
-                if (pp.masterImplicit())
+                if (globalPatchID == -1)
                 {
-                    const vectorField& sf = m.boundary()[patchI].Sf();
+                    //const lduInterface& inter = m.interfaces()[patchI];
 
-                    if (isA<cyclicAMIPolyPatch>(pp))
+                    if (pp.masterImplicit())
                     {
-                        const cyclicAMIPolyPatch& mpp =
-                                refCast<const cyclicAMIPolyPatch>(patches[patchI]);
-
-                        const scalarListList& srcWeight = mpp.AMI().srcWeights();
+                        const vectorField& sf = m.boundary()[patchI].Sf();
 
-                        label subFaceI = 0;
-                        forAll(pp.faceCells(), faceI)
+                        if (isA<cyclicAMIPolyPatch>(pp))
                         {
-                            const scalarList& w = srcWeight[faceI];
+                            const cyclicAMIPolyPatch& mpp =
+                                    refCast<const cyclicAMIPolyPatch>(pp);
+
+                            const scalarListList& srcWeight =
+                                mpp.AMI().srcWeights();
 
-                            for(label j=0; j<w.size(); j++)
+                            label subFaceI = 0;
+                            forAll(pp.faceCells(), faceI)
                             {
-                                const label globalFaceI =
-                                    mesh.faceBoundMap()[i][patchI][subFaceI];
+                                const scalarList& w = srcWeight[faceI];
 
-                                if (globalFaceI != -1)
+                                for(label j=0; j<w.size(); j++)
                                 {
-                                    faceAreas[globalFaceI] = w[j]*sf[faceI];
+                                    const label globalFaceI =
+                                        mesh.faceBoundMap()[i][patchI][subFaceI];
+
+                                    if (globalFaceI != -1)
+                                    {
+                                        faceAreas[globalFaceI] = w[j]*sf[faceI];
+                                    }
+                                    subFaceI++;
                                 }
-                                subFaceI++;
                             }
                         }
-                    }
-                    else
-                    {
-                        forAll(pp.faceCells(), faceI)
+                        else
                         {
-                            const label globalFaceI =
-                                mesh.faceBoundMap()[i][patchI][faceI];
-
-                            if (globalFaceI != -1)
+                            forAll(pp.faceCells(), faceI)
                             {
-                                faceAreas[globalFaceI] = sf[faceI];
+                                const label globalFaceI =
+                                    mesh.faceBoundMap()[i][patchI][faceI];
+
+                                if (globalFaceI != -1)
+                                {
+                                    faceAreas[globalFaceI] = sf[faceI];
+                                }
                             }
                         }
                     }
                 }
             }
-        }
-
-        // Fill cellVolumes
-        const scalarField& vol = m.V();
-        const label cellOffset = mesh.cellOffsets()[i];
 
-        forAll(vol, cellI)
-        {
-            cellVolumes[cellOffset + cellI] = vol[cellI];
+            //// Fill cellVolumes
+            //const scalarField& vol = m.V();
+            //const label cellOffset = mesh.cellOffsets()[i];
+            //
+            //forAll(vol, cellI)
+            //{
+            //    cellVolumes[cellOffset + cellI] = vol[cellI];
+            //}
         }
+
+        agglomerate
+        (
+            mesh,
+            mag
+            (
+                cmptMultiply
+                (
+                    faceAreas/sqrt(mag(faceAreas)),
+                    vector(1, 1.01, 1.02)
+                )
+            )
+        );
     }
+    else
+    {
+        // Revert to faceAreaPairGAMGAgglomeration
+        const fvMesh& fvmesh = refCast<const fvMesh>(matrix.mesh());
 
-    agglomerate
-    (
-        mesh,
-        mag
+        agglomerate
         (
-            cmptMultiply
+            matrix.mesh(),
+            mag
             (
-                *faceAreasPtr_/sqrt(mag(*faceAreasPtr_)),
-                vector(1, 1.01, 1.02)
+                cmptMultiply
+                (
+                    fvmesh.Sf().primitiveField()
+                   /sqrt(fvmesh.magSf().primitiveField()),
+                    vector(1, 1.01, 1.02)
+                    //vector::one
+                )
             )
-        )
-    );
+        );
+    }
 }
 
 

src/regionCoupled/assemblyFaceAreaPairGAMGAgglomeration/assemblyFaceAreaPairGAMGAgglomeration.H

@@ -50,14 +50,6 @@ class assemblyFaceAreaPairGAMGAgglomeration
 :
     public pairGAMGAgglomeration
 {
-private:
-
-        //- Face areas of the assembled mesh
-        vectorField* faceAreasPtr_;
-
-        //- Cell volumes of the assembled mesh
-        scalarField* cellVolumesPtr_;
-
 public:
 
     //- Runtime type information