From a0597834a9199f076e032e83226d8807f5e02eac Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 21 Feb 2012 13:00:35 +0000
Subject: [PATCH] ENH: renumberMesh: added Sloan renumbering

---
 .../mesh/manipulation/renumberMesh/Allwmake   |  14 +
 .../manipulation/renumberMesh/Make/options    |   2 +
 .../manipulation/renumberMesh/renumberMesh.C  |  12 +-
 .../renumberMesh/renumberMeshDict             |   2 +-
 src/renumber/Allwmake                         |  10 +
 src/renumber/SloanRenumber/Make/files         |   3 +
 src/renumber/SloanRenumber/Make/options       |  12 +
 src/renumber/SloanRenumber/SloanRenumber.C    | 263 ++++++++++++++++++
 src/renumber/SloanRenumber/SloanRenumber.H    | 123 ++++++++
 9 files changed, 430 insertions(+), 11 deletions(-)
 create mode 100755 applications/utilities/mesh/manipulation/renumberMesh/Allwmake
 create mode 100644 src/renumber/SloanRenumber/Make/files
 create mode 100644 src/renumber/SloanRenumber/Make/options
 create mode 100644 src/renumber/SloanRenumber/SloanRenumber.C
 create mode 100644 src/renumber/SloanRenumber/SloanRenumber.H

diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Allwmake b/applications/utilities/mesh/manipulation/renumberMesh/Allwmake
new file mode 100755
index 0000000000..f5f12d5993
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/renumberMesh/Allwmake
@@ -0,0 +1,14 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+export SLOAN_LINK_FLAGS=''
+
+if [ -f "${FOAM_LIBBIN}/libSloanRenumber.so" ]
+then
+    echo "Found libSloanRenumber.so  --  enabling Sloan renumbering support."
+    export SLOAN_LINK_FLAGS="-lSloanRenumber"
+fi
+
+wmake
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Make/options b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
index 6a2c95738d..72eeafc9d1 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
@@ -1,4 +1,5 @@
 EXE_INC = \
+    -DFULLDEBUG -g -O0 \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -11,4 +12,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lgenericPatchFields \
     -lrenumberMethods \
+    $(SLOAN_LINK_FLAGS) \
     -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 487fde704a..4120f8c67a 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -131,9 +131,9 @@ void getBand
 
     forAll(nIntersect, cellI)
     {
-        for (label rowI = cellI-cellBandwidth[cellI]; rowI < cellI; rowI++)
+        for (label colI = cellI-cellBandwidth[cellI]; colI <= cellI; colI++)
         {
-            nIntersect[rowI]++;
+            nIntersect[colI]++;
         }
     }
 
@@ -598,14 +598,6 @@ int main(int argc, char *argv[])
         sumSqrIntersect
     );
 
-    if (band != getBand(mesh.faceOwner(), mesh.faceNeighbour()))
-    {
-        FatalErrorIn(args.executable())
-            << "band:" << band
-            << abort(FatalError);
-    }
-
-
     reduce(band, maxOp<label>());
     reduce(profile, sumOp<scalar>());
     scalar rmsFrontwidth = Foam::sqrt
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
index bc7084616d..8878469010 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
@@ -37,6 +37,7 @@ sortCoupledFaceCells false;
 
 
 method          CuthillMcKee;
+//method          Sloan;
 //method          manual;
 //method          random;
 //method          spring;
@@ -47,7 +48,6 @@ method          CuthillMcKee;
 //    reverse true;
 //}
 
-
 manualCoeffs
 {
     // In system directory: new-to-original (i.e. order) labelIOList
diff --git a/src/renumber/Allwmake b/src/renumber/Allwmake
index 36c057dafe..3e28994c25 100755
--- a/src/renumber/Allwmake
+++ b/src/renumber/Allwmake
@@ -18,6 +18,16 @@ set -x
 
 wmake $makeType renumberMethods
 
+if [ -n "$BOOST_ARCH_PATH" ]
+then
+    wmake $makeType SloanRenumber
+else
+    echo
+    echo "Skipping SloanRenumber"
+    echo
+fi
+
+
 #if [ -n "$ZOLTAN_ARCH_PATH" ]
 #then
 #    wmake $makeType zoltanRenumber
diff --git a/src/renumber/SloanRenumber/Make/files b/src/renumber/SloanRenumber/Make/files
new file mode 100644
index 0000000000..f627ab1cb3
--- /dev/null
+++ b/src/renumber/SloanRenumber/Make/files
@@ -0,0 +1,3 @@
+SloanRenumber.C
+
+LIB = $(FOAM_LIBBIN)/libSloanRenumber
diff --git a/src/renumber/SloanRenumber/Make/options b/src/renumber/SloanRenumber/Make/options
new file mode 100644
index 0000000000..431c883479
--- /dev/null
+++ b/src/renumber/SloanRenumber/Make/options
@@ -0,0 +1,12 @@
+EXE_INC = \
+    /* -DFULLDEBUG -g -O0 */ \
+    -I$BOOST_ARCH_PATH/include \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/renumber/renumberMethods/lnInclude
+
+LIB_LIBS = \
+    -L$(BOOST_ARCH_PATH)/lib -lboost_thread \
+    -lmeshTools \
+    -ldecompositionMethods \
+    -lrenumberMethods
diff --git a/src/renumber/SloanRenumber/SloanRenumber.C b/src/renumber/SloanRenumber/SloanRenumber.C
new file mode 100644
index 0000000000..2490478940
--- /dev/null
+++ b/src/renumber/SloanRenumber/SloanRenumber.C
@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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/>.
+
+    Adapted from Boost graph/example/sloan_ordering.cpp
+
+\*---------------------------------------------------------------------------*/
+
+#include "SloanRenumber.H"
+#include "addToRunTimeSelectionTable.H"
+#include "decompositionMethod.H"
+#include "processorPolyPatch.H"
+#include "syncTools.H"
+
+#include <boost/config.hpp>
+#include <vector>
+#include <iostream>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/sloan_ordering.hpp>
+#include <boost/graph/properties.hpp>
+#include <boost/graph/bandwidth.hpp>
+#include <boost/graph/profile.hpp>
+#include <boost/graph/wavefront.hpp>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace boost;
+using namespace std;
+
+//Defining the graph type
+typedef adjacency_list
+<
+    setS,
+    vecS,
+    undirectedS,
+    property
+    <
+        vertex_color_t,
+        default_color_type,
+        property
+        <
+            vertex_degree_t,
+            Foam::label,
+            property
+            <
+                vertex_priority_t,
+                Foam::scalar
+            >
+        >
+    >
+> Graph;
+
+typedef graph_traits<Graph>::vertex_descriptor Vertex;
+typedef graph_traits<Graph>::vertices_size_type size_type;
+
+namespace Foam
+{
+    defineTypeNameAndDebug(SloanRenumber, 0);
+
+    addToRunTimeSelectionTable
+    (
+        renumberMethod,
+        SloanRenumber,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::SloanRenumber::SloanRenumber(const dictionary& renumberDict)
+:
+    renumberMethod(renumberDict),
+    reverse_
+    (
+        renumberDict.found(typeName + "Coeffs")
+      ? Switch(renumberDict.subDict(typeName + "Coeffs").lookup("reverse"))
+      : Switch(false)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::SloanRenumber::renumber
+(
+    const polyMesh& mesh,
+    const pointField& points
+) const
+{
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    // Construct graph : faceOwner + connections across cyclics.
+
+    // Determine neighbour cell
+    labelList nbr(mesh.nFaces()-mesh.nInternalFaces(), -1);
+    forAll(pbm, patchI)
+    {
+        if (pbm[patchI].coupled() && !isA<processorPolyPatch>(pbm[patchI]))
+        {
+            SubList<label>
+            (
+                nbr,
+                pbm[patchI].size(),
+                pbm[patchI].start()-mesh.nInternalFaces()
+            ).assign(pbm[patchI].faceCells());
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh, nbr);
+
+
+    Graph G(mesh.nCells());
+
+    // Add internal faces
+    forAll(mesh.faceNeighbour(), faceI)
+    {
+        add_edge(mesh.faceOwner()[faceI], mesh.faceNeighbour()[faceI], G);
+    }
+    // Add cyclics
+    forAll(pbm, patchI)
+    {
+        if
+        (
+            pbm[patchI].coupled()
+        && !isA<processorPolyPatch>(pbm[patchI])
+        &&  refCast<const coupledPolyPatch>(pbm[patchI]).owner()
+        )
+        {
+            const labelUList& faceCells = pbm[patchI].faceCells();
+            forAll(faceCells, i)
+            {
+                label bFaceI = pbm[patchI].start()+i-mesh.nInternalFaces();
+                label nbrCellI = nbr[bFaceI];
+
+                if (faceCells[i] < nbrCellI)
+                {
+                    add_edge(faceCells[i], nbrCellI, G);
+                }
+                else
+                {
+                    add_edge(nbrCellI, faceCells[i], G);
+                }
+            }
+        }
+    }
+
+
+    //Creating two iterators over the vertices
+    graph_traits<Graph>::vertex_iterator ui, ui_end;
+
+    //Creating a property_map with the degrees of the degrees of each vertex
+    property_map<Graph,vertex_degree_t>::type deg = get(vertex_degree, G);
+    for (boost::tie(ui, ui_end) = vertices(G); ui != ui_end; ++ui)
+        deg[*ui] = degree(*ui, G);
+
+    //Creating a property_map for the indices of a vertex
+    property_map<Graph, vertex_index_t>::type index_map = get(vertex_index, G);
+
+    //Creating a vector of vertices
+    std::vector<Vertex> sloan_order(num_vertices(G));
+
+    sloan_ordering
+    (
+        G,
+        sloan_order.begin(),
+        get(vertex_color, G),
+        make_degree_map(G),
+        get(vertex_priority, G)
+    );
+
+    labelList orderedToOld(sloan_order.size());
+    forAll(orderedToOld, c)
+    {
+        orderedToOld[c] = index_map[sloan_order[c]];
+    }
+
+    if (reverse_)
+    {
+        reverse(orderedToOld);
+    }
+
+    return orderedToOld;
+}
+
+
+Foam::labelList Foam::SloanRenumber::renumber
+(
+    const labelListList& cellCells,
+    const pointField& points
+) const
+{
+    Graph G(cellCells.size());
+
+    forAll(cellCells, cellI)
+    {
+        const labelList& nbrs = cellCells[cellI];
+        forAll(nbrs, i)
+        {
+            if (nbrs[i] > cellI)
+            {
+                add_edge(cellI, nbrs[i], G);
+            }
+        }
+    }
+
+    //Creating two iterators over the vertices
+    graph_traits<Graph>::vertex_iterator ui, ui_end;
+
+    //Creating a property_map with the degrees of the degrees of each vertex
+    property_map<Graph,vertex_degree_t>::type deg = get(vertex_degree, G);
+    for (boost::tie(ui, ui_end) = vertices(G); ui != ui_end; ++ui)
+        deg[*ui] = degree(*ui, G);
+
+    //Creating a property_map for the indices of a vertex
+    property_map<Graph, vertex_index_t>::type index_map = get(vertex_index, G);
+
+    //Creating a vector of vertices
+    std::vector<Vertex> sloan_order(num_vertices(G));
+
+    sloan_ordering
+    (
+        G,
+        sloan_order.begin(),
+        get(vertex_color, G),
+        make_degree_map(G),
+        get(vertex_priority, G)
+    );
+
+    labelList orderedToOld(sloan_order.size());
+    forAll(orderedToOld, c)
+    {
+        orderedToOld[c] = index_map[sloan_order[c]];
+    }
+
+    if (reverse_)
+    {
+        reverse(orderedToOld);
+    }
+
+    return orderedToOld;
+}
+
+
+// ************************************************************************* //
diff --git a/src/renumber/SloanRenumber/SloanRenumber.H b/src/renumber/SloanRenumber/SloanRenumber.H
new file mode 100644
index 0000000000..ddd42f70fc
--- /dev/null
+++ b/src/renumber/SloanRenumber/SloanRenumber.H
@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     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::SloanRenumber
+
+Description
+    Sloan renumbering algorithm
+
+    E.g. http://www.apav.it/sito_ratio/file_pdf/ratio_4/capitolo_1.pdf
+
+SourceFiles
+    SloanRenumber.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SloanRenumber_H
+#define SloanRenumber_H
+
+#include "renumberMethod.H"
+#include "Switch.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class SloanRenumber Declaration
+\*---------------------------------------------------------------------------*/
+
+class SloanRenumber
+:
+    public renumberMethod
+{
+    // Private data
+
+        const Switch reverse_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        void operator=(const SloanRenumber&);
+        SloanRenumber(const SloanRenumber&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Sloan");
+
+
+    // Constructors
+
+        //- Construct given the renumber dictionary
+        SloanRenumber(const dictionary& renumberDict);
+
+
+    //- Destructor
+    virtual ~SloanRenumber()
+    {}
+
+
+    // Member Functions
+
+        //- Return the order in which cells need to be visited, i.e.
+        //  from ordered back to original cell label.
+        //  This is only defined for geometric renumberMethods.
+        virtual labelList renumber(const pointField&) const
+        {
+            notImplemented("SloanRenumber::renumber(const pointField&)");
+            return labelList(0);
+        }
+
+        //- Return the order in which cells need to be visited, i.e.
+        //  from ordered back to original cell label.
+        //  Use the mesh connectivity (if needed)
+        virtual labelList renumber
+        (
+            const polyMesh& mesh,
+            const pointField& cc
+        ) const;
+
+        //- Return the order in which cells need to be visited, i.e.
+        //  from ordered back to original cell label.
+        //  The connectivity is equal to mesh.cellCells() except
+        //  - the connections are across coupled patches
+        virtual labelList renumber
+        (
+            const labelListList& cellCells,
+            const pointField& cc
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //