From a909be7d9b8c2ee5aa96e09b9eccc4bc6bd179e5 Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Fri, 4 Feb 2011 17:57:52 +0000
Subject: [PATCH] ENH: pairPatchAgglomeratnio

---
 src/fvAgglomerationMethods/Allwmake           |   2 +
 .../pairPatchAgglomeration/Make/files         |   3 +
 .../pairPatchAgglomeration/Make/options       |  10 +
 .../pairPatchAgglomeration.C                  | 528 ++++++++++++++++++
 .../pairPatchAgglomeration.H                  | 221 ++++++++
 .../pairPatchAgglomerationTemplates.C         |  79 +++
 6 files changed, 843 insertions(+)
 create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/Make/files
 create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
 create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
 create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.H
 create mode 100644 src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomerationTemplates.C

diff --git a/src/fvAgglomerationMethods/Allwmake b/src/fvAgglomerationMethods/Allwmake
index 61d951c8c6..8328d915fc 100755
--- a/src/fvAgglomerationMethods/Allwmake
+++ b/src/fvAgglomerationMethods/Allwmake
@@ -10,4 +10,6 @@ then
     wmake $makeOption MGridGenGamgAgglomeration
 fi
 
+wmake libso pairPatchAgglomeration
+
 # ----------------------------------------------------------------- end-of-file
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/files b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/files
new file mode 100644
index 0000000000..275c4795f3
--- /dev/null
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/files
@@ -0,0 +1,3 @@
+pairPatchAgglomeration.C
+
+LIB = $(FOAM_LIBBIN)/libpairPatchAgglomeration
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
new file mode 100644
index 0000000000..6dddff38af
--- /dev/null
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/Make/options
@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/OpenFOAM/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lOpenFOAM \
+    -lmeshTools
+
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
new file mode 100644
index 0000000000..d4657ef5e8
--- /dev/null
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.C
@@ -0,0 +1,528 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pairPatchAgglomeration.H"
+//#include "OFstream.H"
+#include "meshTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::pairPatchAgglomeration::compactLevels(const label nCreatedLevels)
+{
+    nFaces_.setSize(nCreatedLevels);
+    restrictAddressing_.setSize(nCreatedLevels);
+    patchLevels_.setSize(nCreatedLevels);
+}
+
+
+bool Foam::pairPatchAgglomeration::continueAgglomerating
+(
+    const label nCoarseFaces
+)
+{
+    // Check the need for further agglomeration on all processors
+    bool contAgg = nCoarseFaces >= nFacesInCoarsestLevel_;
+    reduce(contAgg, andOp<bool>());
+    return contAgg;
+}
+
+
+void Foam::pairPatchAgglomeration::setBasedEdgeWeights()
+{
+    const bPatch& coarsePatch = patchLevels_[0];
+    forAll(coarsePatch.edges(), i)
+    {
+        if(coarsePatch.isInternalEdge(i))
+        {
+            scalar edgeLength =
+                coarsePatch.edges()[i].mag(coarsePatch.localPoints());
+
+            const labelList& eFaces = coarsePatch.edgeFaces()[i];
+
+            if (eFaces.size() == 2)
+            {
+                scalar cosI =
+                        coarsePatch.faceNormals()[eFaces[0]] &
+                        coarsePatch.faceNormals()[eFaces[1]];
+
+                const edge edgeCommon = edge(eFaces[0], eFaces[1]);
+
+                if(facePairWeight_.found(edgeCommon))
+                {
+                    facePairWeight_[edgeCommon] += edgeLength;
+                }
+                else
+                {
+                    facePairWeight_.insert(edgeCommon, edgeLength);
+                }
+
+                if
+                (
+                    cosI <
+                    Foam::cos
+                    (
+                        featureAngle_*constant::mathematical::pi/180.0
+                    )
+                )
+                {
+                    facePairWeight_[edgeCommon] = -1.0;
+                }
+            }
+            else if (eFaces.size() == 3)
+            {
+                facePairWeight_.insert(edge(eFaces[1], eFaces[0]), -1.0);
+                facePairWeight_.insert(edge(eFaces[2], eFaces[0]), -1.0);
+                facePairWeight_.insert(edge(eFaces[1], eFaces[2]), -1.0);
+            }
+        }
+    }
+}
+
+
+void Foam::pairPatchAgglomeration::setEdgeWeights
+(
+    const label fineLevelIndex
+)
+{
+
+    const bPatch& coarsePatch = patchLevels_[fineLevelIndex];
+
+    const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
+    const label nCoarseI =  max(fineToCoarse)+1;
+    labelListList coarseToFine
+    (
+        invertOneToMany(nCoarseI, fineToCoarse)
+    );
+
+    HashSet<edge, Hash<edge> > fineFeaturedFaces(coarsePatch.nEdges()/10);
+
+    // Map fine faces with featured edge into coarse faces
+    forAllConstIter(EdgeMap<scalar>, facePairWeight_, iter)
+    {
+        if (iter() == -1.0)
+        {
+            const edge e = iter.key();
+            const edge edgeFeatured
+            (
+                fineToCoarse[e[0]],
+                fineToCoarse[e[1]]
+            );
+            fineFeaturedFaces.insert(edgeFeatured);
+        }
+    }
+
+    // Clean old weitghs
+    facePairWeight_.clear();
+    facePairWeight_.resize(coarsePatch.nEdges());
+
+    forAll(coarsePatch.edges(), i)
+    {
+        if(coarsePatch.isInternalEdge(i))
+        {
+            scalar edgeLength =
+                coarsePatch.edges()[i].mag(coarsePatch.localPoints());
+
+            const labelList& eFaces = coarsePatch.edgeFaces()[i];
+
+            if (eFaces.size() == 2)
+            {
+                const edge edgeCommon = edge(eFaces[0], eFaces[1]);
+                if(facePairWeight_.found(edgeCommon))
+                {
+                    facePairWeight_[edgeCommon] += edgeLength;
+                }
+                else
+                {
+                    facePairWeight_.insert(edgeCommon, edgeLength);
+                }
+                // If the fine 'pair' faces was featured edge so it is
+                // the coarse 'pair'
+                if (fineFeaturedFaces.found(edgeCommon))
+                {
+                    facePairWeight_[edgeCommon] = -1.0;
+                }
+            }
+            else if (eFaces.size() == 3)
+            {
+                facePairWeight_.insert(edge(eFaces[1], eFaces[0]), -1.0);
+                facePairWeight_.insert(edge(eFaces[2], eFaces[0]), -1.0);
+                facePairWeight_.insert(edge(eFaces[1], eFaces[2]), -1.0);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pairPatchAgglomeration::pairPatchAgglomeration
+(
+    const polyPatch& patch,
+    const dictionary& controlDict,
+    const bool additionalWeights
+)
+:
+    mergeLevels_
+    (
+        controlDict.lookupOrDefault<label>("mergeLevels", 2)
+    ),
+    maxLevels_(50),
+    nFacesInCoarsestLevel_
+    (
+        readLabel(controlDict.lookup("nFacesInCoarsestLevel"))
+    ),
+    featureAngle_
+    (
+        controlDict.lookupOrDefault<scalar>("featureAngle", 0)
+    ),
+    nFaces_(maxLevels_),
+    restrictAddressing_(maxLevels_),
+    restrictTopBottomAddressing_(identity(patch.size())),
+    patchLevels_(maxLevels_),
+    facePairWeight_(patch.size())
+{
+    // Set base fine patch
+    patchLevels_.set
+    (
+        0,
+        new bPatch
+        (
+            patch.localFaces(),
+            patch.localPoints()
+        )
+    );
+
+    // Set number of faces for the base patch
+    nFaces_[0] = patch.size();
+
+    // Set edge weights for level 0
+    setBasedEdgeWeights();
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::pairPatchAgglomeration::~pairPatchAgglomeration()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::bPatch& Foam::pairPatchAgglomeration::patchLevel
+(
+    const label i
+) const
+{
+    return patchLevels_[i];
+}
+
+
+void Foam::pairPatchAgglomeration::mapBaseToTopAgglom
+(
+    const label fineLevelIndex
+)
+{
+    const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
+    forAll (restrictTopBottomAddressing_, i)
+    {
+        restrictTopBottomAddressing_[i] =
+            fineToCoarse[restrictTopBottomAddressing_[i]];
+    }
+}
+
+
+bool Foam::pairPatchAgglomeration::agglomeratePatch
+(
+    const bPatch& patch,
+    const labelList& fineToCoarse,
+    const label fineLevelIndex
+)
+{
+    if (min(fineToCoarse) == -1)
+    {
+        FatalErrorIn("pairPatchAgglomeration::agglomeratePatch")
+            << "min(fineToCoarse) == -1" << exit(FatalError);
+    }
+
+    if (fineToCoarse.size() != patch.size())
+    {
+        FatalErrorIn
+        (
+            "pairPatchAgglomeration::agglomeratePatch"
+            "(const label fineLevelIndex)"
+        )   << "restrict map does not correspond to fine level. " << endl
+            << " Sizes: restrictMap: " << fineToCoarse.size()
+            << " nEqns: " << patch.size()
+            << abort(FatalError);
+    }
+
+    const label nCoarseI =  max(fineToCoarse)+1;
+    List<face> patchFaces(nCoarseI);
+
+    // Patch faces per agglomeration
+    labelListList coarseToFine
+    (
+        invertOneToMany(nCoarseI, fineToCoarse)
+    );
+
+    for (label coarseI = 0; coarseI < nCoarseI; coarseI++)
+    {
+        const labelList& fineFaces = coarseToFine[coarseI];
+
+        // Construct single face
+        indirectPrimitivePatch upp
+        (
+            IndirectList<face>(patch, fineFaces),
+            patch.points()
+        );
+
+        if (upp.edgeLoops().size() != 1)
+        {
+            if (fineFaces.size() == 2)
+            {
+                const edge e(fineFaces[0], fineFaces[1]);
+                facePairWeight_[e] = -1.0;
+            }
+            else if (fineFaces.size() == 3)
+            {
+                const edge e(fineFaces[0], fineFaces[1]);
+                const edge e1(fineFaces[0], fineFaces[2]);
+                const edge e2(fineFaces[2], fineFaces[1]);
+                facePairWeight_[e] = -1.0;
+                facePairWeight_[e1] = -1.0;
+                facePairWeight_[e2] = -1.0;
+            }
+
+            return false;
+        }
+
+        patchFaces[coarseI] = face
+        (
+            renumber
+            (
+                upp.meshPoints(),
+                upp.edgeLoops()[0]
+            )
+        );
+    }
+
+    patchLevels_.set
+    (
+        fineLevelIndex,
+        new bPatch
+        (
+            SubList<face>(patchFaces, nCoarseI, 0),
+            patch.points()
+        )
+    );
+    return true;
+}
+
+
+void Foam::pairPatchAgglomeration:: agglomerate()
+{
+    label nPairLevels = 0;
+    label nCreatedLevels = 1; //0 level is the base patch
+
+    while (nCreatedLevels < maxLevels_)
+    {
+        label nCoarseCells = -1;
+
+        const bPatch& patch = patchLevels_[nCreatedLevels - 1];
+        tmp<labelField> finalAgglomPtr(new labelField(patch.size()));
+
+        bool agglomOK = false;
+        while (!agglomOK)
+        {
+            finalAgglomPtr = agglomerateOneLevel
+            (
+                nCoarseCells,
+                patch
+            );
+
+            agglomOK = agglomeratePatch
+            (
+                patch,
+                finalAgglomPtr,
+                nCreatedLevels
+            );
+        }
+
+        nFaces_[nCreatedLevels] = nCoarseCells;
+        restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
+        mapBaseToTopAgglom(nCreatedLevels);
+
+        if(!continueAgglomerating(nCoarseCells))
+        {
+            break;
+        }
+
+        setEdgeWeights(nCreatedLevels);
+
+        if (nPairLevels % mergeLevels_)
+        {
+            combineLevels(nCreatedLevels);
+        }
+        else
+        {
+            nCreatedLevels++;
+        }
+
+        nPairLevels++;
+    }
+}
+
+
+Foam::tmp<Foam::labelField> Foam::pairPatchAgglomeration::agglomerateOneLevel
+(
+    label& nCoarseCells,
+    const bPatch& patch
+)
+{
+    const label nFineFaces = patch.size();
+
+    tmp<labelField> tcoarseCellMap(new labelField(nFineFaces, -1));
+    labelField& coarseCellMap = tcoarseCellMap();
+
+    const labelListList& faceFaces = patch.faceFaces();
+
+    nCoarseCells = 0;
+
+    forAll (faceFaces, facei)
+    {
+        const labelList& fFaces = faceFaces[facei];
+
+        if (coarseCellMap[facei] < 0)
+        {
+            label matchFaceNo = -1;
+            label matchFaceNeibNo = -1;
+            scalar maxFaceWeight = -GREAT;
+
+            // check faces to find ungrouped neighbour with largest face weight
+            forAll(fFaces, i)
+            {
+                label faceNeig = fFaces[i];
+                const edge edgeCommon = edge(facei, faceNeig);
+                if
+                (
+                    facePairWeight_[edgeCommon] > maxFaceWeight
+                    && coarseCellMap[faceNeig] < 0
+                    && facePairWeight_[edgeCommon] != -1.0
+                )
+                {
+                    // Match found. Pick up all the necessary data
+                    matchFaceNo = facei;
+                    matchFaceNeibNo = faceNeig;
+                    maxFaceWeight = facePairWeight_[edgeCommon];
+                }
+            }
+
+            if (matchFaceNo >= 0)
+            {
+                // Make a new group
+                coarseCellMap[matchFaceNo] = nCoarseCells;
+                coarseCellMap[matchFaceNeibNo] = nCoarseCells;
+                nCoarseCells++;
+            }
+            else
+            {
+                // No match. Find the best neighbouring cluster and
+                // put the cell there
+                label clusterMatchFaceNo = -1;
+                scalar clusterMaxFaceCoeff = -GREAT;
+
+                forAll(fFaces, i)
+                {
+                    label faceNeig = fFaces[i];
+                    const edge edgeCommon = edge(facei, faceNeig);
+                    if
+                    (
+                        facePairWeight_[edgeCommon] > clusterMaxFaceCoeff
+                        && facePairWeight_[edgeCommon] != -1.0
+                        && coarseCellMap[faceNeig] > 0
+                    )
+                    {
+                        clusterMatchFaceNo = faceNeig;
+                        clusterMaxFaceCoeff = facePairWeight_[edgeCommon];
+                    }
+                }
+
+                if (clusterMatchFaceNo >= 0)
+                {
+                    // Add the cell to the best cluster
+                    coarseCellMap[facei] = coarseCellMap[clusterMatchFaceNo];
+                }
+                else
+                {// if not create single-cell "clusters" for each
+                    coarseCellMap[facei] = nCoarseCells;
+                    nCoarseCells ++;
+                }
+            }
+        }
+    }
+
+    // Check that all faces are part of clusters,
+
+    for (label facei=0; facei<nFineFaces; facei++)
+    {
+        if (coarseCellMap[facei] < 0)
+        {
+            FatalErrorIn
+            (
+                "pairPatchAgglomeration::agglomerateOneLevel "
+                "(label&, const bPatch&) "
+            ) << " face " << facei
+            << " is not part of a cluster"
+            << exit(FatalError);
+            //coarseCellMap[facei] = nCoarseCells;
+            //nCoarseCells++;
+        }
+    }
+
+    return tcoarseCellMap;
+}
+
+void Foam::pairPatchAgglomeration::combineLevels(const label curLevel)
+{
+    label prevLevel = curLevel - 1;
+
+    // Set the previous level nCells to the current
+    nFaces_[prevLevel] = nFaces_[curLevel];
+
+    // Map the restrictAddressing from the coarser level into the previous
+    // finer level
+
+    const labelList& curResAddr = restrictAddressing_[curLevel];
+    labelList& prevResAddr = restrictAddressing_[prevLevel];
+
+    forAll(prevResAddr, i)
+    {
+        prevResAddr[i] = curResAddr[prevResAddr[i]];
+    }
+
+    // Delete the restrictAddressing for the coarser level
+    restrictAddressing_.set(curLevel, NULL);
+
+    patchLevels_.set(prevLevel, patchLevels_.set(curLevel, NULL));
+}
+// ************************************************************************* //
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.H b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.H
new file mode 100644
index 0000000000..ed66673bd7
--- /dev/null
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomeration.H
@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::pairPatchAgglomeration
+
+Description
+    Primitive patch pair agglomerate method.
+
+SourceFiles
+    pairPatchAgglomeration.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pairPatchAgglomeration_H
+#define pairPatchAgglomeration_H
+
+#include "mathematicalConstants.H"
+#include "polyPatch.H"
+#include "indirectPrimitivePatch.H"
+#include "List.H"
+#include "EdgeMap.H"
+
+namespace Foam
+{
+    typedef PrimitivePatch<face, List, const pointField> bPatch;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class pairPatchAgglomeration Declaration
+\*---------------------------------------------------------------------------*/
+
+class pairPatchAgglomeration
+{
+protected:
+
+    // Protected data
+
+        //- Number of levels to merge, 1 = don't merge, 2 = merge pairs etc.
+        label mergeLevels_;
+
+        //- Max number of levels
+        label maxLevels_;
+
+        //- Number of faces in coarsest level
+        label nFacesInCoarsestLevel_;
+
+        //- Feature angle
+        scalar featureAngle_;
+
+        //- The number of faces in each level
+        labelList nFaces_;
+
+        //- Cell restriction addressing array.
+        //  Maps from the finer to coarse level
+        PtrList<labelField> restrictAddressing_;
+
+        //- Maps from finest to coarsest
+        labelList restrictTopBottomAddressing_;
+
+        //- Hierarchy of patch addressing
+        PtrList<bPatch> patchLevels_;
+
+        //- Edge weights
+        EdgeMap<scalar> facePairWeight_;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Assemble coarse patch
+        bool agglomeratePatch
+        (
+            const bPatch& patch,
+            const labelList& fineToCoarse,
+            const label fineLevelIndex
+        );
+
+        //- Agglomerate one level
+        tmp<labelField> agglomerateOneLevel
+        (
+            label& nCoarseCells,
+            const bPatch& patchLevel
+        );
+
+        //- Combine leves
+        void combineLevels(const label curLevel);
+
+        //- Shrink the number of levels to that specified
+        void compactLevels(const label fineLevelIndex);
+
+        //- Check the need for further agglomeration
+        bool continueAgglomerating(const label fineLevelIndex);
+
+        //- Set edge weights
+        void setEdgeWeights(const label indexLevel);
+
+        //- Set base patch edge weights
+        void setBasedEdgeWeights();
+
+        //- Maps current level with base patch
+        void mapBaseToTopAgglom(const label fineLevelIndex);
+
+        //- Disallow default bitwise copy construct
+        pairPatchAgglomeration(const pairPatchAgglomeration&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const pairPatchAgglomeration&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given mesh and controls
+        pairPatchAgglomeration
+        (
+            const polyPatch& patch,
+            const dictionary& controlDict,
+            const bool additionalWeights = false
+        );
+
+
+    // Destructor
+        ~pairPatchAgglomeration();
+
+
+    // Member Functions
+
+        //- Agglomerate patch
+            void agglomerate();
+
+
+        // Access
+
+            //- Return size
+            label size() const
+            {
+                return patchLevels_.size();
+            }
+
+            //- Return restriction from top level to bottom level
+            const labelList& restrictTopBottomAddressing() const
+            {
+                return restrictTopBottomAddressing_;
+            }
+
+            //- Return primitivePatch of given level
+            const bPatch& patchLevel(const label leveli) const;
+
+            //- Return cell restrict addressing of given level
+            const labelField& restrictAddressing(const label leveli) const
+            {
+                return restrictAddressing_[leveli];
+            }
+
+
+        // Restriction and prolongation
+
+            //- Restrict (integrate by summation) cell field
+            template<class Type>
+            void restrictField
+            (
+                Field<Type>& cf,
+                const Field<Type>& ff,
+                const label fineLevelIndex
+            ) const;
+
+            //- Prolong (interpolate by injection) cell field
+            template<class Type>
+            void prolongField
+            (
+                Field<Type>& ff,
+                const Field<Type>& cf,
+                const label coarseLevelIndex
+            ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "pairPatchAgglomerationTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomerationTemplates.C b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomerationTemplates.C
new file mode 100644
index 0000000000..0ac62aa3ad
--- /dev/null
+++ b/src/fvAgglomerationMethods/pairPatchAgglomeration/pairPatchAgglomerationTemplates.C
@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pairPatchAgglomeration.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::pairPatchAgglomeration::restrictField
+(
+    Field<Type>& cf,
+    const Field<Type>& ff,
+    const label fineLevelIndex
+) const
+{
+    const labelList& fineToCoarse = restrictAddressing_[fineLevelIndex];
+
+    if (ff.size() != fineToCoarse.size())
+    {
+        FatalErrorIn
+        (
+            "void pairPatchAgglomeration::restrictField"
+            "(Field<Type>& cf, const Field<Type>& ff, "
+            "const label fineLevelIndex) const"
+        )   << "field does not correspond to level " << fineLevelIndex
+            << " sizes: field = " << ff.size()
+            << " level = " << fineToCoarse.size()
+            << abort(FatalError);
+    }
+
+    cf = pTraits<Type>::zero;
+
+    forAll(ff, i)
+    {
+        cf[fineToCoarse[i]] += ff[i];
+    }
+}
+
+
+template<class Type>
+void Foam::pairPatchAgglomeration::prolongField
+(
+    Field<Type>& ff,
+    const Field<Type>& cf,
+    const label coarseLevelIndex
+) const
+{
+    const labelList& fineToCoarse = restrictAddressing_[coarseLevelIndex];
+
+    forAll(fineToCoarse, i)
+    {
+        ff[i] = cf[fineToCoarse[i]];
+    }
+}
+
+
+// ************************************************************************* //