ENH: have ptscotch

ptscotch - compiles into ptscotchDecomp. All thirdparty decompositionMethods
now moved out of decompositionMethods so add them explicitly to link line
for programs that need them (decomposePar, snappyHexMesh etc.)

applications/utilities/mesh/generation/snappyHexMesh/Make/options

@@ -11,6 +11,7 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -ldecompositionMethods \
+    -L$(FOAM_MPI_LIBBIN) -lptscotchDecomp -lparMetisDecomp -lmetisDecomp \
     -lmeshTools \
     -ldynamicMesh \
     -lautoMesh

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -347,7 +347,7 @@ int main(int argc, char *argv[])
             << "You have selected decomposition method "
             << decomposer.typeName
             << " which is not parallel aware." << endl
-            << "Please select one that is (hierarchical, parMetis)"
+            << "Please select one that is (hierarchical, ptscotch, parMetis)"
             << exit(FatalError);
     }
 

applications/utilities/parallelProcessing/decomposePar/Make/options

@@ -7,6 +7,6 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -lgenericPatchFields \
-    -ldecompositionMethods \
+    -ldecompositionMethods -lmetisDecomp -lscotchDecomp \
     -llagrangian \
     -lmeshTools

applications/utilities/parallelProcessing/redistributeMeshPar/Make/options

@@ -7,5 +7,6 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -ldecompositionMethods \
+    -L$(FOAM_MPI_LIBBIN) -lptscotchDecomp -lparMetisDecomp -lmetisDecomp \
     -lmeshTools \
     -ldynamicMesh

src/dummyThirdParty/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso scotchDecomp
+wmake libso ptscotchDecomp
 wmake libso metisDecomp
 wmake libso parMetisDecomp
 wmake libso MGridGen

src/dummyThirdParty/parMetisDecomp/dummyParMetisDecomp.C

@@ -163,23 +163,4 @@ Foam::labelList Foam::parMetisDecomp::decompose
 }
 
 
-void Foam::parMetisDecomp::calcMetisDistributedCSR
-(
-    const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    FatalErrorIn
-    (
-        "void parMetisDecomp::calcMetisDistributedCSR"
-        "("
-            "const polyMesh&, "
-            "List<int>&, "
-            "List<int>&"
-        ")"
-    )   << notImplementedMessage << exit(FatalError);
-}
-
-
 // ************************************************************************* //

src/dummyThirdParty/ptscotchDecomp/Make/files

@@ -0,0 +1,3 @@
+dummyPtscotchDecomp.C
+
+LIB = $(FOAM_LIBBIN)/dummy/libptscotchDecomp

src/dummyThirdParty/ptscotchDecomp/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(FOAM_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+    -I$(FOAM_SRC)/parallel/decompose/ptscotchDecomp/lnInclude
+
+LIB_LIBS =

src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "ptscotchDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+
+static const char* notImplementedMessage =
+"You are trying to use ptscotch but do not have the "
+"ptscotchDecomp library loaded."
+"\nThis message is from the dummy ptscotchDecomp stub library instead.\n"
+"\n"
+"Please install ptscotch and make sure that libptscotch.so is in your "
+"LD_LIBRARY_PATH.\n"
+"The ptscotchDecomp library can then be built in "
+"$FOAM_SRC/parallel/decompose/ptscotchDecomp\n";
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(ptscotchDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        ptscotchDecomp,
+        dictionaryMesh
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::ptscotchDecomp::check(const int retVal, const char* str)
+{}
+
+
+Foam::label Foam::ptscotchDecomp::decompose
+(
+    List<int>& adjncy,
+    List<int>& xadj,
+    const scalarField& cWeights,
+    List<int>& finalDecomp
+)
+{
+    FatalErrorIn
+    (
+        "label ptscotchDecomp::decompose"
+        "("
+            "const List<int>&, "
+            "const List<int>&, "
+            "const scalarField&, "
+            "List<int>&"
+        ")"
+    )   << notImplementedMessage << exit(FatalError);
+
+    return -1;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ptscotchDecomp::ptscotchDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    decompositionMethod(decompositionDict),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
+{
+    FatalErrorIn
+    (
+        "labelList ptscotchDecomp::decompose"
+        "("
+            "const pointField&, "
+            "const scalarField&"
+        ")"
+    )   << notImplementedMessage << exit(FatalError);
+
+    return labelList::null();
+}
+
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const labelList& agglom,
+    const pointField& agglomPoints,
+    const scalarField& pointWeights
+)
+{
+    FatalErrorIn
+    (
+        "labelList ptscotchDecomp::decompose"
+        "("
+            "const labelList&, "
+            "const pointField&, "
+            "const scalarField&"
+        ")"
+    )   << notImplementedMessage << exit(FatalError);
+
+    return labelList::null();
+}
+
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cellCentres,
+    const scalarField& cWeights
+)
+{
+    FatalErrorIn
+    (
+        "labelList ptscotchDecomp::decompose"
+        "("
+            "const labelListList&, "
+            "const pointField&, "
+            "const scalarField&"
+        ")"
+    )   << notImplementedMessage << exit(FatalError);
+
+    return labelList::null();
+}
+
+
+// ************************************************************************* //

src/dummyThirdParty/scotchDecomp/dummyScotchDecomp.C

@@ -22,87 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-    From scotch forum:
-
-    By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ]
-    2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order
-    not to be confused, you must have a clear view of how they are built.
-    Here are some rules:
-
-    1- Strategies are made up of "methods" which are combined by means of
-    "operators".
-
-    2- A method is of the form "m{param=value,param=value,...}", where "m"
-    is a single character (this is your first error: "f" is a method name,
-    not a parameter name).
-
-    3- There exist different sort of strategies : bipartitioning strategies,
-    mapping strategies, ordering strategies, which cannot be mixed. For
-    instance, you cannot build a bipartitioning strategy and feed it to a
-    mapping method (this is your second error).
-
-    To use the "mapCompute" routine, you must create a mapping strategy, not
-    a bipartitioning one, and so use stratGraphMap() and not
-    stratGraphBipart(). Your mapping strategy should however be based on the
-    "recursive bipartitioning" method ("b"). For instance, a simple (and
-    hence not very efficient) mapping strategy can be :
-
-    "b{sep=f}"
-
-    which computes mappings with the recursive bipartitioning method "b",
-    this latter using the Fiduccia-Mattheyses method "f" to compute its
-    separators.
-
-    If you want an exact partition (see your previous post), try
-    "b{sep=fx}".
-
-    However, these strategies are not the most efficient, as they do not
-    make use of the multi-level framework.
-
-    To use the multi-level framework, try for instance:
-
-    "b{sep=m{vert=100,low=h,asc=f}x}"
-
-    The current default mapping strategy in Scotch can be seen by using the
-    "-vs" option of program gmap. It is, to date:
-
-    b
-    {
-        job=t,
-        map=t,
-        poli=S,
-        sep=
-        (
-            m
-            {
-                asc=b
-                {
-                    bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
-                    org=f{move=80,pass=-1,bal=0.005},
-                    width=3
-                },
-                low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
-                type=h,
-                vert=80,
-                rat=0.8
-            }
-          | m
-            {
-                asc=b
-                {
-                    bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
-                    org=f{move=80,pass=-1,bal=0.005},
-                    width=3
-                },
-                low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
-                type=h,
-                vert=80,
-                rat=0.8
-            }
-        )
-    }
-
-
 \*---------------------------------------------------------------------------*/
 
 #include "scotchDecomp.H"
@@ -239,42 +158,4 @@ Foam::labelList Foam::scotchDecomp::decompose
 }
 
 
-void Foam::scotchDecomp::calcCSR
-(
-    const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    FatalErrorIn
-    (
-        "labelList scotchDecomp::decompose"
-        "("
-            "const polyMesh&, "
-            "const List<int>&, "
-            "const List<int>&"
-        ")"
-    )   << notImplementedMessage << exit(FatalError);
-}
-
-
-void Foam::scotchDecomp::calcCSR
-(
-    const labelListList& cellCells,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    FatalErrorIn
-    (
-        "labelList scotchDecomp::decompose"
-        "("
-            "const labelListList&, "
-            "const List<int>&, "
-            "const List<int>&"
-        ")"
-    )   << notImplementedMessage << exit(FatalError);
-}
-
-
 // ************************************************************************* //

src/parallel/decompose/Allwmake

@@ -9,7 +9,7 @@ wmake libso metisDecomp
 
 if [ -d "$FOAM_MPI_LIBBIN" ]
 then
-    ( WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB; wmake libso parMetisDecomp )
+    ( WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB; wmake libso ptscotchDecomp && wmake libso parMetisDecomp )
 fi
 
 wmake libso decompositionMethods

src/parallel/decompose/AllwmakeLnInclude

@@ -6,5 +6,6 @@ wmakeLnInclude decompositionMethods
 wmakeLnInclude metisDecomp
 wmakeLnInclude parMetisDecomp
 wmakeLnInclude scotchDecomp
+wmakeLnInclude ptscotchDecomp
 
 # ----------------------------------------------------------------- end-of-file

src/parallel/decompose/decompositionMethods/Make/options

@@ -1,9 +1,9 @@
 EXE_INC =
 
 LIB_LIBS = \
-    -L$(FOAM_LIBBIN)/dummy \
-    -L$(FOAM_MPI_LIBBIN) \
-    -lscotchDecomp \
-    -lmetisDecomp \
-    -lparMetisDecomp
+    /* -L$(FOAM_LIBBIN)/dummy */ \
+    /* -L$(FOAM_MPI_LIBBIN) */ \
+    /* -lscotchDecomp */ \
+    /* -lmetisDecomp */ \
+    /* -lparMetisDecomp */
 

src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C

@@ -28,6 +28,9 @@ InClass
 \*---------------------------------------------------------------------------*/
 
 #include "decompositionMethod.H"
+#include "globalIndex.H"
+#include "cyclicPolyPatch.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -156,6 +159,18 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cc
+)
+{
+    scalarField cWeights(0);
+
+    return decompose(globalCellCells, cc, cWeights);
+}
+
+
 void Foam::decompositionMethod::calcCellCells
 (
     const polyMesh& mesh,
@@ -201,15 +216,284 @@ void Foam::decompositionMethod::calcCellCells
 }
 
 
-Foam::labelList Foam::decompositionMethod::decompose
+void Foam::decompositionMethod::calcCSR
 (
-    const labelListList& globalCellCells,
-    const pointField& cc
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
 )
 {
-    scalarField cWeights(0);
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
 
-    return decompose(globalCellCells, cc, cWeights);
+    xadj.setSize(mesh.nCells()+1);
+
+    // Initialise the number of internal faces of the cells to twice the
+    // number of internal faces
+    label nInternalFaces = 2*mesh.nInternalFaces();
+
+    // Check the boundary for coupled patches and add to the number of
+    // internal faces
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            nInternalFaces += pbm[patchi].size();
+        }
+    }
+
+    // Create the adjncy array the size of the total number of internal and
+    // coupled faces
+    adjncy.setSize(nInternalFaces);
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        const labelList& cFaces = mesh.cells()[cellI];
+
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if
+            (
+                mesh.isInternalFace(faceI)
+             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
+            )
+            {
+                freeAdj++;
+            }
+        }
+    }
+    xadj[mesh.nCells()] = freeAdj;
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    labelList nFacesPerCell(mesh.nCells(), 0);
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = mesh.faceOwner()[faceI];
+        label nei = mesh.faceNeighbour()[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+    }
+
+    // Coupled faces. Only cyclics done.
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            const unallocLabelList& faceCells = pbm[patchi].faceCells();
+
+            label sizeby2 = faceCells.size()/2;
+
+            for (label facei=0; facei<sizeby2; facei++)
+            {
+                label own = faceCells[facei];
+                label nei = faceCells[facei + sizeby2];
+
+                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+            }
+        }
+    }
+}
+
+
+// From cell-cell connections to Metis format (like CompactListList)
+void Foam::decompositionMethod::calcCSR
+(
+    const labelListList& cellCells,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Count number of internal faces
+    label nConnections = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        nConnections += cellCells[coarseI].size();
+    }
+
+    // Create the adjncy array as twice the size of the total number of
+    // internal faces
+    adjncy.setSize(nConnections);
+
+    xadj.setSize(cellCells.size()+1);
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        xadj[coarseI] = freeAdj;
+
+        const labelList& cCells = cellCells[coarseI];
+
+        forAll(cCells, i)
+        {
+            adjncy[freeAdj++] = cCells[i];
+        }
+    }
+    xadj[cellCells.size()] = freeAdj;
+}
+
+
+void Foam::decompositionMethod::calcDistributedCSR
+(
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Create global cell numbers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalIndex globalCells(mesh.nCells());
+
+
+    //
+    // Make Metis Distributed CSR (Compressed Storage Format) storage
+    //   adjncy      : contains cellCells (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    //
+
+
+    const labelList& faceOwner = mesh.faceOwner();
+    const labelList& faceNeighbour = mesh.faceNeighbour();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+
+    // Get renumbered owner on other side of coupled faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start() - mesh.nInternalFaces();
+
+            forAll(pp, i)
+            {
+                globalNeighbour[bFaceI++] = globalCells.toGlobal
+                (
+                    faceOwner[faceI++]
+                );
+            }
+        }
+    }
+
+    // Get the cell on the other side of coupled patches
+    syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
+
+
+    // Count number of faces (internal + coupled)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Number of faces per cell
+    List<int> nFacesPerCell(mesh.nCells(), 0);
+
+    // Number of coupled faces
+    label nCoupledFaces = 0;
+
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        nFacesPerCell[faceOwner[faceI]]++;
+        nFacesPerCell[faceNeighbour[faceI]]++;
+    }
+    // Handle coupled faces
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                nCoupledFaces++;
+                nFacesPerCell[faceOwner[faceI++]]++;
+            }
+        }
+    }
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+
+    xadj.setSize(mesh.nCells()+1);
+
+    int freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        freeAdj += nFacesPerCell[cellI];
+    }
+    xadj[mesh.nCells()] = freeAdj;
+
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
+
+    nFacesPerCell = 0;
+
+    // For internal faces is just offsetted owner and neighbour
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = faceOwner[faceI];
+        label nei = faceNeighbour[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
+    }
+    // For boundary faces is offsetted coupled neighbour
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start()-mesh.nInternalFaces();
+
+            forAll(pp, i)
+            {
+                label own = faceOwner[faceI];
+                adjncy[xadj[own] + nFacesPerCell[own]++] =
+                    globalNeighbour[bFaceI];
+
+                faceI++;
+                bFaceI++;
+            }
+        }
+    }
 }
 
 

src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.H

@@ -66,6 +66,31 @@ protected:
             labelListList& cellCells
         );
 
+        // From mesh to compact row storage format
+        // (like CompactListList)
+        static void calcCSR
+        (
+            const polyMesh& mesh,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        // From cell-cell connections to compact row storage format
+        // (like CompactListList)
+        static void calcCSR
+        (
+            const labelListList& cellCells,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        static void calcDistributedCSR
+        (
+            const polyMesh& mesh,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
 private:
 
     // Private Member Functions

src/parallel/decompose/metisDecomp/Make/options

@@ -1,7 +1,6 @@
 EXE_INC = \
     -I$(WM_THIRD_PARTY_DIR)/metis-5.0pre2/include \
-    -I../decompositionMethods/lnInclude \
-    -I../scotchDecomp/lnInclude
+    -I../decompositionMethods/lnInclude
 
 LIB_LIBS = \
     -lmetis \

src/parallel/decompose/metisDecomp/metisDecomp.C

@@ -28,7 +28,6 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
 #include "Time.H"
-#include "scotchDecomp.H"
 
 extern "C"
 {
@@ -340,7 +339,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
     List<int> adjncy;
     List<int> xadj;
-    scotchDecomp::calcCSR(mesh_, adjncy, xadj);
+    calcCSR(mesh_, adjncy, xadj);
 
     // Decompose using default weights
     List<int> finalDecomp;
@@ -390,7 +389,7 @@ Foam::labelList Foam::metisDecomp::decompose
             cellCells
         );
 
-        scotchDecomp::calcCSR(cellCells, adjncy, xadj);
+        calcCSR(cellCells, adjncy, xadj);
     }
 
     // Decompose using default weights
@@ -435,7 +434,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
     List<int> adjncy;
     List<int> xadj;
-    scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
+    calcCSR(globalCellCells, adjncy, xadj);
 
 
     // Decompose using default weights

src/parallel/decompose/parMetisDecomp/parMetisDecomp.C

@@ -26,7 +26,6 @@ License
 
 #include "parMetisDecomp.H"
 #include "metisDecomp.H"
-#include "scotchDecomp.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
@@ -406,7 +405,7 @@ Foam::labelList Foam::parMetisDecomp::decompose
     Field<int> adjncy;
     // Offsets into adjncy
     Field<int> xadj;
-    calcMetisDistributedCSR
+    calcDistributedCSR
     (
         mesh_,
         adjncy,
@@ -774,7 +773,7 @@ Foam::labelList Foam::parMetisDecomp::decompose
     Field<int> adjncy;
     // Offsets into adjncy
     Field<int> xadj;
-    scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
+    calcCSR(globalCellCells, adjncy, xadj);
 
     // decomposition options. 0 = use defaults
     List<int> options(3, 0);
@@ -871,146 +870,4 @@ Foam::labelList Foam::parMetisDecomp::decompose
 }
 
 
-void Foam::parMetisDecomp::calcMetisDistributedCSR
-(
-    const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    // Create global cell numbers
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    globalIndex globalCells(mesh.nCells());
-
-
-    //
-    // Make Metis Distributed CSR (Compressed Storage Format) storage
-    //   adjncy      : contains cellCells (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-    //
-
-
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeighbour = mesh.faceNeighbour();
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-
-    // Get renumbered owner on other side of coupled faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            label bFaceI = pp.start() - mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                globalNeighbour[bFaceI++] = globalCells.toGlobal
-                (
-                    faceOwner[faceI++]
-                );
-            }
-        }
-    }
-
-    // Get the cell on the other side of coupled patches
-    syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
-
-
-    // Count number of faces (internal + coupled)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Number of faces per cell
-    List<int> nFacesPerCell(mesh.nCells(), 0);
-
-    // Number of coupled faces
-    label nCoupledFaces = 0;
-
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        nFacesPerCell[faceOwner[faceI]]++;
-        nFacesPerCell[faceNeighbour[faceI]]++;
-    }
-    // Handle coupled faces
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                nCoupledFaces++;
-                nFacesPerCell[faceOwner[faceI++]]++;
-            }
-        }
-    }
-
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-
-    xadj.setSize(mesh.nCells()+1);
-
-    int freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        freeAdj += nFacesPerCell[cellI];
-    }
-    xadj[mesh.nCells()] = freeAdj;
-
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
-
-    adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
-
-    nFacesPerCell = 0;
-
-    // For internal faces is just offsetted owner and neighbour
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label own = faceOwner[faceI];
-        label nei = faceNeighbour[faceI];
-
-        adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
-    }
-    // For boundary faces is offsetted coupled neighbour
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            label bFaceI = pp.start()-mesh.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                label own = faceOwner[faceI];
-                adjncy[xadj[own] + nFacesPerCell[own]++] =
-                    globalNeighbour[bFaceI];
-
-                faceI++;
-                bFaceI++;
-            }
-        }
-    }
-}
-
-
 // ************************************************************************* //

src/parallel/decompose/ptscotchDecomp/Make/files

@@ -0,0 +1,3 @@
+ptscotchDecomp.C
+
+LIB = $(FOAM_MPI_LIBBIN)/libptscotchDecomp

src/parallel/decompose/ptscotchDecomp/Make/options

@@ -0,0 +1,12 @@
+include $(RULES)/mplib$(WM_MPLIB)
+
+EXE_INC = \
+    $(PFLAGS) $(PINC) \
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
+    -I/usr/include/scotch \
+    /* -I$(LIB_SRC)/parallel/decompose/parMetisDecomp */ \
+    /* -I$(LIB_SRC)/parallel/decompose/scotchDecomp */ \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
+
+LIB_LIBS = \
+    -L$(FOAM_MPI_LIBBIN) -lptscotch -lptscotcherrexit

src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C

@@ -0,0 +1,592 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+    From scotch forum:
+ 	
+    By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ]  
+    2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order
+    not to be confused, you must have a clear view of how they are built.
+    Here are some rules:
+
+    1- Strategies are made up of "methods" which are combined by means of
+    "operators".
+
+    2- A method is of the form "m{param=value,param=value,...}", where "m"
+    is a single character (this is your first error: "f" is a method name,
+    not a parameter name).
+
+    3- There exist different sort of strategies : bipartitioning strategies,
+    mapping strategies, ordering strategies, which cannot be mixed. For
+    instance, you cannot build a bipartitioning strategy and feed it to a
+    mapping method (this is your second error).
+
+    To use the "mapCompute" routine, you must create a mapping strategy, not
+    a bipartitioning one, and so use stratGraphMap() and not
+    stratGraphBipart(). Your mapping strategy should however be based on the
+    "recursive bipartitioning" method ("b"). For instance, a simple (and
+    hence not very efficient) mapping strategy can be :
+
+    "b{sep=f}"
+
+    which computes mappings with the recursive bipartitioning method "b",
+    this latter using the Fiduccia-Mattheyses method "f" to compute its
+    separators.
+
+    If you want an exact partition (see your previous post), try
+    "b{sep=fx}".
+
+    However, these strategies are not the most efficient, as they do not
+    make use of the multi-level framework.
+
+    To use the multi-level framework, try for instance:
+
+    "b{sep=m{vert=100,low=h,asc=f}x}"
+
+    The current default mapping strategy in Scotch can be seen by using the
+    "-vs" option of program gmap. It is, to date:
+
+    b
+    {
+        job=t,
+        map=t,
+        poli=S,
+        sep=
+        (
+            m
+            {
+                asc=b
+                {
+                    bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
+                    org=f{move=80,pass=-1,bal=0.005},
+                    width=3
+                },
+                low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
+                type=h,
+                vert=80,
+                rat=0.8
+            }
+          | m
+            {
+                asc=b
+                {
+                    bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
+                    org=f{move=80,pass=-1,bal=0.005},
+                    width=3
+                },
+                low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
+                type=h,
+                vert=80,
+                rat=0.8
+            }
+        )
+    }
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "ptscotchDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "OFstream.H"
+
+extern "C"
+{
+#include <stdio.h>
+#include "mpi.h"
+#include "ptscotch.h"
+}
+
+
+// Hack: scotch generates floating point errors so need to switch of error
+//       trapping!
+#if defined(linux) || defined(linuxAMD64) || defined(linuxIA64)
+#    define LINUX
+#endif
+
+#if defined(LINUX) && defined(__GNUC__)
+#    define LINUX_GNUC
+#endif
+
+#ifdef LINUX_GNUC
+#   ifndef __USE_GNU
+#       define __USE_GNU
+#   endif
+#   include <fenv.h>
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(ptscotchDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        ptscotchDecomp,
+        dictionaryMesh
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::ptscotchDecomp::check(const int retVal, const char* str)
+{
+    if (retVal)
+    {
+        FatalErrorIn("ptscotchDecomp::decompose(..)")
+            << "Call to scotch routine " << str << " failed."
+            << exit(FatalError);
+    }
+}
+
+
+// Call scotch with options from dictionary.
+Foam::label Foam::ptscotchDecomp::decompose
+(
+    List<int>& adjncy,
+    List<int>& xadj,
+    const scalarField& cWeights,
+
+    List<int>& finalDecomp
+)
+{
+//    // Dump graph
+//    if (decompositionDict_.found("ptscotchCoeffs"))
+//    {
+//        const dictionary& scotchCoeffs =
+//            decompositionDict_.subDict("ptscotchCoeffs");
+//
+//        if (scotchCoeffs.found("writeGraph"))
+//        {
+//            Switch writeGraph(scotchCoeffs.lookup("writeGraph"));
+//
+//            if (writeGraph)
+//            {
+//                OFstream str(mesh_.time().path() / mesh_.name() + ".grf");
+//
+//                Info<< "Dumping Scotch graph file to " << str.name() << endl
+//                    << "Use this in combination with gpart." << endl;
+//
+//                label version = 0;
+//                str << version << nl;
+//                // Numer of vertices
+//                str << xadj.size()-1 << ' ' << adjncy.size() << nl;
+//                // Numbering starts from 0
+//                label baseval = 0;
+//                // Has weights?
+//                label hasEdgeWeights = 0;
+//                label hasVertexWeights = 0;
+//                label numericflag = 10*hasEdgeWeights+hasVertexWeights;
+//                str << baseval << ' ' << numericflag << nl;
+//                for (label cellI = 0; cellI < xadj.size()-1; cellI++)
+//                {
+//                    label start = xadj[cellI];
+//                    label end = xadj[cellI+1];
+//                    str << end-start;
+//
+//                    for (label i = start; i < end; i++)
+//                    {
+//                        str << ' ' << adjncy[i];
+//                    }
+//                    str << nl;
+//                }
+//            }
+//        }
+//    }
+
+
+    // Strategy
+    // ~~~~~~~~
+
+    // Default.
+    SCOTCH_Strat stradat;
+    check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
+
+    if (decompositionDict_.found("scotchCoeffs"))
+    {
+        const dictionary& scotchCoeffs =
+            decompositionDict_.subDict("scotchCoeffs");
+
+
+        string strategy;
+        if (scotchCoeffs.readIfPresent("strategy", strategy))
+        {
+            if (debug)
+            {
+                Info<< "ptscotchDecomp : Using strategy " << strategy << endl;
+            }
+            SCOTCH_stratDgraphMap(&stradat, strategy.c_str());
+            //fprintf(stdout, "S\tStrat=");
+            //SCOTCH_stratSave(&stradat, stdout);
+            //fprintf(stdout, "\n");
+        }
+    }
+
+
+    // Graph
+    // ~~~~~
+
+    List<int> velotab;
+
+
+    // Check for externally provided cellweights and if so initialise weights
+    scalar minWeights = gMin(cWeights);
+    if (cWeights.size() > 0)
+    {
+        if (minWeights <= 0)
+        {
+            WarningIn
+            (
+                "ptscotchDecomp::decompose"
+                "(const pointField&, const scalarField&)"
+            )   << "Illegal minimum weight " << minWeights
+                << endl;
+        }
+
+        if (cWeights.size() != xadj.size()-1)
+        {
+            FatalErrorIn
+            (
+                "ptscotchDecomp::decompose"
+                "(const pointField&, const scalarField&)"
+            )   << "Number of cell weights " << cWeights.size()
+                << " does not equal number of cells " << xadj.size()-1
+                << exit(FatalError);
+        }
+
+        // Convert to integers.
+        velotab.setSize(cWeights.size());
+        forAll(velotab, i)
+        {
+            velotab[i] = int(cWeights[i]/minWeights);
+        }
+    }
+
+
+
+    SCOTCH_Dgraph grafdat;
+    check(SCOTCH_dgraphInit(&grafdat, MPI_COMM_WORLD), "SCOTCH_dgraphInit");
+    check
+    (
+        SCOTCH_dgraphBuild
+        (
+            &grafdat,               // grafdat
+            0,                      // baseval, c-style numbering
+            xadj.size()-1,          // vertlocnbr, nCells
+            xadj.size()-1,          // vertlocmax
+            const_cast<SCOTCH_Num*>(xadj.begin()),  // vertloctab, start index per cell into
+                                    // adjncy
+            &xadj[1],               // vendloctab, end index  ,,
+
+            velotab.begin(),        // veloloctab, vertex weights
+            NULL,                   // vlblloctab
+
+            adjncy.size(),          // edgelocnbr, number of arcs
+            adjncy.size(),          // edgelocsiz
+            adjncy.begin(),         // edgeloctab
+            NULL,                   // edgegsttab
+            NULL                    // edlotab, edge weights
+        ),
+        "SCOTCH_dgraphBuild"
+    );
+    check(SCOTCH_dgraphCheck(&grafdat), "SCOTCH_dgraphCheck");
+
+
+    // Architecture
+    // ~~~~~~~~~~~~
+    // (fully connected network topology since using switch)
+
+    SCOTCH_Arch archdat;
+    check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
+
+    List<label> processorWeights;
+    if (decompositionDict_.found("scotchCoeffs"))
+    {
+        const dictionary& scotchCoeffs =
+            decompositionDict_.subDict("scotchCoeffs");
+
+        scotchCoeffs.readIfPresent("processorWeights", processorWeights);
+    }
+    if (processorWeights.size())
+    {
+        if (debug)
+        {
+            Info<< "ptscotchDecomp : Using procesor weights " << processorWeights
+                << endl;
+        }
+        check
+        (
+            SCOTCH_archCmpltw(&archdat, nProcessors_, processorWeights.begin()),
+            "SCOTCH_archCmpltw"
+        );
+    }
+    else
+    {
+        check
+        (
+            SCOTCH_archCmplt(&archdat, nProcessors_),
+            "SCOTCH_archCmplt"
+        );
+    }
+
+
+    //SCOTCH_Mapping mapdat;
+    //SCOTCH_dgraphMapInit(&grafdat, &mapdat, &archdat, NULL);
+    //SCOTCH_dgraphMapCompute(&grafdat, &mapdat, &stradat); /*Perform mapping*/
+    //SCOTCHdgraphMapExit(&grafdat, &mapdat);
+
+
+    // Hack:switch off fpu error trapping
+#   ifdef LINUX_GNUC
+    int oldExcepts = fedisableexcept
+    (
+        FE_DIVBYZERO
+      | FE_INVALID
+      | FE_OVERFLOW
+    );
+#   endif
+
+    finalDecomp.setSize(xadj.size()-1);
+    finalDecomp = 0;
+    check
+    (
+        SCOTCH_dgraphMap
+        (
+            &grafdat,
+            &archdat,
+            &stradat,           // const SCOTCH_Strat *
+            finalDecomp.begin() // parttab
+        ),
+        "SCOTCH_graphMap"
+    );
+
+#   ifdef LINUX_GNUC
+    feenableexcept(oldExcepts);
+#   endif
+
+
+
+    //finalDecomp.setSize(xadj.size()-1);
+    //check
+    //(
+    //    SCOTCH_dgraphPart
+    //    (
+    //        &grafdat,
+    //        nProcessors_,       // partnbr
+    //        &stradat,           // const SCOTCH_Strat *
+    //        finalDecomp.begin() // parttab
+    //    ),
+    //    "SCOTCH_graphPart"
+    //);
+
+    // Release storage for graph
+    SCOTCH_dgraphExit(&grafdat);
+    // Release storage for strategy
+    SCOTCH_stratExit(&stradat);
+    // Release storage for network topology
+    SCOTCH_archExit(&archdat);
+
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ptscotchDecomp::ptscotchDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    decompositionMethod(decompositionDict),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
+{
+    if (points.size() != mesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "ptscotchDecomp::decompose(const pointField&, const scalarField&)"
+        )
+            << "Can use this decomposition method only for the whole mesh"
+            << endl
+            << "and supply one coordinate (cellCentre) for every cell." << endl
+            << "The number of coordinates " << points.size() << endl
+            << "The number of cells in the mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+//    // For running sequential ...
+//    if (Pstream::nProcs() <= 1)
+//    {
+//        return scotchDecomp(decompositionDict_, mesh_)
+//            .decompose(points, pointWeights);
+//    }
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    // Connections
+    Field<int> adjncy;
+    // Offsets into adjncy
+    Field<int> xadj;
+    calcDistributedCSR
+    (
+        mesh_,
+        adjncy,
+        xadj
+    );
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, pointWeights, finalDecomp);
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const labelList& agglom,
+    const pointField& agglomPoints,
+    const scalarField& pointWeights
+)
+{
+    if (agglom.size() != mesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "ptscotchDecomp::decompose(const labelList&, const pointField&)"
+        )   << "Size of cell-to-coarse map " << agglom.size()
+            << " differs from number of cells in mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+//    // For running sequential ...
+//    if (Pstream::nProcs() <= 1)
+//    {
+//        return scotchDecomp(decompositionDict_, mesh_)
+//            .decompose(agglom, agglomPoints, pointWeights);
+//    }
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    List<int> adjncy;
+    List<int> xadj;
+    {
+        // Get cellCells on coarse mesh.
+        labelListList cellCells;
+        calcCellCells
+        (
+            mesh_,
+            agglom,
+            agglomPoints.size(),
+            cellCells
+        );
+
+        calcCSR(cellCells, adjncy, xadj);
+    }
+
+    // Decompose using weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, pointWeights, finalDecomp);
+
+    // Rework back into decomposition for original mesh_
+    labelList fineDistribution(agglom.size());
+
+    forAll(fineDistribution, i)
+    {
+        fineDistribution[i] = finalDecomp[agglom[i]];
+    }
+
+    return fineDistribution;
+}
+
+
+Foam::labelList Foam::ptscotchDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cellCentres,
+    const scalarField& cWeights
+)
+{
+    if (cellCentres.size() != globalCellCells.size())
+    {
+        FatalErrorIn
+        (
+            "ptscotchDecomp::decompose(const pointField&, const labelListList&)"
+        )   << "Inconsistent number of cells (" << globalCellCells.size()
+            << ") and number of cell centres (" << cellCentres.size()
+            << ")." << exit(FatalError);
+    }
+
+//    // For running sequential ...
+//    if (Pstream::nProcs() <= 1)
+//    {
+//        return scotchDecomp(decompositionDict_, mesh_)
+//            .decompose(globalCellCells, cellCentres, cWeights);
+//    }
+
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    List<int> adjncy;
+    List<int> xadj;
+    calcCSR(globalCellCells, adjncy, xadj);
+
+    // Decompose using weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, cWeights, finalDecomp);
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
+// ************************************************************************* //

src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::ptscotchDecomp
+
+Description
+    PTScotch domain decomposition
+
+SourceFiles
+    ptscotchDecomp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ptscotchDecomp_H
+#define ptscotchDecomp_H
+
+#include "decompositionMethod.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class ptscotchDecomp Declaration
+\*---------------------------------------------------------------------------*/
+
+class ptscotchDecomp
+:
+    public decompositionMethod
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+
+    // Private Member Functions
+
+        //- Check and print error message
+        static void check(const int, const char*);
+
+        label decompose
+        (
+            List<int>& adjncy,
+            List<int>& xadj,
+            const scalarField& cWeights,
+            List<int>& finalDecomp
+        );
+
+        //- Disallow default bitwise copy construct and assignment
+        void operator=(const ptscotchDecomp&);
+        ptscotchDecomp(const ptscotchDecomp&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ptscotch");
+
+
+    // Constructors
+
+        //- Construct given the decomposition dictionary and mesh
+        ptscotchDecomp
+        (
+            const dictionary& decompositionDict,
+            const polyMesh& mesh
+        );
+
+
+    // Destructor
+
+        virtual ~ptscotchDecomp()
+        {}
+
+
+    // Member Functions
+
+        virtual bool parallelAware() const
+        {
+            // ptscotch does not know about proc boundaries
+            return true;
+        }
+
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
+
+        //- Return for every coordinate the wanted processor number. Gets
+        //  passed agglomeration map (from fine to coarse cells) and coarse cell
+        //  location. Can be overridden by decomposers that provide this
+        //  functionality natively.
+        virtual labelList decompose
+        (
+            const labelList& agglom,
+            const pointField& regionPoints,
+            const scalarField& regionWeights
+        );
+
+        //- Return for every coordinate the wanted processor number. Explicitly
+        //  provided mesh connectivity.
+        //  The connectivity is equal to mesh.cellCells() except for
+        //  - in parallel the cell numbers are global cell numbers (starting
+        //    from 0 at processor0 and then incrementing all through the
+        //    processors)
+        //  - the connections are across coupled patches
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights
+        );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/parallel/decompose/scotchDecomp/scotchDecomp.C

@@ -556,144 +556,4 @@ Foam::labelList Foam::scotchDecomp::decompose
 }
 
 
-void Foam::scotchDecomp::calcCSR
-(
-    const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    // Make Metis CSR (Compressed Storage Format) storage
-    //   adjncy      : contains neighbours (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-
-    xadj.setSize(mesh.nCells()+1);
-
-    // Initialise the number of internal faces of the cells to twice the
-    // number of internal faces
-    label nInternalFaces = 2*mesh.nInternalFaces();
-
-    // Check the boundary for coupled patches and add to the number of
-    // internal faces
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            nInternalFaces += pbm[patchi].size();
-        }
-    }
-
-    // Create the adjncy array the size of the total number of internal and
-    // coupled faces
-    adjncy.setSize(nInternalFaces);
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        const labelList& cFaces = mesh.cells()[cellI];
-
-        forAll(cFaces, i)
-        {
-            label faceI = cFaces[i];
-
-            if
-            (
-                mesh.isInternalFace(faceI)
-             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
-            )
-            {
-                freeAdj++;
-            }
-        }
-    }
-    xadj[mesh.nCells()] = freeAdj;
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
-
-    labelList nFacesPerCell(mesh.nCells(), 0);
-
-    // Internal faces
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label own = mesh.faceOwner()[faceI];
-        label nei = mesh.faceNeighbour()[faceI];
-
-        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-    }
-
-    // Coupled faces. Only cyclics done.
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            const unallocLabelList& faceCells = pbm[patchi].faceCells();
-
-            label sizeby2 = faceCells.size()/2;
-
-            for (label facei=0; facei<sizeby2; facei++)
-            {
-                label own = faceCells[facei];
-                label nei = faceCells[facei + sizeby2];
-
-                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-            }
-        }
-    }
-}
-
-
-// From cell-cell connections to Metis format (like CompactListList)
-void Foam::scotchDecomp::calcCSR
-(
-    const labelListList& cellCells,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    // Count number of internal faces
-    label nConnections = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        nConnections += cellCells[coarseI].size();
-    }
-
-    // Create the adjncy array as twice the size of the total number of
-    // internal faces
-    adjncy.setSize(nConnections);
-
-    xadj.setSize(cellCells.size()+1);
-
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        xadj[coarseI] = freeAdj;
-
-        const labelList& cCells = cellCells[coarseI];
-
-        forAll(cCells, i)
-        {
-            adjncy[freeAdj++] = cCells[i];
-        }
-    }
-    xadj[cellCells.size()] = freeAdj;
-}
-
-
-
 // ************************************************************************* //

src/parallel/decompose/scotchDecomp/scotchDecomp.H

@@ -139,26 +139,6 @@ public:
             const pointField& cc,
             const scalarField& cWeights
         );
-
-
-        //- Helper to convert local connectivity (supplied as owner,neighbour)
-        //  into CSR (Metis,scotch) storage. Does cyclics but not processor
-        //  patches
-        static void calcCSR
-        (
-            const polyMesh& mesh,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
-
-        //- Helper to convert connectivity (supplied as cellcells) into
-        //  CSR (Metis,scotch) storage
-        static void calcCSR
-        (
-            const labelListList& globalCellCells,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
 };