Parallel topological changes and internal combustion engine simulations

applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C

@@ -36,7 +36,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-
 #   include "setRootCase.H"
 
 #   include "createTime.H"

applications/solvers/engine/sonicTurbDyMEngineFoam/Make/options

@@ -4,10 +4,12 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/engine/lnInclude
 
 EXE_LIBS = \
     -lengine \
+    -lfiniteVolume \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -ltopoChangerFvMesh \
@@ -17,6 +19,5 @@ EXE_LIBS = \
     -lbasicThermophysicalModels \
     -lspecie \
     -lmeshTools \
-    -lfiniteVolume \
     $(WM_DECOMP_LIBS) \
     -llduSolvers

applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H

@@ -6,6 +6,13 @@
       + turbulence->divDevRhoReff(U)
     );
 
-    UEqn.relax();
+    if (oCorr == nOuterCorr - 1)
+    {
+        UEqn.relax(1);
+    }
+    else
+    {
+        UEqn.relax();
+    }
 
     solve(UEqn == -fvc::grad(p));

applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H

@@ -8,6 +8,8 @@
     );
     basicPsiThermo& thermo = pThermo();
 
+    // Make density field with zero gradient boundary conditions to handle
+    // attach-detach cases.  HJ, 20/Mar/2011
     volScalarField rho
     (
         IOobject
@@ -18,7 +20,8 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo.rho()
+        thermo.rho(),
+        zeroGradientFvPatchScalarField::typeName
     );
     rho.oldTime();
 
@@ -73,6 +76,7 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
+    volScalarField dpdt = fvc::ddt(p);
 
     volScalarField rUA
     (

applications/solvers/engine/sonicTurbDyMEngineFoam/infoDataOutput.H

@@ -1,17 +1,17 @@
 Pav << runTime.theta() << "   " << p.weightedAverage(mesh.V()).value() << endl;
 Tav << runTime.theta() << "   " << T.weightedAverage(mesh.V()).value() << endl;
 
-Info << "Max T = " << max(T) << " K" << endl;
+Info << "Max T = " << max(T).value() << " K" << endl;
-Info << "Min T = " << min(T) << " K" << endl;
+Info << "Min T = " << min(T).value() << " K" << endl;
-Info << "Max p = " << max(p)/1.0e5 << " bar" << endl;
+Info << "Max p = " << max(p).value()/1.0e5 << " bar" << endl;
 
 massBal << runTime.theta() << "   "
         << fvc::domainIntegrate(rho).value() << endl;
 volume << runTime.theta() << "   " << sum(mesh.V()) << endl;
 
 debugT << nl << "Crank angle: " << runTime.theta() << endl;
-debugT << "Max T = " << max(T) << " K" << endl;
+debugT << "Max T = " << max(T).value() << " K" << endl;
-debugT << "Min T = " << min(T) << " K" << endl;
+debugT << "Min T = " << min(T).value() << " K" << endl;
 
 kav << runTime.theta() << "   "
     << (turbulence->k())().weightedAverage(mesh.V()).value() << endl;

applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H

@@ -2,45 +2,71 @@
     rho = thermo.rho();
 
     rUA = 1.0/UEqn.A();
-    H = UEqn.H();
     U = rUA*UEqn.H();
 
-    phi = fvc::interpolate(rho)
+    if (nOuterCorr != 1)
-         *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
-
-    // Store pressure for under-relaxation
-    p.storePrevIter();
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        p.storePrevIter();
-        (
+    }
-            fvm::ddt(psi, p)
-          + fvc::div(phi)
-          - fvm::laplacian(rho*rUA, p)
-        );
 
-        pEqn.solve();
+    if (transonic)
+    {
+        surfaceScalarField phid =
+            fvc::interpolate(thermo.psi())*
+            ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-        if (nonOrth == nNonOrthCorr)
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
-            phi += pEqn.flux();
+            fvScalarMatrix pEqn
+            (
+                fvm::ddt(psi, p)
+              + fvm::div(phid, p, "div(phid,p)")
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            pEqn.solve();
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi == pEqn.flux();
+            }
         }
     }
+    else
+    {
+        phi = fvc::interpolate(rho)*
+            ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvm::ddt(psi, p)
+              + fvc::div(phi)
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            pEqn.solve();
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+
+    // Explicitly relax pressure except for last corrector
+    if (oCorr != nOuterCorr - 1)
+    {
+        p.relax();
+    }
 
 #   include "rhoEqn.H"
 #   include "compressibleContinuityErrs.H"
 
-    // Warning:
-    // rho does not carry working boundary conditions and needs to be updated
-    // strictly according to the thermodynamics package
-    // HJ, 22/Aug/2007
-    thermo.correct();
-    rho = thermo.rho();
-
-    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-
     U -= rUA*fvc::grad(p);
     U.correctBoundaryConditions();
-}
 
+    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    dpdt = fvc::ddt(p);
+}

applications/solvers/engine/sonicTurbDyMEngineFoam/readControls.H

@@ -1,14 +0,0 @@
-#   include "readTimeControls.H"
-#   include "readPISOControls.H"
-
-    bool correctPhi = false;
-    if (piso.found("correctPhi"))
-    {
-        correctPhi = Switch(piso.lookup("correctPhi"));
-    }
-
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }

applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C

@@ -49,8 +49,8 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 
 #   include "createEngineTime.H"
-#   include "createDynamicFvMesh.H"
+#   include "createEngineDynamicMesh.H"
-#   include "readPISOControls.H"
+#   include "readPIMPLEControls.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 #   include "readEngineTimeControls.H"
@@ -67,7 +67,9 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readControls.H"
+#       include "readPIMPLEControls.H"
+#       include "checkTotalVolume.H"
+#       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
 
@@ -75,39 +77,48 @@ int main(int argc, char *argv[])
 
         Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
 
-//      make phi relative
+        // Make flux absolute
-
         phi += meshFlux;
 
         bool meshChanged = mesh.update();
-        reduce(meshChanged, orOp<bool>());
+
+#       include "volContinuity.H"
+
+        mesh.setBoundaryVelocity(U);
 
         if (meshChanged)
         {
             thermo.correct();
-
+            rho = thermo.rho();
-#           include "checkTotalVolume.H"
+            rho.correctBoundaryConditions();
-#           include "compressibleCorrectPhi.H"
-#           include "CourantNo.H"
         }
 
         meshFlux = fvc::interpolate(rho)*fvc::meshPhi(rho, U);
 
-        // Make phi absolute
+        phi = fvc::interpolate(rho)
+            *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-        phi -= meshFlux;
+        DpDt = dpdt + fvc::div(phi/fvc::interpolate(rho), p)
+            - fvc::div(phi/fvc::interpolate(rho) + fvc::meshPhi(U))*p;
 
-#       include "rhoEqn.H"
-
-        // --- SIMPLE loop
-        for (int corr=1; corr<=nCorr; corr++)
         {
+#           include "compressibleCourantNo.H"
+        }
+
+        // Pressure-velocity corrector
+        int oCorr = 0;
+        do
+        {
+#           include "rhoEqn.H"
 #           include "UEqn.H"
 
-#           include "hEqn.H"
+            // --- PISO loop
-
+            for (int corr = 1; corr <= nCorr; corr++)
-#           include "pEqn.H"
+            {
-        }
+#               include "pEqn.H"
+#               include "hEqn.H"
+            }
+        } while (++oCorr < nOuterCorr);
 
         turbulence->correct();
 

applications/utilities/errorEstimation/simpleFoamResidual/simpleFoamResidual.C

@@ -103,12 +103,18 @@ int main(int argc, char *argv[])
                     mesh,
                     IOobject::NO_READ
                 ),
-                mag(fvc::div(phi, U) + fvc::div(turbulence->R()) + fvc::grad(p))
+                mag
+                (
+                    fvc::div(phi, U)
+                  + fvc::div(turbulence->R())
+                  + fvc::grad(p)
+                )
             );
 
             Info<< "uResidual max: " << max(uResidual.internalField())
                 << " mean: "
-                << sum(uResidual.internalField()*mesh.V())/sum(mesh.V()).value()
+                << sum(uResidual.internalField()*mesh.V())/
+                       sum(mesh.V()).value()
                 << endl;
 
             uResidual.write();
@@ -119,7 +125,6 @@ int main(int argc, char *argv[])
         }
 
         Info<< endl;
-            
     }
 
     Info<< "End\n" << endl;

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

@@ -4,6 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
@@ -12,4 +13,5 @@ EXE_LIBS = \
     -ldecompositionMethods \
     -lmeshTools \
     -ldynamicMesh \
+    -ldynamicFvMesh \
     -lautoMesh

applications/utilities/mesh/manipulation/moveDyMEngineMesh/moveDyMEngineMesh.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
 #   include "createEngineTime.H"
 #   include "createEngineDynamicMesh.H"
 
+
     fileName path = runTime.caseName();
     OFstream volFile(path+"/totVol.Cyl");
 

applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C

@@ -54,7 +54,8 @@ int main(int argc, char *argv[])
     fadScalarField b = tf - ten;
     fadScalarField c = ten - tf;
     fadScalarField d = sqr(tf - ten);
-    fadScalarField e = (sqr(tf - bs + 10) + 22.3*sqr(cs - bs + 10));
+    fadScalarField e =
+        (sqr(tf - bs + fadScalar(10)) + 22.3*sqr(cs - bs + fadScalar(10)));
     fadScalarField f = d/e;
 //     fadScalarField g = Foam::atan(f);    // HJ, not prepared
 //     fadScalarField h = Foam::log(f);     // HJ, not prepared
@@ -62,7 +63,11 @@ int main(int argc, char *argv[])
 
     Info << "new tf: " << tf << endl;
 
-    Field<Vector<fadScalar> > tfVector(3, Vector<fadScalar>(1, 2, 3));
+    Field<Vector<fadScalar> > tfVector
+    (
+        3,
+        Vector<fadScalar>(fadScalar(1), fadScalar(2), fadScalar(3))
+    );
 
     Info << "tfVector: " << tfVector << endl;
 

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -144,9 +144,51 @@ int main(int argc, char *argv[])
 
     // Mesh write will be controlled by hand
     meshPtr->write();
+    procMeshes.writeAddressing();
     meshPtr->setMotionWriteOpt(IOobject::NO_WRITE);
     meshPtr->setTopoWriteOpt(IOobject::NO_WRITE);
 
+    // Write cell decomposition
+    if (writeCellDist)
+    {
+        // Write as volScalarField for post-processing
+        Info<< "Writing cellDist to time " << runTime.timeName()
+            << endl;
+        volScalarField cellDist
+        (
+            IOobject
+            (
+                "cellDist",
+                runTime.timeName(),
+                meshPtr(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            meshPtr(),
+            dimensionedScalar("cellDist", dimless, 0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        scalarField& cellDistIn = cellDist.internalField();
+
+        label cellI = 0;
+
+        forAll (procMeshes.meshes(), procI)
+        {
+            for
+            (
+                label i = 0;
+                i < procMeshes.meshes()[procI].nCells();
+                i++
+            )
+            {
+                cellDistIn[cellI] = procI;
+                cellI++;
+            }
+        }
+
+        cellDist.write();
+    }
+
     // Get region prefix for lagrangian
     fileName regionPrefix = "";
     if (regionName != fvMesh::defaultRegion)
@@ -222,6 +264,8 @@ int main(int argc, char *argv[])
             if (writeCellDist)
             {
                 // Write as volScalarField for post-processing
+                Info<< "Writing cellDist to time " << runTime.timeName()
+                    << endl;
                 volScalarField cellDist
                 (
                     IOobject

applications/utilities/preProcessing/setFields/setFields.C

@@ -77,13 +77,13 @@ void setFieldType
         }
         else
         {
-            forAll(selectedCells, celli)
+            forAll (selectedCells, celli)
             {
                 field[selectedCells[celli]] = value;
             }
         }
 
-        forAll(field.boundaryField(), patchi)
+        forAll (field.boundaryField(), patchi)
         {
             // Forced patch assignment.  HJ, 1/Aug/2010
             field.boundaryField()[patchi] ==
@@ -176,13 +176,10 @@ public:
 
 int main(int argc, char *argv[])
 {
-    timeSelector::addOptions();
-
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    // Get times list
+    Info<< "Time = " << runTime.timeName() << endl;
-    instantList timeDirs = timeSelector::select0(runTime, args);
 
 #   include "createMesh.H"
 
@@ -216,7 +213,7 @@ int main(int argc, char *argv[])
 
     PtrList<entry> regions(setFieldsDict.lookup("regions"));
 
-    forAll(regions, regionI)
+    forAll (regions, regionI)
     {
         const entry& region = regions[regionI];
 

etc/controlDict

@@ -40,6 +40,8 @@ DebugSwitches
     mixingPlane                           0;
     MixingPlaneInterpolation              0;
 
+    tetFemVectorMatrix        0;
+
     overlapGgi          0;
     cyclicGgi           0;
     coupledLduMatrix    1;

src/OpenFOAM/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C

@@ -466,9 +466,12 @@ initEvaluate
     const Pstream::commsTypes commsType
 )
 {
-    if (this->isPointField())
+    if (Pstream::parRun())
     {
-        initAddFieldTempl(Pstream::blocking, this->internalField());
+        if (this->isPointField())
+        {
+            initAddFieldTempl(Pstream::blocking, this->internalField());
+        }
     }
 }
 
@@ -491,27 +494,30 @@ evaluate
     const Pstream::commsTypes commsType
 )
 {
-    if (this->isPointField())
+    if (Pstream::parRun())
     {
-        // Get the neighbour side values
+        if (this->isPointField())
-        tmp<Field<Type> > tpNeighbour = receivePointField<Type>(commsType);
-        Field<Type>& tpn = tpNeighbour();
-
-        if (doTransform())
         {
-            const processorPolyPatch& ppp = procPatch_.procPolyPatch();
+            // Get the neighbour side values
-            const tensorField& forwardT = ppp.forwardT();
+            tmp<Field<Type> > tpNeighbour = receivePointField<Type>(commsType);
+            Field<Type>& tpn = tpNeighbour();
 
-            transform(tpn, forwardT[0], tpn);
+            if (doTransform())
+            {
+                const processorPolyPatch& ppp = procPatch_.procPolyPatch();
+                const tensorField& forwardT = ppp.forwardT();
+
+                transform(tpn, forwardT[0], tpn);
+            }
+
+            // Average over two sides
+            tpn = 0.5*(patchInternalField(this->internalField()) + tpn);
+
+            // Get internal field to insert values into
+            Field<Type>& iF = const_cast<Field<Type>&>(this->internalField());
+
+            this->setInInternalField(iF, tpn);
         }
-
-        // Average over two sides
-        tpn = 0.5*(patchInternalField(this->internalField()) + tpn);
-
-        // Get internal field to insert values into
-        Field<Type>& iF = const_cast<Field<Type>&>(this->internalField());
-
-        this->setInInternalField(iF, tpn);
     }
 }
 
@@ -697,7 +703,7 @@ initAddUpperLower
 {
     // Gather the data from the given field and send it. Note that the
     // size of the field is equal to the number of edges on the field
-    // and NOT the number of points.  
+    // and NOT the number of points.  HJ, 14/Nov/2001
 
     // Get the addressing
     const labelList& me = procPatch_.localEdgeIndices();

src/OpenFOAM/graph/curve/curve.H

@@ -144,6 +144,7 @@ public:
             const scalarField& y
         );
 
+        //- Create and return a clone
         autoPtr<curve> clone() const
         {
             return autoPtr<curve>(new curve(*this));

src/OpenFOAM/memory/autoPtr/autoPtr.H

@@ -95,7 +95,8 @@ public:
             //  If object pointer already set issue a FatalError.
             inline void set(T*);
 
-            //- If object pointer already set, delete object and set to given pointer
+            //- If object pointer already set, delete object and set to
+            //  given pointer
             inline void reset(T* = 0);
 
             //- If object pointer points to valid object:

src/autoMesh/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/decompositionMethods/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
@@ -10,6 +11,7 @@ EXE_INC = \
 LIB_LIBS = \
     -ldecompositionMethods \
     -ldynamicMesh \
+    -ldynamicFvMesh \
     -lfiniteVolume \
     -llagrangian \
     -lmeshTools \

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -2975,7 +2975,8 @@ void Foam::autoLayerDriver::addLayers
             );
 
             const_cast<Time&>(mesh.time())++;
-            Info<< "Writing shrunk mesh to " << meshRefiner_.timeName() << endl;
+            Info<< "Writing shrunk mesh to " << meshRefiner_.timeName()
+                << endl;
 
             // See comment in autoSnapDriver why we should not remove meshPhi
             // using mesh.clearPout().

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -1613,7 +1613,7 @@ Foam::label Foam::meshRefinement::addPatch
     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
 
     // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
+    polyPatches.setSize(sz + 1);
     polyPatches.set
     (
         sz,

src/decompositionMethods/decompositionMethods/Make/files

@@ -1,7 +1,8 @@
 decompositionMethod/decompositionMethod.C
+manualDecomp/manualDecomp.C
 geomDecomp/geomDecomp.C
 simpleGeomDecomp/simpleGeomDecomp.C
 hierarchGeomDecomp/hierarchGeomDecomp.C
-manualDecomp/manualDecomp.C
+patchConstrainedDecomp/patchConstrainedDecomp.C
 
 LIB = $(FOAM_LIBBIN)/libdecompositionMethods

src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C

@@ -28,6 +28,9 @@ InClass
 \*---------------------------------------------------------------------------*/
 
 #include "decompositionMethod.H"
+#include "cyclicPolyPatch.H"
+#include "syncTools.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -38,7 +41,336 @@ namespace Foam
     defineRunTimeSelectionTable(decompositionMethod, dictionaryMesh);
 }
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+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::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;
+            }
+        }
+    }
+}
+
+
+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++;
+            }
+        }
+    }
+}
+
+
+void Foam::decompositionMethod::calcCellCells
+(
+    const polyMesh& mesh,
+    const labelList& fineToCoarse,
+    const label nCoarse,
+    labelListList& cellCells
+)
+{
+    if (fineToCoarse.size() != mesh.nCells())
+    {
+        FatalErrorIn
+        (
+            "decompositionMethod::calcCellCells"
+            "(const labelList&, labelListList&) const"
+        )   << "Only valid for mesh agglomeration." << exit(FatalError);
+    }
+
+    List<DynamicList<label> > dynCellCells(nCoarse);
+
+    forAll(mesh.faceNeighbour(), faceI)
+    {
+        label own = fineToCoarse[mesh.faceOwner()[faceI]];
+        label nei = fineToCoarse[mesh.faceNeighbour()[faceI]];
+
+        if (own != nei)
+        {
+            if (findIndex(dynCellCells[own], nei) == -1)
+            {
+                dynCellCells[own].append(nei);
+            }
+            if (findIndex(dynCellCells[nei], own) == -1)
+            {
+                dynCellCells[nei].append(own);
+            }
+        }
+    }
+
+    cellCells.setSize(dynCellCells.size());
+    forAll(dynCellCells, coarseI)
+    {
+        cellCells[coarseI].transfer(dynCellCells[coarseI]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 (
@@ -102,17 +434,40 @@ Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const pointField& points
 )
 {
-    scalarField weights(0);
+    scalarField weights(points.size(), 1);
 
     return decompose(points, weights);
 }
 
 
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const labelList& fineToCoarse,
+    const pointField& coarsePoints
+)
+{
+    // Decompose based on agglomerated points
+    labelList coarseDistribution(decompose(coarsePoints));
+
+    // Rework back into decomposition for original mesh_
+    labelList fineDistribution(fineToCoarse.size());
+
+    forAll(fineDistribution, i)
+    {
+        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
+    }
+
+    return fineDistribution;
+}
+
+
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const labelList& fineToCoarse,
@@ -135,82 +490,4 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const labelList& fineToCoarse,
-    const pointField& coarsePoints
-)
-{
-    // Decompose based on agglomerated points
-    labelList coarseDistribution(decompose(coarsePoints));
-
-    // Rework back into decomposition for original mesh_
-    labelList fineDistribution(fineToCoarse.size());
-
-    forAll(fineDistribution, i)
-    {
-        fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
-    }
-
-    return fineDistribution;
-}
-
-
-void Foam::decompositionMethod::calcCellCells
-(
-    const polyMesh& mesh,
-    const labelList& fineToCoarse,
-    const label nCoarse,
-    labelListList& cellCells
-)
-{
-    if (fineToCoarse.size() != mesh.nCells())
-    {
-        FatalErrorIn
-        (
-            "decompositionMethod::calcCellCells"
-            "(const labelList&, labelListList&) const"
-        )   << "Only valid for mesh agglomeration." << exit(FatalError);
-    }
-
-    List<DynamicList<label> > dynCellCells(nCoarse);
-
-    forAll(mesh.faceNeighbour(), faceI)
-    {
-        label own = fineToCoarse[mesh.faceOwner()[faceI]];
-        label nei = fineToCoarse[mesh.faceNeighbour()[faceI]];
-
-        if (own != nei)
-        {
-            if (findIndex(dynCellCells[own], nei) == -1)
-            {
-                dynCellCells[own].append(nei);
-            }
-            if (findIndex(dynCellCells[nei], own) == -1)
-            {
-                dynCellCells[nei].append(own);
-            }
-        }
-    }
-
-    cellCells.setSize(dynCellCells.size());
-    forAll(dynCellCells, coarseI)
-    {
-        cellCells[coarseI].transfer(dynCellCells[coarseI]);
-    }
-}
-
-
-Foam::labelList Foam::decompositionMethod::decompose
-(
-    const labelListList& globalCellCells,
-    const pointField& cc
-)
-{
-    scalarField cWeights(0);
-
-    return decompose(globalCellCells, cc, cWeights);
-}
-
-
 // ************************************************************************* //

src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H

@@ -48,30 +48,64 @@ namespace Foam
 
 class decompositionMethod
 {
-
 protected:
 
     // Protected data
 
+        //- Decomposition dictionary
         const dictionary& decompositionDict_;
+
+        //- Number of processors in decomposition
         label nProcessors_;
 
 
-        //- Helper: determine (non-parallel) cellCells from mesh agglomeration.
+        //- Helper to convert connectivity supplied as cellCells into
+        //  simple CSR (Metis, scotch) storage
+        static void calcCSR
+        (
+            const labelListList& globalCellCells,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper: convert local connectivity from the mesh
+        //  into CSR (Metis, scotch) storage
+        //  Treats cyclics as coupled, but not processor patches
+        static void calcCSR
+        (
+            const polyMesh& mesh,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper: convert mesh connectivity into distributed CSR
+        //  Very dubious coding.  HJ, 1/Mar/2011
+        static void calcDistributedCSR
+        (
+            const polyMesh&,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper: determine (non-parallel) cellCells from mesh
+        // with agglomeration
         static void calcCellCells
         (
             const polyMesh& mesh,
-            const labelList& agglom,
+            const labelList& fineToCoarse,
             const label nCoarse,
             labelListList& cellCells
         );
 
+
 private:
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct and assignment
+        //- Disallow default bitwise copy construct
         decompositionMethod(const decompositionMethod&);
+
+        //- Disallow default bitwise assignment
         void operator=(const decompositionMethod&);
 
 
@@ -109,7 +143,7 @@ public:
 
     // Selectors
 
-        //- Return a reference to the selected decomposition method
+        //- Return a pointer to the selected decomposition method
         static autoPtr<decompositionMethod> New
         (
             const dictionary& decompositionDict
@@ -148,46 +182,37 @@ public:
         //  proc boundaries)
         virtual bool parallelAware() const = 0;
 
-        //- Return for every coordinate the wanted processor number. Use the
+        //- Decompose cells
-        //  mesh connectivity (if needed)
+        //  If needed, use connectivity directly from the mesh
+        //  Calls decompose (below) with uniform weights
+        virtual labelList decompose(const pointField&);
+
+        //- Decompose cells with weights
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         ) = 0;
 
-        //- Like decompose but with uniform weights on the points
-        virtual labelList decompose(const pointField&);
 
-
+        //- Decompose cell clusters
-        //- Return for every coordinate the wanted processor number. Gets
+        //  Calls decompose (below) with uniform weights
-        //  passed agglomeration map (from fine to coarse cells) and coarse cell
-        //  location. Can be overridden by decomposers that provide this
-        //  functionality natively. Coarse cells are local to the processor
-        //  (if in parallel). If you want to have coarse cells spanning
-        //  processors use the globalCellCells instead.
         virtual labelList decompose
         (
-            const labelList& cellToRegion,
+            const labelList& fineToCoarse,
-            const pointField& regionPoints,
+            const pointField& coarsePoints
-            const scalarField& regionWeights
         );
 
-        //- Like decompose but with uniform weights on the regions
+        //- Decompose cell clusters with weights on clusters
         virtual labelList decompose
         (
-            const labelList& cellToRegion,
+            const labelList& fineToCoarse,
-            const pointField& regionPoints
+            const pointField& coarsePoints,
+            const scalarField& coarseWeights
         );
 
 
-        //- Return for every coordinate the wanted processor number. Explicitly
+        //- Decompose cells with weights with explicitly provided connectivity
-        //  provided connectivity - does not use mesh_.
-        //  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,
@@ -195,13 +220,6 @@ public:
             const scalarField& cWeights
         ) = 0;
 
-        //- Like decompose but with uniform weights on the cells
-        virtual labelList decompose
-        (
-            const labelListList& globalCellCells,
-            const pointField& cc
-        );
-
 };
 
 

src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.C

@@ -40,8 +40,7 @@ Foam::geomDecomp::geomDecomp
     delta_(readScalar(geomDecomDict_.lookup("delta"))),
     rotDelta_(I)
 {
-    // check that the case makes sense :
+    // Check that the case makes sense
-
     if (nProcessors_ != n_.x()*n_.y()*n_.z())
     {
         FatalErrorIn

src/decompositionMethods/decompositionMethods/geomDecomp/geomDecomp.H

@@ -37,7 +37,7 @@ SourceFiles
 #define geomDecomp_H
 
 #include "decompositionMethod.H"
-#include "Vector.H"
+#include "labelVector.H"
 
 namespace Foam
 {
@@ -50,17 +50,32 @@ class geomDecomp
 :
     public decompositionMethod
 {
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        geomDecomp(const geomDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const geomDecomp&);
+
 
 protected:
 
     // Protected data
 
+        //- Geometric decomposition dictionary
         const dictionary& geomDecomDict_;
 
-        Vector<label> n_;
+        //- Number of splits in the (x y z) direction
+        labelVector n_;
+
+        //- Small rotation to achieve smooth geometric decomposition
         scalar delta_;
+
+        //- Small rotation tensor, calculated from delta
         tensor rotDelta_;
 
+
 public:
 
     // Constructors

src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C

@@ -50,6 +50,7 @@ namespace Foam
     );
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 void Foam::hierarchGeomDecomp::setDecompOrder()
@@ -157,7 +158,7 @@ void Foam::hierarchGeomDecomp::calculateSortedWeightedSizes
 {
     // Evaluate cumulative weights.
     sortedWeightedSizes[0] = 0;
-    forAll(current, i)
+    forAll (current, i)
     {
         label pointI = current[indices[i]];
         sortedWeightedSizes[i + 1] = sortedWeightedSizes[i] + weights[pointI];
@@ -328,7 +329,7 @@ void Foam::hierarchGeomDecomp::findBinary
 }
 
 
-// Sort points into bins according to one component. Recurses to next component.
+// Sort points into bins according to one component. Recurses to next component
 void Foam::hierarchGeomDecomp::sortComponent
 (
     const label sizeTol,
@@ -351,7 +352,7 @@ void Foam::hierarchGeomDecomp::sortComponent
     // Storage for sorted component compI
     SortableList<scalar> sortedCoord(current.size());
 
-    forAll(current, i)
+    forAll (current, i)
     {
         label pointI = current[i];
 
@@ -458,7 +459,7 @@ void Foam::hierarchGeomDecomp::sortComponent
         // Copy localSize elements starting from leftIndex.
         labelList slice(localSize);
 
-        forAll(slice, i)
+        forAll (slice, i)
         {
             label pointI = current[sortedCoord.indices()[leftIndex+i]];
 
@@ -502,7 +503,7 @@ void Foam::hierarchGeomDecomp::sortComponent
 }
 
 
-// Sort points into bins according to one component. Recurses to next component.
+// Sort points into bins according to one component. Recurses to next component
 void Foam::hierarchGeomDecomp::sortComponent
 (
     const label sizeTol,
@@ -526,7 +527,7 @@ void Foam::hierarchGeomDecomp::sortComponent
     // Storage for sorted component compI
     SortableList<scalar> sortedCoord(current.size());
 
-    forAll(current, i)
+    forAll (current, i)
     {
         label pointI = current[i];
 
@@ -641,7 +642,7 @@ void Foam::hierarchGeomDecomp::sortComponent
         // Copy localSize elements starting from leftIndex.
         labelList slice(localSize);
 
-        forAll(slice, i)
+        forAll (slice, i)
         {
             label pointI = current[sortedCoord.indices()[leftIndex+i]];
 
@@ -725,7 +726,8 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 
     // Start off with every point sorted onto itself.
     labelList slice(points.size());
-    forAll(slice, i)
+
+    forAll (slice, i)
     {
         slice[i] = i;
     }
@@ -767,7 +769,8 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 
     // Start off with every point sorted onto itself.
     labelList slice(points.size());
-    forAll(slice, i)
+
+    forAll (slice, i)
     {
         slice[i] = i;
     }

src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H

@@ -77,6 +77,13 @@ class hierarchGeomDecomp
 
     // Private Member Functions
 
+        //- Disallow default bitwise copy construct
+        hierarchGeomDecomp(const hierarchGeomDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const hierarchGeomDecomp&);
+
+
         //- Convert ordering string ("xyz") into list of components.
         void setDecompOrder();
 
@@ -105,13 +112,13 @@ class hierarchGeomDecomp
         //  wantedSize. Binary search.
         static void findBinary
         (
-            const label sizeTol,        // size difference considered acceptible
+            const label sizeTol,        // Acceptable size difference
             const List<scalar>&,
             const label leftIndex,      // index of previous value
             const scalar leftValue,     // value at leftIndex
             const scalar maxValue,      // global max of values
             const scalar wantedSize,    // wanted size
-            label& mid,                 // index where size of bin is wantedSize
+            label& mid,                 // index where bin size is wantedSize
             scalar& midValue            // value at mid
         );
 
@@ -120,14 +127,14 @@ class hierarchGeomDecomp
         //  wantedSize. Binary search.
         static void findBinary
         (
-            const label sizeTol,        // size difference considered acceptible
+            const label sizeTol,        // Acceptable size difference
             const List<scalar>& sortedWeightedSizes,
             const List<scalar>&,
             const label leftIndex,      // index of previous value
             const scalar leftValue,     // value at leftIndex
             const scalar maxValue,      // global max of values
             const scalar wantedSize,    // wanted size
-            label& mid,                 // index where size of bin is wantedSize
+            label& mid,                 // index where bin size is wantedSize
             scalar& midValue            // value at mid
         );
 
@@ -156,11 +163,6 @@ class hierarchGeomDecomp
         );
 
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const hierarchGeomDecomp&);
-        hierarchGeomDecomp(const hierarchGeomDecomp&);
-
-
 public:
 
     //- Runtime type information
@@ -194,25 +196,19 @@ public:
             return true;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
+        //- Decompose cells without weights. Code for weighted decomposition
-        //  mesh connectivity (if needed)
+        //  is a bit complex so it is kept separate for now
+        virtual labelList decompose(const pointField&);
+
+        //- Decompose cells with weights
         virtual labelList decompose
         (
-            const pointField&,
+            const pointField& points,
             const scalarField& weights
         );
 
-        //- Without weights. Code for weighted decomposition is a bit complex
+        //- Decompose cells with weights with explicitly provided connectivity
-        //  so kept separate for now.
+        //  Does not use mesh for connectivity
-        virtual labelList decompose(const pointField&);
-
-        //- Return for every coordinate the wanted processor number. Explicitly
-        //  provided connectivity - does not use mesh_.
-        //  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,
@@ -223,6 +219,7 @@ public:
             return decompose(cc, cWeights);
         }
 
+        //- Decompose cells with weights with explicitly provided connectivity
         virtual labelList decompose
         (
             const labelListList& globalCellCells,

src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C

@@ -38,13 +38,6 @@ namespace Foam
 {
     defineTypeNameAndDebug(manualDecomp, 0);
 
-    addToRunTimeSelectionTable
-    (
-        decompositionMethod,
-        manualDecomp,
-        dictionary
-    );
-
     addToRunTimeSelectionTable
     (
         decompositionMethod,
@@ -56,14 +49,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::manualDecomp::manualDecomp(const dictionary& decompositionDict)
-:
-    decompositionMethod(decompositionDict)
-{
-    notImplemented("manualDecomp(const dictionary&)");
-}
-
-
 Foam::manualDecomp::manualDecomp
 (
     const dictionary& decompositionDict,
@@ -71,11 +56,14 @@ Foam::manualDecomp::manualDecomp
 )
 :
     decompositionMethod(decompositionDict),
-    meshPtr_(&mesh),
+    mesh_(mesh),
     decompDataFile_
     (
-        decompositionDict.subDict(word(decompositionDict.lookup("method"))
+        decompositionDict.subDict
-      + "Coeffs").lookup("dataFile")
+        (
+            word(decompositionDict.lookup("method"))
+          + "Coeffs"
+        ).lookup("dataFile")
     )
 {}
 
@@ -88,22 +76,20 @@ Foam::labelList Foam::manualDecomp::decompose
     const scalarField& pointWeights
 )
 {
-    const polyMesh& mesh = *meshPtr_;
-
     labelIOList finalDecomp
     (
         IOobject
         (
             decompDataFile_,
-            mesh.facesInstance(),
+            mesh_.facesInstance(),
-            mesh,
+            mesh_,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE,
             false
         )
     );
 
-    // check if the final decomposition is OK
+    // Check if the final decomposition is OK
 
     if (finalDecomp.size() != points.size())
     {

src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H

@@ -42,7 +42,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class manualDecomp Declaration
+                        Class manualDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
 class manualDecomp
@@ -51,17 +51,19 @@ class manualDecomp
 {
     // Private data
 
-        const polyMesh* meshPtr_;
+        const polyMesh& mesh_;
 
         fileName decompDataFile_;
 
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct and assignment
+        //- Disallow default bitwise copy construct
-        void operator=(const manualDecomp&);
         manualDecomp(const manualDecomp&);
 
+        //- Disallow default bitwise assignment
+        void operator=(const manualDecomp&);
+
 
 public:
 
@@ -97,21 +99,14 @@ public:
             return true;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
+        //- Decompose cells with weights
-        //  mesh connectivity (if needed)
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Explicitly
+        //- Decompose cells with weights with explicitly provided connectivity
-        //  provided connectivity - does not use mesh_.
-        //  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,

src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+
+Description
+    Decomposition given a cell-to-processor association in a file
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchConstrainedDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+#include "IFstream.H"
+#include "labelIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchConstrainedDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        patchConstrainedDecomp,
+        dictionaryMesh
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchConstrainedDecomp::patchConstrainedDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    decompositionMethod(decompositionDict),
+    mesh_(mesh),
+    dict_
+    (
+        decompositionDict.subDict
+        (
+            typeName + "Coeffs"
+        )
+    ),
+    baseDecompPtr_
+    (
+        decompositionMethod::New(dict_, mesh)
+    ),
+    patchConstraints_(dict_.lookup("patchConstraints"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::patchConstrainedDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
+{
+
+    labelList finalDecomp = baseDecompPtr_->decompose(points, pointWeights);
+
+    // Impose the decomposition along patches
+    forAll (patchConstraints_, i)
+    {
+        const label patchID =
+            mesh_.boundaryMesh().findPatchID(patchConstraints_[i].first());
+
+        const label procID = patchConstraints_[i].second();
+
+        if (patchID < 0 || procID < 0 || procID > nProcessors_ - 1)
+        {
+            FatalErrorIn
+            (
+                "labelList patchConstrainedDecomp::decompose\n"
+                "(\n"
+                "    const pointField& points,\n"
+                "    const scalarField& pointWeights\n"
+                ")"
+            )   << "Incorrect patch constraint definition for "
+                << patchConstraints_[i]
+                << abort(FatalError);
+        }
+
+        const labelList fc = mesh_.boundaryMesh()[patchID].faceCells();
+
+        forAll (fc, fcI)
+        {
+            finalDecomp[fc[fcI]] = procID;
+        }
+    }
+
+    return finalDecomp;
+}
+
+
+// ************************************************************************* //

src/decompositionMethods/decompositionMethods/patchConstrainedDecomp/patchConstrainedDecomp.H

@@ -22,114 +22,112 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Namespace
-    Foam::solidBodyMotionFunctions
-
-Description
-    Namespace for solid-body motions
-
-
 Class
-    Foam::solidBodyMotionFunction
+    Foam::patchConstrainedDecomp
 
 Description
-    Base class for defining solid-body motions
+    Base decomposition type is modified by assigning patches to given
+    processors.  This is controlled by the patch entry, giving patch names
+    and processor index for them.
 
 SourceFiles
-    solidBodyMotionFunction.C
+    patchConstrainedDecomp.C
-    newDynamicFvMesh.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef solidBodyMotionFunction_H
+#ifndef patchConstrainedDecomp_H
-#define solidBodyMotionFunction_H
+#define patchConstrainedDecomp_H
 
-#include "Time.H"
+#include "decompositionMethod.H"
-#include "dictionary.H"
+#include "Tuple2.H"
-#include "septernion.H"
-#include "autoPtr.H"
-#include "runTimeSelectionTables.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class solidBodyMotionFunction Declaration
+                   Class patchConstrainedDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
-class solidBodyMotionFunction
+class patchConstrainedDecomp
+:
+    public decompositionMethod
 {
-protected:
+    // Private data types
 
-    // Protected data
+    typedef List<Tuple2<word, label> > procWordList;
-
-        dictionary SBMFCoeffs_;
-        const Time& time_;
 
 
-private:
+    // Private data
+
+        //- Mesh
+        const polyMesh& mesh_;
+
+        //- Decomposition dictionary
+        dictionary dict_;
+
+        //- Base decomposition method
+        autoPtr<decompositionMethod> baseDecompPtr_;
+
+        //- Patch names to processor association
+        procWordList patchConstraints_;
+
 
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        solidBodyMotionFunction(const solidBodyMotionFunction&);
+        patchConstrainedDecomp(const patchConstrainedDecomp&);
 
-        //- Disallow default bitwise assignment
+        //- Disallow default bitwise copy construct and assignment
-        void operator=(const solidBodyMotionFunction&);
+        void operator=(const patchConstrainedDecomp&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("solidBodyMotionFunction");
+    TypeName("patchConstrained");
-
-
-    // Declare run-time constructor selection table
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            solidBodyMotionFunction,
-            dictionary,
-            (const dictionary& SBMFCoeffs, const Time& runTime),
-            (SBMFCoeffs, runTime)
-        );
 
 
     // Constructors
 
-        //- Construct from the SBMFCoeffs dictionary and Time
+        //- Construct given the decomposition dictionary and mesh
-        solidBodyMotionFunction
+        patchConstrainedDecomp
         (
-            const dictionary& SBMFCoeffs,
+            const dictionary& decompositionDict,
-            const Time& runTime
+            const polyMesh& mesh
-        );
-
-
-    // Selectors
-
-        //- Select constructed from the SBMFCoeffs dictionary and Time
-        static autoPtr<solidBodyMotionFunction> New
-        (
-            const dictionary& SBMFCoeffs,
-            const Time& runTime
         );
 
 
     // Destructor
 
-        virtual ~solidBodyMotionFunction();
+        virtual ~patchConstrainedDecomp()
+        {}
 
 
     // Member Functions
 
-        //- Return the solid-body motion transformation septernion
+        //- Return parallel aware of base method
-        virtual septernion transformation() const = 0;
+        virtual bool parallelAware() const
+        {
+            return baseDecompPtr_->parallelAware();
+        }
 
-        //- Update properties from given dictionary
+        //- Decompose cells with weights
-        virtual bool read(const dictionary& SBMFCoeffs) = 0;
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
+
+        //- Decompose cells with weights with explicitly provided connectivity
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights
+        )
+        {
+            return decompose(cc, cWeights);
+        }
 };
 
 

src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C

@@ -72,9 +72,9 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 
     // assign cells to the first few processor groups (those with
     // one extra cell each
-    for (j=0; j<fstProcessorGroup; j++)
+    for (j = 0; j < fstProcessorGroup; j++)
     {
-        for (register label k=0; k<jumpb; k++)
+        for (register label k = 0; k < jumpb; k++)
         {
             processorGroup[ind++] = j;
         }
@@ -83,7 +83,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     // and now to the `normal' processor groups
     for (; j<nProcGroup; j++)
     {
-        for (register label k=0; k<jump; k++)
+        for (register label k = 0; k < jump; k++)
         {
             processorGroup[ind++] = j;
         }
@@ -316,4 +316,6 @@ Foam::labelList Foam::simpleGeomDecomp::decompose
 
     return finalDecomp;
 }
+
+
 // ************************************************************************* //

src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H

@@ -41,7 +41,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class simpleGeomDecomp Declaration
+                      Class simpleGeomDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
 class simpleGeomDecomp
@@ -61,10 +61,12 @@ class simpleGeomDecomp
             const scalar summedWeights
         );
 
-        //- Disallow default bitwise copy construct and assignment
+        //- Disallow default bitwise copy construct
-        void operator=(const simpleGeomDecomp&);
         simpleGeomDecomp(const simpleGeomDecomp&);
 
+        //- Disallow default bitwise assignment
+        void operator=(const simpleGeomDecomp&);
+
 
 public:
 
@@ -75,7 +77,7 @@ public:
     // Constructors
 
         //- Construct given the decomposition dictionary
-        simpleGeomDecomp(const dictionary& decompositionDict);
+        explicit simpleGeomDecomp(const dictionary& decompositionDict);
 
         //- Construct given the decomposition dictionary and mesh
         simpleGeomDecomp

src/decompositionMethods/metisDecomp/metisDecomp.C

@@ -54,13 +54,11 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Call Metis with options from dictionary.
 Foam::label Foam::metisDecomp::decompose
 (
     const List<int>& adjncy,
     const List<int>& xadj,
     const scalarField& cWeights,
-
     List<int>& finalDecomp
 )
 {
@@ -133,7 +131,8 @@ Foam::label Foam::metisDecomp::decompose
             if (method != "recursive" && method != "k-way")
             {
                 FatalErrorIn("metisDecomp::decompose()")
-                    << "Method " << method << " in metisCoeffs in dictionary : "
+                    << "Method " << method
+                    << " in metisCoeffs in dictionary : "
                     << decompositionDict_.name()
                     << " should be 'recursive' or 'k-way'"
                     << exit(FatalError);
@@ -363,18 +362,18 @@ Foam::labelList Foam::metisDecomp::decompose
 
 Foam::labelList Foam::metisDecomp::decompose
 (
-    const labelList& agglom,
+    const labelList& fineToCoarse,
-    const pointField& agglomPoints,
+    const pointField& coarsePoints,
-    const scalarField& agglomWeights
+    const scalarField& coarseWeights
 )
 {
-    if (agglom.size() != mesh_.nCells())
+    if (fineToCoarse.size() != mesh_.nCells())
     {
         FatalErrorIn
         (
             "metisDecomp::decompose"
             "(const labelList&, const pointField&, const scalarField&)"
-        )   << "Size of cell-to-coarse map " << agglom.size()
+        )   << "Size of cell-to-coarse map " << fineToCoarse.size()
             << " differs from number of cells in mesh " << mesh_.nCells()
             << exit(FatalError);
     }
@@ -390,8 +389,8 @@ Foam::labelList Foam::metisDecomp::decompose
         calcCellCells
         (
             mesh_,
-            agglom,
+            fineToCoarse,
-            agglomPoints.size(),
+            coarsePoints.size(),
             cellCells
         );
 
@@ -400,15 +399,15 @@ Foam::labelList Foam::metisDecomp::decompose
 
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, agglomWeights, finalDecomp);
+    decompose(adjncy, xadj, coarseWeights, finalDecomp);
 
 
     // Rework back into decomposition for original mesh_
-    labelList fineDistribution(agglom.size());
+    labelList fineDistribution(fineToCoarse.size());
 
     forAll(fineDistribution, i)
     {
-        fineDistribution[i] = finalDecomp[agglom[i]];
+        fineDistribution[i] = finalDecomp[fineToCoarse[i]];
     }
 
     return fineDistribution;
@@ -418,18 +417,18 @@ Foam::labelList Foam::metisDecomp::decompose
 Foam::labelList Foam::metisDecomp::decompose
 (
     const labelListList& globalCellCells,
-    const pointField& cellCentres,
+    const pointField& cc,
-    const scalarField& cellWeights
+    const scalarField& cWeights
 )
 {
-    if (cellCentres.size() != globalCellCells.size())
+    if (cc.size() != globalCellCells.size())
     {
         FatalErrorIn
         (
             "metisDecomp::decompose"
             "(const pointField&, const labelListList&, const scalarField&)"
         )   << "Inconsistent number of cells (" << globalCellCells.size()
-            << ") and number of cell centres (" << cellCentres.size()
+            << ") and number of cell centres (" << cc.size()
             << ")." << exit(FatalError);
     }
 
@@ -442,17 +441,18 @@ Foam::labelList Foam::metisDecomp::decompose
     List<int> xadj;
     scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
 
-
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, cellWeights, finalDecomp);
+    decompose(adjncy, xadj, cWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
+
     forAll(decomp, i)
     {
         decomp[i] = finalDecomp[i];
     }
+
     return decomp;
 }
 

src/decompositionMethods/metisDecomp/metisDecomp.H

@@ -56,6 +56,12 @@ class metisDecomp
 
     // Private Member Functions
 
+        //- Disallow default bitwise copy construct
+        metisDecomp(const metisDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const metisDecomp&);
+
         label decompose
         (
             const List<int>& adjncy,
@@ -64,10 +70,6 @@ class metisDecomp
             List<int>& finalDecomp
         );
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const metisDecomp&);
-        metisDecomp(const metisDecomp&);
-
 
 public:
 
@@ -99,44 +101,28 @@ public:
             return false;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
+        //- Decompose cells with weights
-        //  mesh connectivity (if needed)
-        //  Weights get normalised so the minimum value is 1 before truncation
-        //  to an integer so the weights should be multiples of the minimum
-        //  value. The overall sum of weights might otherwise overflow.
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Gets
+        //- Decompose cell clusters with weights on clusters
-        //  passed agglomeration map (from fine to coarse cells) and coarse cell
-        //  location. Can be overridden by decomposers that provide this
-        //  functionality natively.
-        //  See note on weights above.
         virtual labelList decompose
         (
-            const labelList& agglom,
+            const labelList& fineToCoarse,
-            const pointField& regionPoints,
+            const pointField& coarsePoints,
-            const scalarField& regionWeights
+            const scalarField& coarseWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Explicitly
+        //- Decompose cells with weights with explicitly provided connectivity
-        //  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
-        //  See note on weights above.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
             const pointField& cc,
             const scalarField& cWeights
         );
-
 };
 
 

src/decompositionMethods/parMetisDecomp/parMetisDecomp.C

@@ -26,13 +26,11 @@ License
 
 #include "parMetisDecomp.H"
 #include "metisDecomp.H"
-#include "scotchDecomp.H"
-#include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
-#include "polyMesh.H"
 #include "Time.H"
 #include "labelIOField.H"
+#include "syncTools.H"
 #include "globalIndex.H"
 
 #include <mpi.h>
@@ -68,7 +66,6 @@ Foam::label Foam::parMetisDecomp::decompose
     Field<int>& cellWeights,
     Field<int>& faceWeights,
     const List<int>& options,
-
     List<int>& finalDecomp
 )
 {
@@ -406,7 +403,7 @@ Foam::labelList Foam::parMetisDecomp::decompose
     Field<int> adjncy;
     // Offsets into adjncy
     Field<int> xadj;
-    calcMetisDistributedCSR
+    calcDistributedCSR
     (
         mesh_,
         adjncy,
@@ -594,7 +591,6 @@ Foam::labelList Foam::parMetisDecomp::decompose
         cellWeights,
         faceWeights,
         options,
-
         finalDecomp
     );
 
@@ -706,7 +702,11 @@ Foam::labelList Foam::parMetisDecomp::decompose
                     label globalNei = globalNeighbour[bFaceI++];
                     faceI++;
 
-                    if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
+                    if
+                    (
+                        findIndex(dynRegionRegions[ownRegion], globalNei)
+                     == -1
+                    )
                     {
                         dynRegionRegions[ownRegion].append(globalNei);
                     }
@@ -738,6 +738,7 @@ Foam::labelList Foam::parMetisDecomp::decompose
     {
         cellDistribution[cellI] = regionDecomp[cellToRegion[cellI]];
     }
+
     return cellDistribution;
 }
 
@@ -763,8 +764,12 @@ Foam::labelList Foam::parMetisDecomp::decompose
     // For running sequential ...
     if (Pstream::nProcs() <= 1)
     {
-        return metisDecomp(decompositionDict_, mesh_)
+        return metisDecomp(decompositionDict_, mesh_).decompose
-            .decompose(globalCellCells, cellCentres, cWeights);
+        (
+            globalCellCells,
+            cellCentres,
+            cWeights
+        );
     }
 
 
@@ -772,9 +777,11 @@ Foam::labelList Foam::parMetisDecomp::decompose
 
     // Connections
     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);
@@ -857,7 +864,6 @@ Foam::labelList Foam::parMetisDecomp::decompose
         cellWeights,
         faceWeights,
         options,
-
         finalDecomp
     );
 
@@ -871,146 +877,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/decompositionMethods/parMetisDecomp/parMetisDecomp.H

@@ -41,7 +41,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class parMetisDecomp Declaration
+                        Class parMetisDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
 class parMetisDecomp
@@ -55,10 +55,18 @@ class parMetisDecomp
 
     // Private Member Functions
 
-        //- Insert list in front of list.
+        //- Disallow default bitwise copy construct
+        parMetisDecomp(const parMetisDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const parMetisDecomp&);
+
+
+        //- Insert list in front of list
         template<class Type>
         static void prepend(const UList<Type>&, List<Type>&);
-        //- Insert list at end of list.
+
+        //- Insert list at end of list
         template<class Type>
         static void append(const UList<Type>&, List<Type>&);
 
@@ -75,11 +83,6 @@ class parMetisDecomp
         );
 
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const parMetisDecomp&);
-        parMetisDecomp(const parMetisDecomp&);
-
-
 public:
 
     //- Runtime type information
@@ -110,22 +113,14 @@ public:
             return true;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
+        //- Decompose cells with weights
-        //  mesh connectivity (if needed)
-        //  Weights get normalised so the minimum value is 1 before truncation
-        //  to an integer so the weights should be multiples of the minimum
-        //  value. The overall sum of weights might otherwise overflow.
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Gets
+        //- Decompose cell clusters with weights on clusters
-        //  passed agglomeration map (from fine to coarse cells) and coarse cell
-        //  location. Can be overridden by decomposers that provide this
-        //  functionality natively.
-        //  See note on weights above.
         virtual labelList decompose
         (
             const labelList& cellToRegion,
@@ -133,28 +128,13 @@ public:
             const scalarField& regionWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Explicitly
+        //- Decompose cells with weights with explicitly provided connectivity
-        //  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
-        //  See note on weights above.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
             const pointField& cc,
             const scalarField& cWeights
         );
-
-        //- Helper to convert mesh connectivity into distributed CSR
-        static void calcMetisDistributedCSR
-        (
-            const polyMesh&,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
 };
 
 

src/decompositionMethods/scotchDecomp/Make/files

@@ -1,3 +1,4 @@
-scotchDecomp.C
+scotchDecomp/scotchDecomp.C
+engineScotchDecomp/engineScotchDecomp.C
 
 LIB = $(FOAM_LIBBIN)/libscotchDecomp

src/decompositionMethods/scotchDecomp/engineScotchDecomp/engineScotchDecomp.C

@@ -0,0 +1,459 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 "engineScotchDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(engineScotchDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        scotchDecomp,
+        engineScotchDecomp,
+        dictionaryMesh
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::engineScotchDecomp::engineScotchDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    scotchDecomp(decompositionDict, mesh),
+    dict_(decompositionDict.subDict(typeName + "Coeffs")),
+    slidingPatchPairs_(dict_.lookup("slidingPatchPairs")),
+    expandSliding_(dict_.lookup("expandSliding"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::engineScotchDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
+{
+    if (points.size() != mesh().nCells())
+    {
+        FatalErrorIn
+        (
+            "engineScotchDecomp::decompose\n"
+            "(\n"
+            "    const pointField&,\n"
+            "    const scalarField&\n"
+            ")"
+        )   << "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);
+    }
+
+    // Create clustering to coarse level
+
+    labelList fineToCoarse(mesh().nCells(), -1);
+
+    // Mask all cells in the cylinder layering region with true
+    // Used in two-pass decomposition later
+    boolList cylinderMask(mesh().nCells(), false);
+    label nClusters = 0;
+
+    // Locate piston patch and cluster up colums, using opposite faces
+    label pistonPatchID = mesh().boundaryMesh().findPatchID("piston");
+
+    // Go through the sliding pairs and mark clustering
+    forAll (slidingPatchPairs_, pairI)
+    {
+        // Locate patch and cluster up colums, using opposite faces
+        label firstPatchID = mesh().boundaryMesh().findPatchID
+        (
+            slidingPatchPairs_[pairI].first()
+        );
+
+        label secondPatchID = mesh().boundaryMesh().findPatchID
+        (
+            slidingPatchPairs_[pairI].second()
+        );
+
+        if (firstPatchID == -1 || secondPatchID == -1)
+        {
+            FatalErrorIn
+            (
+                "labelList engineScotchDecomp::decompose\n"
+                "(\n"
+                "    const pointField& points,\n"
+                "    const scalarField& pointWeights\n"
+                ")"
+            )   << "Cannot find sliding patch pair "
+                << slidingPatchPairs_[pairI]
+                << abort(FatalError);
+        }
+
+        // Put all cells next to the patch into a cluster
+        if (expandSliding_)
+        {
+            // Use point-cell addressing from patch points
+            const labelListList& pc = mesh().pointCells();
+
+            // First side
+            const labelList& mp1 =
+                mesh().boundaryMesh()[firstPatchID].meshPoints();
+
+            forAll (mp1, pointI)
+            {
+                const labelList& curCells = pc[mp1[pointI]];
+
+                forAll (curCells, cellI)
+                {
+                    fineToCoarse[curCells[cellI]] = nClusters;
+                    cylinderMask[curCells[cellI]] = true;
+                }
+            }
+
+            // Second side
+            {
+                const labelList& mp2 =
+                    mesh().boundaryMesh()[secondPatchID].meshPoints();
+
+                forAll (mp2, pointI)
+                {
+                    const labelList& curCells = pc[mp2[pointI]];
+
+                    forAll (curCells, cellI)
+                    {
+                        fineToCoarse[curCells[cellI]] = nClusters;
+                        cylinderMask[curCells[cellI]] = true;
+                    }
+                }
+            }
+        }
+        else
+        {
+            // First side
+            {
+                const labelList& fc1 =
+                    mesh().boundaryMesh()[firstPatchID].faceCells();
+
+                forAll (fc1, fcI)
+                {
+                    fineToCoarse[fc1[fcI]] = nClusters;
+                    cylinderMask[fc1[fcI]] = true;
+                }
+            }
+
+            // Second side
+            {
+                const labelList& fc2 =
+                    mesh().boundaryMesh()[secondPatchID].faceCells();
+
+                forAll (fc2, fcI)
+                {
+                    fineToCoarse[fc2[fcI]] = nClusters;
+                    cylinderMask[fc2[fcI]] = true;
+                }
+            }
+        }
+
+        nClusters++;
+    }
+
+    if (pistonPatchID > -1)
+    {
+        // Found piston patch
+        const faceList& f = mesh().allFaces();
+        const cellList& c = mesh().cells();
+        const labelList& owner = mesh().faceOwner();
+        const labelList& neighbour = mesh().faceNeighbour();
+
+        const labelList& faceCells =
+            mesh().boundaryMesh()[pistonPatchID].faceCells();
+
+        forAll (faceCells, faceI)
+        {
+            // Get face index
+            label curFaceNo = faceI
+                + mesh().boundaryMesh()[pistonPatchID].start();
+
+            // Get cell index
+            label curCellNo = faceCells[faceI];
+
+            // Mark cell to cluster
+            if (fineToCoarse[curCellNo] < 0)
+            {
+                // New cluster
+                fineToCoarse[curCellNo] = nClusters;
+                cylinderMask[curCellNo] = true;
+
+                for(;;)
+                {
+                    // Attempt to find the next face and cell
+                    curFaceNo = c[curCellNo].opposingFaceLabel(curFaceNo, f);
+
+                    if (curFaceNo > -1)
+                    {
+                        // Face found, try for a cell
+                        if (curFaceNo < mesh().nInternalFaces())
+                        {
+                            if (owner[curFaceNo] == curCellNo)
+                            {
+                                curCellNo = neighbour[curFaceNo];
+                            }
+                            else if (neighbour[curFaceNo] == curCellNo)
+                            {
+                                curCellNo = owner[curFaceNo];
+                            }
+                            else
+                            {
+                                // Error in layering.  Should never happen
+                                break;
+                            }
+
+                            // Mark cell to cluster
+                            fineToCoarse[curCellNo] = nClusters;
+                            cylinderMask[curCellNo] = true;
+                        }
+                        else
+                        {
+                            // Hit boundary face
+                            break;
+                        }
+                    }
+                    else
+                    {
+                        // Cannot find opposing face: out of prismatic region
+                        break;
+                    }
+                }
+
+                // Next cluster
+                nClusters++;
+            }
+        }
+    }
+
+    // Count cylinder cells from mask
+    label nCylinderCells = 0;
+
+    forAll (cylinderMask, cellI)
+    {
+        if (cylinderMask[cellI]) nCylinderCells++;
+    }
+
+    label nStaticCells = mesh().nCells() - nCylinderCells;
+    label nCylClusters = nClusters;
+
+    // Visit all unmarked cells and put them into single clusters
+    forAll (fineToCoarse, cellI)
+    {
+        if (fineToCoarse[cellI] == -1)
+        {
+            fineToCoarse[cellI] = nClusters;
+            nClusters++;
+        }
+    }
+
+    label nStaticClusters = nClusters - nCylClusters;
+
+    Info<< "Number of cells: " << mesh().nCells()
+        << ", in cylinder + sliding: " << nCylinderCells
+        << ", in static part: " << nStaticCells << nl
+        << "Number of cylinder clusters " << nCylClusters
+        << ", static clusters " << nStaticClusters
+        << ", total clusters " << nClusters << endl;
+
+    // Mark-up complete.  Create point centres and weights for all clusters
+    vectorField clusterCentres(nClusters, vector::zero);
+
+    // Stabilise cluster volumes in case a cluster ends up empty
+    // It is possible to have empty clusters without connectivity
+    scalarField clusterVols(nClusters, SMALL);
+    scalarField clusterWeights(nClusters, 0);
+
+    const vectorField& centres = mesh().cellCentres();
+    const scalarField& vols = mesh().cellVolumes();
+
+    forAll (fineToCoarse, cellI)
+    {
+        const label& curCoarse = fineToCoarse[cellI];
+
+        clusterCentres[curCoarse] += centres[cellI]*vols[cellI];
+        clusterVols[curCoarse] += vols[cellI];
+        clusterWeights[curCoarse] += 1;
+    }
+
+    clusterCentres /= clusterVols;
+
+    // Execute decomposition, first on cylinder layering zone, then on the rest
+
+    // Collect cell-cells in cylinder layering zone and the rest
+
+    // Create and cellCells hash lookup on two pieces
+
+    // Note: cell clusters come first and they will be done
+    // on a shortened list.  Static clusters need to be renumbered by
+    // throwing away the first part of the list
+
+    List<labelHashSet> cylCellCellsHash(nCylClusters);
+    List<labelHashSet> staticCellCellsHash(nStaticClusters);
+
+    const labelListList cellCells = mesh().cellCells();
+
+    forAll (cellCells, cellI)
+    {
+        const labelList& curCC = cellCells[cellI];
+
+        label curCluster = fineToCoarse[cellI];
+
+        if (cylinderMask[cellI])
+        {
+            labelHashSet& curCylAddr = cylCellCellsHash[curCluster];
+
+            // Collect neighbour cluster addressing
+            forAll (curCC, neiI)
+            {
+                // Add neighbour if marked
+                if (cylinderMask[curCC[neiI]])
+                {
+                    label nbrCluster = fineToCoarse[curCC[neiI]];
+
+                    if (nbrCluster != curCluster)
+                    {
+                        if (!curCylAddr.found(nbrCluster))
+                        {
+                            curCylAddr.insert(nbrCluster);
+                        }
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Offset index
+            curCluster -= nCylClusters;
+
+            labelHashSet& curStaticAddr = staticCellCellsHash[curCluster];
+
+            forAll (curCC, neiI)
+            {
+                // Add neighbour if marked
+                if (!cylinderMask[curCC[neiI]])
+                {
+                    label nbrCluster = fineToCoarse[curCC[neiI]]
+                        - nCylClusters;
+
+                    if (nbrCluster != curCluster)
+                    {
+                        if (!curStaticAddr.found(nbrCluster))
+                        {
+                            curStaticAddr.insert(nbrCluster);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // Pack cellCells on the cylinder
+    labelListList cylCellCells(nCylClusters);
+
+    forAll (cylCellCellsHash, clusterI)
+    {
+        cylCellCells[clusterI] = cylCellCellsHash[clusterI].toc();
+    }
+
+    // Decompose cylinder: size of list equals the number of clusters
+    // in the cylinder region
+    vectorField clusterCentresCyl = clusterCentres;
+    scalarField clusterWeightsCyl = clusterWeights;
+
+    clusterCentresCyl.setSize(nCylClusters);
+    clusterWeightsCyl.setSize(nCylClusters);
+
+    labelList cylDecomp = scotchDecomp::decompose
+    (
+        cylCellCells,
+        clusterCentresCyl,
+        clusterWeightsCyl
+    );
+
+    // Decompose static: size of list equals the number of clusters
+
+    labelListList staticCellCells(nStaticClusters);
+
+    forAll (staticCellCellsHash, clusterI)
+    {
+        staticCellCells[clusterI] = staticCellCellsHash[clusterI].toc();
+    }
+
+    vectorField clusterCentresStatic(nStaticClusters);
+    scalarField clusterWeightsStatic(nStaticClusters);
+
+    forAll (clusterCentresStatic, i)
+    {
+        clusterCentresStatic[i] = clusterCentres[nCylClusters + i];
+        clusterWeightsStatic[i] = clusterWeights[nCylClusters + i];
+    }
+
+    labelList staticDecomp = scotchDecomp::decompose
+    (
+        staticCellCells,
+        clusterCentresStatic,
+        clusterWeightsStatic
+    );
+
+    // Reconstruct final decomposition
+
+    labelList finalDecomp(mesh().nCells(), -1);
+
+    forAll (cylinderMask, cellI)
+    {
+        if (cylinderMask[cellI])
+        {
+            // Cylinder cell
+            finalDecomp[cellI] = cylDecomp[fineToCoarse[cellI]];
+        }
+        else
+        {
+            // Static cell
+            finalDecomp[cellI] =
+                staticDecomp[fineToCoarse[cellI] - nCylClusters];
+        }
+    }
+
+    return finalDecomp;
+}
+
+
+// ************************************************************************* //

src/decompositionMethods/scotchDecomp/engineScotchDecomp/engineScotchDecomp.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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::engineScotchDecomp
+
+Description
+    Domain decomposition for internal combustion engine simulations with
+    topological changes.
+
+    The algorithm attempts to be as general as possible and operates in the
+    following steps:
+    - identify piston region and peform a decomposition biased towards the
+    piston axis.  This aims to preserve load balancing when cell layers
+    are added
+    - identify static regions and decompose them separately
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    engineScotchDecomp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef engineScotchDecomp_H
+#define engineScotchDecomp_H
+
+#include "scotchDecomp.H"
+#include "Switch.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class engineScotchDecomp Declaration
+\*---------------------------------------------------------------------------*/
+
+class engineScotchDecomp
+:
+    public scotchDecomp
+{
+    // Private data
+
+        //- Dictionary
+        dictionary dict_;
+
+        //- Sliding patch pairs
+        List<Pair<word> > slidingPatchPairs_;
+
+        //- Use expanded area around sliding patches
+        Switch expandSliding_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        engineScotchDecomp(const engineScotchDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const engineScotchDecomp&);
+
+
+        //- Check and print error message
+        static void check(const int, const char*);
+
+        label decompose
+        (
+            const List<int>& adjncy,
+            const List<int>& xadj,
+            const scalarField& cWeights,
+            List<int>& finalDecomp
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("engineScotch");
+
+
+    // Constructors
+
+        //- Construct given the decomposition dictionary and mesh
+        engineScotchDecomp
+        (
+            const dictionary& decompositionDict,
+            const polyMesh& mesh
+        );
+
+
+    // Destructor
+
+        virtual ~engineScotchDecomp()
+        {}
+
+
+    // Member Functions
+
+        virtual bool parallelAware() const
+        {
+            return false;
+        }
+
+        //- Decompose cells with weights
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
+
+        //- Decompose cells with weights with explicitly provided connectivity
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights
+        )
+        {
+            return decompose(cc, cWeights);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/decompositionMethods/scotchDecomp/scotchDecomp/scotchDecomp.C

@@ -109,7 +109,6 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
 #include "Time.H"
-#include "cyclicPolyPatch.H"
 #include "OFstream.H"
 
 extern "C"
@@ -465,17 +464,17 @@ Foam::labelList Foam::scotchDecomp::decompose
 
 Foam::labelList Foam::scotchDecomp::decompose
 (
-    const labelList& agglom,
+    const labelList& fineToCoarse,
-    const pointField& agglomPoints,
+    const pointField& coarsePoints,
-    const scalarField& pointWeights
+    const scalarField& coarseWeights
 )
 {
-    if (agglom.size() != mesh_.nCells())
+    if (fineToCoarse.size() != mesh_.nCells())
     {
         FatalErrorIn
         (
             "parMetisDecomp::decompose(const labelList&, const pointField&)"
-        )   << "Size of cell-to-coarse map " << agglom.size()
+        )   << "Size of cell-to-coarse map " << fineToCoarse.size()
             << " differs from number of cells in mesh " << mesh_.nCells()
             << exit(FatalError);
     }
@@ -488,11 +487,12 @@ Foam::labelList Foam::scotchDecomp::decompose
     {
         // Get cellCells on coarse mesh.
         labelListList cellCells;
+
         calcCellCells
         (
             mesh_,
-            agglom,
+            fineToCoarse,
-            agglomPoints.size(),
+            coarsePoints.size(),
             cellCells
         );
 
@@ -501,14 +501,14 @@ Foam::labelList Foam::scotchDecomp::decompose
 
     // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, pointWeights, finalDecomp);
+    decompose(adjncy, xadj, coarseWeights, finalDecomp);
 
     // Rework back into decomposition for original mesh_
-    labelList fineDistribution(agglom.size());
+    labelList fineDistribution(fineToCoarse.size());
 
     forAll(fineDistribution, i)
     {
-        fineDistribution[i] = finalDecomp[agglom[i]];
+        fineDistribution[i] = finalDecomp[fineToCoarse[i]];
     }
 
     return fineDistribution;
@@ -518,18 +518,18 @@ Foam::labelList Foam::scotchDecomp::decompose
 Foam::labelList Foam::scotchDecomp::decompose
 (
     const labelListList& globalCellCells,
-    const pointField& cellCentres,
+    const pointField& cc,
     const scalarField& cWeights
 )
 {
-    if (cellCentres.size() != globalCellCells.size())
+    if (cc.size() != globalCellCells.size())
     {
         FatalErrorIn
         (
             "scotchDecomp::decompose"
             "(const labelListList&, const pointField&, const scalarField&)"
         )   << "Inconsistent number of cells (" << globalCellCells.size()
-            << ") and number of cell centres (" << cellCentres.size()
+            << ") and number of cell centres (" << cc.size()
             << ")." << exit(FatalError);
     }
 
@@ -548,152 +548,14 @@ Foam::labelList Foam::scotchDecomp::decompose
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
+
     forAll(decomp, i)
     {
         decomp[i] = finalDecomp[i];
     }
+
     return decomp;
 }
 
 
-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/decompositionMethods/scotchDecomp/scotchDecomp/scotchDecomp.H

@@ -42,7 +42,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class scotchDecomp Declaration
+                        Class scotchDecomp Declaration
 \*---------------------------------------------------------------------------*/
 
 class scotchDecomp
@@ -51,11 +51,19 @@ class scotchDecomp
 {
     // Private data
 
+        //- Mesh reference
         const polyMesh& mesh_;
 
 
     // Private Member Functions
 
+        //- Disallow default bitwise copy construct
+        scotchDecomp(const scotchDecomp&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const scotchDecomp&);
+
+
         //- Check and print error message
         static void check(const int, const char*);
 
@@ -67,10 +75,6 @@ class scotchDecomp
             List<int>& finalDecomp
         );
 
-        //- Disallow default bitwise copy construct and assignment
-        void operator=(const scotchDecomp&);
-        scotchDecomp(const scotchDecomp&);
-
 
 public:
 
@@ -96,69 +100,40 @@ public:
 
     // Member Functions
 
+        //- Return mesh
+        const polyMesh& mesh() const
+        {
+            return mesh_;
+        }
+
         virtual bool parallelAware() const
         {
-            // Metis does not know about proc boundaries
+            // Scotch does not know about proc boundaries
             return false;
         }
 
-        //- Return for every coordinate the wanted processor number. Use the
+        //- Decompose cells with weights
-        //  mesh connectivity (if needed)
-        //  Weights get normalised so the minimum value is 1 before truncation
-        //  to an integer so the weights should be multiples of the minimum
-        //  value. The overall sum of weights might otherwise overflow.
         virtual labelList decompose
         (
             const pointField& points,
             const scalarField& pointWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Gets
+        //- Decompose cell clusters with weights on clusters
-        //  passed agglomeration map (from fine to coarse cells) and coarse cell
-        //  location. Can be overridden by decomposers that provide this
-        //  functionality natively.
-        //  See note on weights above.
         virtual labelList decompose
         (
-            const labelList& agglom,
+            const labelList& fineToCoarse,
-            const pointField& regionPoints,
+            const pointField& coarsePoints,
-            const scalarField& regionWeights
+            const scalarField& coarseWeights
         );
 
-        //- Return for every coordinate the wanted processor number. Explicitly
+        //- Decompose cells with weights with explicitly provided connectivity
-        //  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
-        //  See note on weights above.
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
             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
-        );
 };
 
 

src/dynamicMesh/dynamicFvMesh/Make/files

@@ -1,6 +1,8 @@
 dynamicFvMesh/dynamicFvMesh.C
 dynamicFvMesh/newDynamicFvMesh.C
+
 staticFvMesh/staticFvMesh.C
+solidBodyMotionFvMesh/solidBodyMotionFvMesh.C
 dynamicBoxFvMesh/dynamicBoxFvMesh.C
 movingBoxFvMesh/movingBoxFvMesh.C
 dynamicBodyFvMesh/dynamicBodyFvMesh.C
@@ -9,17 +11,12 @@ subsetMotionSolverFvMesh/subsetMotionSolverFvMesh.C
 dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
 dynamicRefineFvMesh/dynamicRefineFvMesh.C
 
-solidBodyMotionFvMesh/solidBodyMotionFvMesh.C
-solidBodyMotionFunctions = solidBodyMotionFvMesh/solidBodyMotionFunctions
-$(solidBodyMotionFunctions)/solidBodyMotionFunction/solidBodyMotionFunction.C
-$(solidBodyMotionFunctions)/solidBodyMotionFunction/newSolidBodyMotionFunction.C
-$(solidBodyMotionFunctions)/SDA/SDA.C
-$(solidBodyMotionFunctions)/SKA/SKA.C
-$(solidBodyMotionFunctions)/translation/translation.C
-
 mixerGgiFvMesh/mixerGgiFvMesh.C
 turboFvMesh/turboFvMesh.C
 
+fvMeshAdder/fvMeshAdder.C
+fvMeshDistribute/fvMeshDistribute.C
+
 tetMetrics = dynamicTopoFvMesh/tetMetrics
 $(tetMetrics)/tetMetric.C
 $(tetMetrics)/tetMetrics.C

src/dynamicMesh/dynamicFvMesh/Make/options

@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/meshMotion/solidBodyMotion/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/meshMotion/tetDecompositionMotionSolver/lnInclude \
     -I$(LIB_SRC)/tetDecompositionFiniteElement/lnInclude \
     $(WM_DECOMP_INC) \
@@ -15,6 +16,7 @@ LIB_LIBS = \
     -lmeshTools \
     -ldynamicMesh \
     -lfiniteVolume \
+    -lsolidBodyMotion \
     $(WM_DECOMP_LIBS) \
     -lfvMotionSolver \
     -lRBFMotionSolver \

src/dynamicMesh/dynamicFvMesh/fvMeshAdder/fvMeshAdder.H

@@ -65,10 +65,8 @@ class fvMeshAdder
 :
     public polyMeshAdder
 {
-
 private:
 
-
     // Private class
 
     class directFvPatchFieldMapper

src/dynamicMesh/dynamicFvMesh/fvMeshDistribute/fvMeshDistribute.C

@@ -209,7 +209,7 @@ void Foam::fvMeshDistribute::printCoupleInfo
 }
 
 
-// Finds (non-empty) patch that exposed internal and proc faces can be put into.
+// Finds (non-empty) patch that exposed internal and proc faces can be put into
 Foam::label Foam::fvMeshDistribute::findNonEmptyPatch() const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -288,9 +288,11 @@ Foam::label Foam::fvMeshDistribute::addProcPatch
 
     if (polyPatches.findPatchID(patchName) != -1)
     {
-        FatalErrorIn("fvMeshDistribute::addProcPatch(const word&, const label)")
+        FatalErrorIn
-            << "Cannot create patch " << patchName << " since already exists."
+        (
-            << nl
+            "fvMeshDistribute::addProcPatch(const word&, const label)"
+        )   << "Cannot create patch " << patchName
+            << " since it already exists." << nl
             << "Current patch names:" << polyPatches.names()
             << exit(FatalError);
     }
@@ -635,7 +637,7 @@ void Foam::fvMeshDistribute::getNeighbourData
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Get neighbouring meshFace labels and new processor of coupled boundaries.
+    // Get neighbouring meshFace labels and new processor of coupled boundaries
     labelList nbrFaces(nBnd, -1);
     labelList nbrNewProc(nBnd, -1);
 
@@ -931,8 +933,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::doRemoveCells
 }
 
 
-// Delete and add processor patches. Changes mesh and returns per neighbour proc
+// Delete and add processor patches. Changes mesh and returns per
-// the processor patchID.
+// neighbour processor the processor patchID.
 void Foam::fvMeshDistribute::addProcPatches
 (
     const labelList& neighbourNewProc,   // processor that neighbour is on
@@ -1371,7 +1373,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
         FatalErrorIn("fvMeshDistribute::distribute(const labelList&)")
             << "This application requires all non-processor patches"
             << " to be present in the same order on all patches" << nl
-            << "followed by the processor patches (which of course are unique)."
+            << "followed by the processor patches (which are unique)."
             << nl
             << "Local patches:" << mesh_.boundaryMesh().names()
             << abort(FatalError);

src/dynamicMesh/dynamicFvMesh/fvMeshDistribute/fvMeshDistributeTemplates.C

@@ -195,7 +195,8 @@ void Foam::fvMeshDistribute::addPatchFields(const word& patchFieldType)
 
     for
     (
-        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        typename HashTable<const GeoField*>::const_iterator iter =
+            flds.begin();
         iter != flds.end();
         ++iter
     )
@@ -235,7 +236,8 @@ void Foam::fvMeshDistribute::deleteTrailingPatchFields()
 
     for
     (
-        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        typename HashTable<const GeoField*>::const_iterator iter =
+            flds.begin();
         iter != flds.end();
         ++iter
     )
@@ -347,7 +349,11 @@ void Foam::fvMeshDistribute::mapBoundaryFields
                 {
                     label oldLocalI = oldFaceI - oldPatchStarts[oldPatchI];
 
-                    if (oldLocalI >= 0 && oldLocalI < oldBfld[oldPatchI].size())
+                    if
+                    (
+                        oldLocalI >= 0
+                     && oldLocalI < oldBfld[oldPatchI].size()
+                    )
                     {
                         patchFld[i] = oldBfld[oldPatchI][oldLocalI];
                     }
@@ -372,7 +378,8 @@ void Foam::fvMeshDistribute::initPatchFields
 
     for
     (
-        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        typename HashTable<const GeoField*>::const_iterator iter =
+            flds.begin();
         iter != flds.end();
         ++iter
     )
@@ -396,7 +403,7 @@ void Foam::fvMeshDistribute::initPatchFields
 }
 
 
-// Send fields. Note order supplied so we can receive in exactly the same order.
+// Send fields. Note order supplied so we can receive in exactly the same order
 // Note that field gets written as entry in dictionary so we
 // can construct from subdictionary.
 // (since otherwise the reading as-a-dictionary mixes up entries from
@@ -413,7 +420,7 @@ void Foam::fvMeshDistribute::initPatchFields
 //  }
 
 // volVectorField {U {internalField ..; boundaryField ..;}}
-// 
+
 template<class GeoField>
 void Foam::fvMeshDistribute::sendFields
 (
@@ -423,7 +430,8 @@ void Foam::fvMeshDistribute::sendFields
     OSstream& toNbr
 )
 {
-    toNbr << GeoField::typeName << token::NL << token::BEGIN_BLOCK << token::NL;
+    toNbr<< GeoField::typeName << token::NL << token::BEGIN_BLOCK
+        << token::NL;
     forAll(fieldNames, i)
     {
         if (debug)

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/SDA/SDA.C

@@ -1,135 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 "SDA.H"
-#include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
-
-using namespace Foam::mathematicalConstant;
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidBodyMotionFunctions
-{
-    defineTypeNameAndDebug(SDA, 0);
-    addToRunTimeSelectionTable(solidBodyMotionFunction, SDA, dictionary);
-};
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunctions::SDA::SDA
-(
-    const dictionary& SBMFCoeffs,
-    const Time& runTime
-)
-:
-    solidBodyMotionFunction(SBMFCoeffs, runTime),
-    CofG_(SBMFCoeffs_.lookup("CofG"))
-{
-    read(SBMFCoeffs);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunctions::SDA::~SDA()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::septernion Foam::solidBodyMotionFunctions::SDA::transformation() const
-{
-    scalar time = time_.value();
-
-    scalar Tpi = Tp_ + dTp_*(time/dTi_);   // Current roll period [sec]
-    scalar wr = 2*pi/Tpi; // Current Freq [/sec]
-
-    // Current Phase for roll [rad]
-    scalar r = dTp_/dTi_;
-    scalar u = Tp_ + r*time;
-    scalar phr = 2*pi*((Tp_/u - 1) + log(mag(u)) - log(Tp_))/r;
-
-    // Current Phase for Sway [rad]
-    scalar phs = phr + pi;
-
-    // Current Phase for Heave [rad]
-    scalar phh = phr + pi/2;
-
-    scalar rollA = max(rollAmax_*exp(-sqr(Tpi - Tpn_)/(2*Q_)), rollAmin_);
-
-    vector T
-    (
-        0,
-        swayA_*(sin(wr*time + phs) - sin(phs)),
-        heaveA_*(sin(wr*time + phh) - sin(phh))
-    );
-    quaternion R(rollA*sin(wr*time + phr), 0, 0);
-    septernion TR(septernion(CofG_ + T)*R*septernion(-CofG_));
-
-    Info<< "solidBodyMotionFunctions::SDA::transformation(): "
-        << "Time = " << time << " transformation: " << TR << endl;
-
-    return TR;
-}
-
-
-bool Foam::solidBodyMotionFunctions::SDA::read(const dictionary& SBMFCoeffs)
-{
-    solidBodyMotionFunction::read(SBMFCoeffs);
-
-    SBMFCoeffs_.lookup("CofG") >> CofG_;
-    SBMFCoeffs_.lookup("lamda") >> lamda_;
-    SBMFCoeffs_.lookup("rollAmax") >> rollAmax_;
-    SBMFCoeffs_.lookup("rollAmin") >> rollAmin_;
-    SBMFCoeffs_.lookup("heaveA") >> heaveA_;
-    SBMFCoeffs_.lookup("swayA") >> swayA_;
-    SBMFCoeffs_.lookup("Q") >> Q_;
-    SBMFCoeffs_.lookup("Tp") >> Tp_;
-    SBMFCoeffs_.lookup("Tpn") >> Tpn_;
-    SBMFCoeffs_.lookup("dTi") >> dTi_;
-    SBMFCoeffs_.lookup("dTp") >> dTp_;
-
-    // Rescale parameters according to the given scale parameter
-    if (lamda_ > 1 + SMALL)
-    {
-        heaveA_ /= lamda_;
-        swayA_ /= lamda_;
-        Tp_ /= sqrt(lamda_);
-        Tpn_ /= sqrt(lamda_);
-        dTi_ /= sqrt(lamda_);
-        dTp_ /= sqrt(lamda_);
-    }
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/SDA/SDA.H

@@ -1,147 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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::solidBodyMotionFunctions::SDA
-
-Description
-    Ship design analysis (SDA) 3DoF motion function.
-
-    Comprising sinusoidal roll (rotation about x), heave (z-translation)
-    and sway (y-translation) motions with changing amplitude and phase.
-
-See Also
-    SKA (Sea Keeping Analysis) for 6DoF motion.
-
-SourceFiles
-    SDA.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef SDA_H
-#define SDA_H
-
-#include "solidBodyMotionFunction.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidBodyMotionFunctions
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class SDA Declaration
-\*---------------------------------------------------------------------------*/
-
-class SDA
-:
-    public solidBodyMotionFunction
-{
-    // Private data
-
-        //- Center of gravity
-        vector CofG_;
-
-        //- Model scale ratio
-        scalar lamda_;
-
-        //- Max roll amplitude [rad]
-        scalar rollAmax_;
-
-        //- Min roll amplitude [rad]
-        scalar rollAmin_;
-
-        //- Heave amplitude [m]
-        scalar heaveA_;
-
-        //- Sway amplitude [m]
-        scalar swayA_;
-
-        //- Damping Coefficient [-]
-        scalar Q_;
-    
-        //- Time Period for liquid [sec]
-        scalar Tp_;
-
-        //- Natural Period of Ship [sec]
-        scalar Tpn_;
-
-        //- Reference time step [sec]
-        scalar dTi_;
-
-        //- Incr. in Tp/unit 'dTi'[-]
-        scalar dTp_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        SDA(const SDA&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const SDA&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("SDA");
-
-
-    // Constructors
-
-        //- Construct from components
-        SDA
-        (
-            const dictionary& SBMFCoeffs,
-            const Time& runTime
-        );
-
-
-    // Destructor
-
-        virtual ~SDA();
-
-
-    // Member Functions
-
-        //- Return the solid-body motion transformation septernion
-        virtual septernion transformation() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& SBMFCoeffs);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace solidBodyMotionFunctions
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/SKA/SKA.C

@@ -1,162 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 "SKA.H"
-#include "addToRunTimeSelectionTable.H"
-#include "Tuple2.H"
-#include "IFstream.H"
-#include "interpolateXY.H"
-#include "mathematicalConstants.H"
-
-using namespace Foam::mathematicalConstant;
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidBodyMotionFunctions
-{
-    defineTypeNameAndDebug(SKA, 0);
-    addToRunTimeSelectionTable(solidBodyMotionFunction, SKA, dictionary);
-};
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunctions::SKA::SKA
-(
-    const dictionary& SBMFCoeffs,
-    const Time& runTime
-)
-:
-    solidBodyMotionFunction(SBMFCoeffs, runTime)
-{
-    read(SBMFCoeffs);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunctions::SKA::~SKA()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::septernion Foam::solidBodyMotionFunctions::SKA::transformation() const
-{
-    scalar t = time_.value();
-
-    if (t < times_[0])
-    {
-        FatalErrorIn
-        (
-            "solidBodyMotionFunctions::SKA::transformation()"
-        )   << "current time (" << t
-            << ") is less than the minimum in the data table ("
-            << times_[0] << ')'
-            << exit(FatalError);
-    }
-
-    if (t > times_[times_.size()-1])
-    {
-        FatalErrorIn
-        (
-            "solidBodyMotionFunctions::SKA::transformation()"
-        )   << "current time (" << t
-            << ") is greater than the maximum in the data table ("
-            << times_[times_.size()-1] << ')'
-            << exit(FatalError);
-    }
-
-    translationRotationVectors TRV = interpolateXY
-    (
-        t,
-        times_,
-        values_
-    );
-
-    // Convert the rotational motion from deg to rad
-    TRV[1] *= pi/180.0;
-
-    quaternion R(TRV[1].x(), TRV[1].y(), TRV[1].z());
-    septernion TR(septernion(CofG_ + TRV[0])*R*septernion(-CofG_));
-
-    Info<< "solidBodyMotionFunctions::SKA::transformation(): "
-        << "Time = " << t << " transformation: " << TR << endl;
-
-    return TR;
-}
-
-
-bool Foam::solidBodyMotionFunctions::SKA::read(const dictionary& SBMFCoeffs)
-{
-    solidBodyMotionFunction::read(SBMFCoeffs);
-
-    // If the timeDataFileName has changed read the file
-
-    fileName newTimeDataFileName(SBMFCoeffs_.lookup("timeDataFileName"));
-
-    if (newTimeDataFileName != timeDataFileName_)
-    {
-        timeDataFileName_ = newTimeDataFileName;
-
-        IFstream dataStream(timeDataFileName_);
-
-        if (dataStream.good())
-        {
-            List<Tuple2<scalar, translationRotationVectors> > timeValues
-            (
-                dataStream
-            );
-
-            times_.setSize(timeValues.size());
-            values_.setSize(timeValues.size());
-
-            forAll(timeValues, i)
-            {
-                times_[i] = timeValues[i].first();
-                values_[i] = timeValues[i].second();
-            }
-        }
-        else
-        {
-            FatalErrorIn
-            (
-                "solidBodyMotionFunctions::SKA::read(const dictionary&)"
-            )   << "Cannot open time data file " << timeDataFileName_
-                << exit(FatalError);
-        }
-    }
-
-    SBMFCoeffs_.lookup("CofG") >> CofG_;
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/SKA/SKA.H

@@ -1,132 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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::solidBodyMotionFunctions::SKA
-
-Description
-    Sea Keeping Analysis (SKA) 6DoF motion function.
-
-    Obtained by interpolating tabulated data for surge (x-translation),
-    sway (y-translation), heave (z-translation), roll (rotation about x),
-    pitch (rotation about y) and yaw (rotation about z).
-
-See Also
-    SDA (Ship design analysis) for 3DoF motion.
-
-SourceFiles
-    SKA.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef SKA_H
-#define SKA_H
-
-#include "solidBodyMotionFunction.H"
-#include "primitiveFields.H"
-#include "Vector2D.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidBodyMotionFunctions
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class SKA Declaration
-\*---------------------------------------------------------------------------*/
-
-class SKA
-:
-    public solidBodyMotionFunction
-{
-    // Private data
-
-        //- Time data file name read from dictionary
-        fileName timeDataFileName_;
-
-        //- Center of gravity read from dictionary
-        vector CofG_;
-
-        //- Type used to read in the translation and rotation "vectors"
-        typedef Vector2D<vector> translationRotationVectors;
-
-        //- Field of times
-        scalarField times_;
-
-        //- Field of translation and rotation "vectors"
-        Field<translationRotationVectors> values_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        SKA(const SKA&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const SKA&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("SKA");
-
-
-    // Constructors
-
-        //- Construct from components
-        SKA
-        (
-            const dictionary& SBMFCoeffs,
-            const Time& runTime
-        );
-
-
-    // Destructor
-
-        virtual ~SKA();
-
-
-    // Member Functions
-
-        //- Return the solid-body motion transformation septernion
-        virtual septernion transformation() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& SBMFCoeffs);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace solidBodyMotionFunctions
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/solidBodyMotionFunction/newSolidBodyMotionFunction.C

@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 "solidBodyMotionFunction.H"
-
-// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::solidBodyMotionFunction> Foam::solidBodyMotionFunction::New
-(
-    const dictionary& SBMFCoeffs,
-    const Time& runTime
-)
-{
-    word solidBodyMotionFunctionTypeName = 
-        SBMFCoeffs.lookup("solidBodyMotionFunction");
-
-    Info<< "Selecting solid-body motion function "
-        << solidBodyMotionFunctionTypeName << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(solidBodyMotionFunctionTypeName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalErrorIn
-        (
-            "solidBodyMotionFunction::New"
-            "("
-            "    const dictionary& SBMFCoeffs,"
-            "    const Time& runTime"
-            ")"
-        )   << "Unknown solidBodyMotionFunction type "
-            << solidBodyMotionFunctionTypeName << endl << endl
-            << "Valid  solidBodyMotionFunctions are : " << endl
-            << dictionaryConstructorTablePtr_->toc()
-            << exit(FatalError);
-    }
-
-    return autoPtr<solidBodyMotionFunction>(cstrIter()(SBMFCoeffs, runTime));
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/solidBodyMotionFunction/solidBodyMotionFunction.C

@@ -1,70 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 "solidBodyMotionFunction.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(Foam::solidBodyMotionFunction, 0);
-
-defineRunTimeSelectionTable(Foam::solidBodyMotionFunction, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunction::solidBodyMotionFunction
-(
-    const dictionary& SBMFCoeffs,
-    const Time& runTime
-)
-:
-    SBMFCoeffs_
-    (
-        SBMFCoeffs.subDict
-        (
-            word(SBMFCoeffs.lookup("solidBodyMotionFunction")) + "Coeffs"
-        )
-    ),
-    time_(runTime)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunction::~solidBodyMotionFunction()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-bool Foam::solidBodyMotionFunction::read(const dictionary& SBMFCoeffs)
-{
-    SBMFCoeffs_ = SBMFCoeffs.subDict(type() + "Coeffs");
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/translation/translation.C

@@ -1,98 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 "translation.H"
-#include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
-
-using namespace Foam::mathematicalConstant;
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidBodyMotionFunctions
-{
-    defineTypeNameAndDebug(translation, 0);
-    addToRunTimeSelectionTable
-    (
-        solidBodyMotionFunction,
-        translation,
-        dictionary
-    );
-};
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunctions::translation::translation
-(
-    const dictionary& SBMFCoeffs,
-    const Time& runTime
-)
-:
-    solidBodyMotionFunction(SBMFCoeffs, runTime),
-    velocity_(SBMFCoeffs_.lookup("velocity"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::solidBodyMotionFunctions::translation::~translation()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::septernion
-Foam::solidBodyMotionFunctions::translation::transformation() const
-{
-    scalar time = time_.value();
-
-    septernion TR(velocity_*time, quaternion::I);
-
-    Info<< "solidBodyMotionFunctions::translation::transformation(): "
-        << "Time = " << time << " transformation: " << TR << endl;
-
-    return TR;
-}
-
-
-bool Foam::solidBodyMotionFunctions::translation::read
-(
-    const dictionary& SBMFCoeffs
-)
-{
-    solidBodyMotionFunction::read(SBMFCoeffs);
-
-    SBMFCoeffs_.lookup("velocity") >> velocity_;
-
-    return true;
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/translation/translation.H

@@ -1,111 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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::solidBodyMotionFunctions::translation
-
-Description
-    translation motion function
-
-SourceFiles
-    translation.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef translation_H
-#define translation_H
-
-#include "solidBodyMotionFunction.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace solidBodyMotionFunctions
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class translation Declaration
-\*---------------------------------------------------------------------------*/
-
-class translation
-:
-    public solidBodyMotionFunction
-{
-    // Private data
-
-        //- Velocity
-        vector velocity_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        translation(const translation&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const translation&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("translation");
-
-
-    // Constructors
-
-        //- Construct from components
-        translation
-        (
-            const dictionary& SBMFCoeffs,
-            const Time& runTime
-        );
-
-
-    // Destructor
-
-        virtual ~translation();
-
-
-    // Member Functions
-
-        //- Return the solid-body motion transformation septernion
-        virtual septernion transformation() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& SBMFCoeffs);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace solidBodyMotionFunctions
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFvMesh.C

@@ -53,20 +53,20 @@ Foam::solidBodyMotionFvMesh::solidBodyMotionFvMesh(const IOobject& io)
             IOobject
             (
                 "dynamicMeshDict",
-                io.time().constant(),
+                time().constant(),
                 *this,
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE
             )
         ).subDict(typeName + "Coeffs")
     ),
-    SBMFPtr_(solidBodyMotionFunction::New(dynamicMeshCoeffs_, io.time())),
+    SBMFPtr_(solidBodyMotionFunction::New(dynamicMeshCoeffs_, time())),
     undisplacedPoints_
     (
         IOobject
         (
             "points",
-            io.time().constant(),
+            time().constant(),
             meshSubDir,
             *this,
             IOobject::MUST_READ,
@@ -88,8 +88,7 @@ bool Foam::solidBodyMotionFvMesh::update()
 {
     fvMesh::movePoints
     (
-        transform(SBMFPtr_().transformation(),
+        transform(SBMFPtr_().transformation(), undisplacedPoints_)
-        undisplacedPoints_)
     );
 
     return false;

src/dynamicMesh/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFvMesh.H

@@ -60,10 +60,10 @@ class solidBodyMotionFvMesh
         //- Dictionary of motion control parameters
         dictionary dynamicMeshCoeffs_;
 
-        //- The motion control function
+        //- Motion control function
         autoPtr<solidBodyMotionFunction> SBMFPtr_;
 
-        //- The reference points which are transformed
+        //- Reference points which are transformed
         pointIOField undisplacedPoints_;
 
 
@@ -85,7 +85,7 @@ public:
     // Constructors
 
         //- Construct from IOobject
-        solidBodyMotionFvMesh(const IOobject& io);
+        explicit solidBodyMotionFvMesh(const IOobject& io);
 
 
     // Destructor

src/dynamicMesh/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.C

@@ -317,25 +317,10 @@ bool Foam::layerAdditionRemoval::changeTopology() const
     // If the patch is empty on a processor in a parallel simulation,
     // original values will be preserved.  HJ, 7/Mar/2011
 
-    Pout<< "HRVOJE 1: " << minDelta;
+    reduce(minDelta, minOp<scalar>());
-//     reduce(minDelta, minOp<scalar>());
+    reduce(maxDelta, maxOp<scalar>());
-    Foam::sleep(2);
+    reduce(avgDelta, sumOp<scalar>());
-    Pout<< " and " << minDelta << endl;;
+    reduce(nAvg, sumOp<scalar>());
-
-    Pout<< "2: " << maxDelta;
-//     reduce(maxDelta, maxOp<scalar>());
-    Foam::sleep(2);
-    Pout<< " and " << maxDelta << endl;;
-
-    Pout<< "3: " << avgDelta;
-//     reduce(avgDelta, sumOp<scalar>());
-    Foam::sleep(2);
-    Pout<< " and " << avgDelta << endl;;
-
-    Pout<< "4: " << nAvg;
-//     reduce(nAvg, sumOp<scalar>());
-    Foam::sleep(2);
-    Pout<< " and " << nAvg << endl;;
 
     avgDelta /= nAvg;
 

src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C

@@ -1014,8 +1014,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
                 )   << "Face " << faceI << " in the new mesh is not "
                     << "mapped correctly." << nl
                     << "It uses a removed or a non-existing vertex or "
-                    << "has been skipped ." << nl
+                    << "has been skipped." << nl
-                    << "Face before mapping: " << oldFace << nl
+                    << "Face before mapping: " << oldFace << " with points "
+                    << oldFace.points(newPointsZeroVol) << nl
+                    << mesh.allPoints().size() << nl
                     << "Face after mapping: " << renumberedFace << nl
                     << "Max new vertex index: "
                     << newPointsZeroVol.size() - 1 << "." << nl
@@ -1098,6 +1100,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh
 
             if (rotate != 0)
             {
+                Info<< "Rotating face" << endl;
                 newFaces[faceI] = rotateFace(newFaces[faceI], rotate);
             }
         }

src/dynamicMesh/topoChangerFvMesh/Make/files

@@ -3,7 +3,7 @@ topoChangerFvMesh/topoChangerFvMesh.C
 linearValveFvMesh/linearValveFvMesh.C
 attachDetachFvMesh/attachDetachFvMesh.C
 linearValveLayersFvMesh/linearValveLayersFvMesh.C
-movingConeTopoFvMesh/movingConeTopoFvMesh.C
+movingBodyTopoFvMesh/movingBodyTopoFvMesh.C
 mixerFvMesh/mixerFvMesh.C
 multiMixerFvMesh/mixerRotor.C
 multiMixerFvMesh/multiMixerFvMesh.C

src/dynamicMesh/topoChangerFvMesh/Make/options

@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/meshMotion/solidBodyMotion/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/meshMotion/tetDecompositionMotionSolver/lnInclude \
     -I$(LIB_SRC)/tetDecompositionFiniteElement/lnInclude \
     $(WM_DECOMP_INC)
@@ -12,4 +13,5 @@ LIB_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lmeshTools \
+    -lsolidBodyMotion \
     $(WM_DECOMP_LIBS)

src/dynamicMesh/topoChangerFvMesh/mixerFvMesh/mixerFvMesh.C

@@ -46,50 +46,40 @@ namespace Foam
 void Foam::mixerFvMesh::addZonesAndModifiers()
 {
     // Add zones and modifiers for motion action
+    Pout<< "pointZones().size(): " << pointZones().size()
+        << " faceZones().size(): " << faceZones().size()
+        << " cellZones().size(): " << cellZones().size()
+        << " topoChanger_.size(): " << topoChanger_.size() << endl;
 
-    if
+    if (cellZones().size() > 0)
-    (
-        pointZones().size() > 0
-     || faceZones().size() > 0
-     || cellZones().size() > 0
-    )
     {
         Info<< "void mixerFvMesh::addZonesAndModifiers() : "
             << "Zones and modifiers already present.  Skipping."
             << endl;
 
-        if (topoChanger_.size() == 0)
+        // Check definition of the modifier
+        if
+        (
+            pointZones().size() > 0
+         || faceZones().size() > 0
+        )
         {
-            FatalErrorIn
+            if (topoChanger_.size() == 0)
-            (
+            {
-                "void mixerFvMesh::addZonesAndModifiers()"
+                FatalErrorIn
-            )   << "Mesh modifiers not read properly"
+                (
-                << abort(FatalError);
+                    "void mixerFvMesh::addZonesAndModifiers()"
+                )   << "Mesh modifiers not read properly.  "
+                    << "pointZones = " <<  pointZones().size()
+                    << " faceZones = " << faceZones().size()
+                    << abort(FatalError);
+            }
         }
 
         return;
     }
 
-    Info<< "Time = " << time().timeName() << endl
+    // Find patches and check sizes
-        << "Adding zones and modifiers to the mesh" << endl;
-
-    // Add zones
-    List<pointZone*> pz(1);
-
-    // Add an empty zone for cut points
-
-    pz[0] = new pointZone
-    (
-        "cutPointZone",
-        labelList(0),
-        0,
-        pointZones()
-    );
-
-
-    // Do face zones for slider
-
-    List<faceZone*> fz(3);
 
     // Moving slider
     const word movingSliderName(dict_.subDict("slider").lookup("moving"));
@@ -114,51 +104,84 @@ void Foam::mixerFvMesh::addZonesAndModifiers()
     }
 
     const polyPatch& movingSlider = boundaryMesh()[movingSliderIndex];
-
-    labelList isf(movingSlider.size());
-
-    forAll (isf, i)
-    {
-        isf[i] = movingSlider.start() + i;
-    }
-
-    fz[0] = new faceZone
-    (
-        movingSliderName + "Zone",
-        isf,
-        boolList(movingSlider.size(), false),
-        0,
-        faceZones()
-    );
-
-    // Static slider
     const polyPatch& staticSlider = boundaryMesh()[staticSliderIndex];
 
-    labelList osf(staticSlider.size());
+    List<pointZone*> pz;
+    List<faceZone*> fz;
 
-    forAll (osf, i)
+    bool addSlider = false;
+
+    if (movingSlider.size() > 0 && staticSlider.size() > 0)
     {
-        osf[i] = staticSlider.start() + i;
+        addSlider = true;
+
+        pz.setSize(1);
+        fz.setSize(3);
+
+        Info<< "Time = " << time().timeName() << endl
+            << "Adding zones and modifiers to the mesh" << endl;
+
+        // Add zones
+
+        // Add an empty zone for cut points
+
+        pz[0] = new pointZone
+        (
+            "cutPointZone",
+            labelList(0),
+            0,
+            pointZones()
+        );
+
+
+        // Do face zones for slider
+
+        // Moving slider
+        labelList isf(movingSlider.size());
+
+        forAll (isf, i)
+        {
+            isf[i] = movingSlider.start() + i;
+        }
+
+        fz[0] = new faceZone
+        (
+            movingSliderName + "Zone",
+            isf,
+            boolList(movingSlider.size(), false),
+            0,
+            faceZones()
+        );
+
+        // Static slider
+        labelList osf(staticSlider.size());
+
+        forAll (osf, i)
+        {
+            osf[i] = staticSlider.start() + i;
+        }
+
+        fz[1] = new faceZone
+        (
+            staticSliderName + "Zone",
+            osf,
+            boolList(staticSlider.size(), false),
+            1,
+            faceZones()
+        );
+
+        // Add empty zone for cut faces
+        fz[2] = new faceZone
+        (
+            "cutFaceZone",
+            labelList(0),
+            boolList(0, false),
+            2,
+            faceZones()
+        );
     }
 
-    fz[1] = new faceZone
+    // Add cell zone even if slider does not exist
-    (
-        staticSliderName + "Zone",
-        osf,
-        boolList(staticSlider.size(), false),
-        1,
-        faceZones()
-    );
-
-    // Add empty zone for cut faces
-    fz[2] = new faceZone
-    (
-        "cutFaceZone",
-        labelList(0),
-        boolList(0, false),
-        2,
-        faceZones()
-    );
 
     List<cellZone*> cz(1);
 
@@ -194,28 +217,32 @@ void Foam::mixerFvMesh::addZonesAndModifiers()
     Info << "Adding point, face and cell zones" << endl;
     addZones(pz, fz, cz);
 
-    // Add a topology modifier
+    if (addSlider)
-    Info << "Adding topology modifiers" << endl;
+    {
-    topoChanger_.setSize(1);
+        // Add a topology modifier
-    topoChanger_.set
+        Info << "Adding topology modifiers" << endl;
-    (
+        topoChanger_.setSize(1);
-        0,
+
-        new slidingInterface
+        topoChanger_.set
         (
-            "mixerSlider",
             0,
-            topoChanger_,
+            new slidingInterface
-            staticSliderName + "Zone",
+            (
-            movingSliderName + "Zone",
+                "mixerSlider",
-            "cutPointZone",
+                0,
-            "cutFaceZone",
+                topoChanger_,
-            staticSliderName,
+                staticSliderName + "Zone",
-            movingSliderName,
+                movingSliderName + "Zone",
-            slidingInterface::INTEGRAL,   // Edge matching algorithm
+                "cutPointZone",
-            attachDetach_,                // Attach-detach action
+                "cutFaceZone",
-            intersection::VISIBLE         // Projection algorithm
+                staticSliderName,
-        )
+                movingSliderName,
-    );
+                slidingInterface::INTEGRAL,   // Edge matching algorithm
+                attachDetach_,                // Attach-detach action
+                intersection::VISIBLE         // Projection algorithm
+            )
+        );
+    }
 
     // Write mesh and modifiers
     topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
@@ -226,7 +253,16 @@ void Foam::mixerFvMesh::addZonesAndModifiers()
 
 bool Foam::mixerFvMesh::attached() const
 {
-    return refCast<const slidingInterface>(topoChanger_[0]).attached();
+    bool result = false;
+
+    if (topoChanger_.size() > 0)
+    {
+        result = refCast<const slidingInterface>(topoChanger_[0]).attached();
+    }
+
+    reduce(result, orOp<bool>());
+
+    return result;
 }
 
 
@@ -272,41 +308,43 @@ void Foam::mixerFvMesh::calcMovingMask() const
         }
     }
 
-    const word movingSliderZoneName
+    if (topoChanger_.size() > 0)
-    (
-        word(dict_.subDict("slider").lookup("moving"))
-      + "Zone"
-    );
-
-    const labelList& movingSliderAddr =
-        faceZones()[faceZones().findZoneID(movingSliderZoneName)];
-
-    forAll (movingSliderAddr, faceI)
     {
-        const face& curFace = f[movingSliderAddr[faceI]];
+        // Topo changer present.  Use zones for motion
+        const word movingSliderZoneName
+        (
+            word(dict_.subDict("slider").lookup("moving")) + "Zone"
+        );
 
-        forAll (curFace, pointI)
+        const labelList& movingSliderAddr =
+            faceZones()[faceZones().findZoneID(movingSliderZoneName)];
+
+        forAll (movingSliderAddr, faceI)
         {
-            movingPointsMask[curFace[pointI]] = 1;
+            const face& curFace = f[movingSliderAddr[faceI]];
+
+            forAll (curFace, pointI)
+            {
+                movingPointsMask[curFace[pointI]] = 1;
+            }
         }
-    }
 
-    const word staticSliderZoneName
+        const word staticSliderZoneName
-    (
+        (
-        word(dict_.subDict("slider").lookup("static"))
+            word(dict_.subDict("slider").lookup("static")) + "Zone"
-      + "Zone"
+        );
-    );
 
-    const labelList& staticSliderAddr =
+        const labelList& staticSliderAddr =
-        faceZones()[faceZones().findZoneID(staticSliderZoneName)];
+            faceZones()[faceZones().findZoneID(staticSliderZoneName)];
 
-    forAll (staticSliderAddr, faceI)
+        forAll (staticSliderAddr, faceI)
-    {
-        const face& curFace = f[staticSliderAddr[faceI]];
-
-        forAll (curFace, pointI)
         {
-            movingPointsMask[curFace[pointI]] = 0;
+            const face& curFace = f[staticSliderAddr[faceI]];
+
+            forAll (curFace, pointI)
+            {
+                movingPointsMask[curFace[pointI]] = 0;
+            }
         }
     }
 }
@@ -415,22 +453,43 @@ bool Foam::mixerFvMesh::update()
 
     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 
-    if (topoChangeMap->morphing())
+    bool localMorphing = topoChangeMap->morphing();
+    bool globalMorphing = localMorphing;
+
+    reduce(globalMorphing, orOp<bool>());
+
+    if (globalMorphing)
     {
         Info << "Attaching rotors" << endl;
 
         deleteDemandDrivenData(movingPointsMaskPtr_);
 
         // Move the sliding interface points to correct position
-        pointField mappedOldPointsNew(allPoints().size());
+        if (localMorphing)
-        mappedOldPointsNew.map(oldPointsNew, topoChangeMap->pointMap());
+        {
+            pointField mappedOldPointsNew(allPoints().size());
+            mappedOldPointsNew.map(oldPointsNew, topoChangeMap->pointMap());
 
-        movePoints(mappedOldPointsNew);
+            movePoints(mappedOldPointsNew);
-        resetMotion();
+
-        setV0();
+            resetMotion();
+            setV0();
+
+            // Move the sliding interface points to correct position
+            movePoints(topoChangeMap->preMotionPoints());
+        }
+        else
+        {
+            pointField newPoints = allPoints();
+            movePoints(oldPointsNew);
+
+            resetMotion();
+            setV0();
+
+            // Move the sliding interface points to correct position
+            movePoints(newPoints);
+        }
 
-        // Move the sliding interface points to correct position
-        movePoints(topoChangeMap->preMotionPoints());
     }
 
     return topoChangeMap->morphing();

src/dynamicMesh/topoChangerFvMesh/movingBodyTopoFvMesh/movingBodyTopoFvMesh.C

@@ -0,0 +1,306 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 "movingBodyTopoFvMesh.H"
+#include "Time.H"
+#include "mapPolyMesh.H"
+#include "layerAdditionRemoval.H"
+#include "volMesh.H"
+#include "transformField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(movingBodyTopoFvMesh, 0);
+
+    addToRunTimeSelectionTable
+    (
+        topoChangerFvMesh,
+        movingBodyTopoFvMesh,
+        IOobject
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField>
+Foam::movingBodyTopoFvMesh::calcMotionMask() const
+{
+    Info<< "Updating vertex markup" << endl;
+
+    tmp<scalarField> tvertexMarkup(new scalarField(allPoints().size(), 0));
+    scalarField& vertexMarkup = tvertexMarkup();
+
+    cellZoneID movingCellsID(movingCellsName_, cellZones());
+
+    // In order to do a correct update on a mask on processor boundaries,
+    // Detection of moving cells should use patchNeighbourField for
+    // processor (not coupled!) boundaries.  This is done by expanding
+    // a moving cell set into a field and making sure that processor patch
+    // points move in sync.  Not done at the moment, probably best to do
+    // using parallel update of pointFields.  HJ, 19/Feb/2011
+
+    // If moving cells are found, perform mark-up
+    if (movingCellsID.active())
+    {
+        // Get cell-point addressing
+        const labelListList& cp = cellPoints();
+
+        // Get labels of all moving cells
+        const labelList& movingCells = cellZones()[movingCellsID.index()];
+
+        forAll (movingCells, cellI)
+        {
+            const labelList& curCp = cp[movingCells[cellI]];
+
+            forAll (curCp, pointI)
+            {
+                vertexMarkup[curCp[pointI]] = 1;
+            }
+        }
+    }
+
+    faceZoneID frontFacesID(frontFacesName_, faceZones());
+
+    if (frontFacesID.active())
+    {
+        const faceZone& frontFaces = faceZones()[frontFacesID.index()];
+
+        const labelList& mp = frontFaces().meshPoints();
+
+        forAll (mp, mpI)
+        {
+            vertexMarkup[mp[mpI]] = 1;
+        }
+    }
+
+    faceZoneID backFacesID(backFacesName_, faceZones());
+
+    if (backFacesID.active())
+    {
+        const faceZone& backFaces = faceZones()[backFacesID.index()];
+
+        const labelList& mp = backFaces().meshPoints();
+
+        forAll (mp, mpI)
+        {
+            vertexMarkup[mp[mpI]] = 1;
+        }
+    }
+
+    return tvertexMarkup;
+}
+
+
+void Foam::movingBodyTopoFvMesh::addZonesAndModifiers()
+{
+    // Add zones and modifiers for motion action
+
+    if (topoChanger_.size() > 0)
+    {
+        Info<< "void movingBodyTopoFvMesh::addZonesAndModifiers() : "
+            << "Zones and modifiers already present.  Skipping."
+            << endl;
+
+        return;
+    }
+
+    // Add layer addition/removal interfaces
+    topoChanger_.setSize(2);
+    label nMods = 0;
+
+
+    faceZoneID frontFacesID(frontFacesName_, faceZones());
+    faceZoneID backFacesID(backFacesName_, faceZones());
+
+    if (frontFacesID.active())
+    {
+        const faceZone& frontFaces = faceZones()[frontFacesID.index()];
+
+        if (!frontFaces.empty())
+        {
+            topoChanger_.set
+            (
+                nMods,
+                new layerAdditionRemoval
+                (
+                    frontFacesName_ + "Layer",
+                    nMods,
+                    topoChanger_,
+                    frontFacesName_,
+                    readScalar
+                    (
+                        dict_.subDict("front").lookup("minThickness")
+                    ),
+                    readScalar
+                    (
+                        dict_.subDict("front").lookup("maxThickness")
+                    )
+                )
+            );
+
+            nMods++;
+        }
+    }
+
+    if (backFacesID.active())
+    {
+        const faceZone& backFaces = faceZones()[backFacesID.index()];
+
+        if (!backFaces.empty())
+        {
+            topoChanger_.set
+            (
+                nMods,
+                new layerAdditionRemoval
+                (
+                    backFacesName_ + "Layer",
+                    nMods,
+                    topoChanger_,
+                    backFacesName_,
+                    readScalar
+                    (
+                        dict_.subDict("back").lookup("minThickness")
+                    ),
+                    readScalar
+                    (
+                        dict_.subDict("back").lookup("maxThickness")
+                    )
+                )
+            );
+
+            nMods++;
+        }
+    }
+
+    topoChanger_.setSize(nMods);
+
+    reduce(nMods, sumOp<label>());
+
+    Info << "Adding " << nMods << " mesh modifiers" << endl;
+
+    // Write mesh and modifiers
+    topoChanger_.write();
+
+    // No need to write the mesh - only modifiers are added.
+    // HJ, 18/Feb/2011
+//     write();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::movingBodyTopoFvMesh::movingBodyTopoFvMesh(const IOobject& io)
+:
+    topoChangerFvMesh(io),
+    dict_
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+        ).subDict(typeName + "Coeffs")
+    ),
+    movingCellsName_(dict_.lookup("movingCells")),
+    frontFacesName_(dict_.lookup("frontFaces")),
+    backFacesName_(dict_.lookup("backFaces")),
+    SBMFPtr_(solidBodyMotionFunction::New(dict_, time())),
+    motionMask_()
+{
+    addZonesAndModifiers();
+    motionMask_ = calcMotionMask();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::movingBodyTopoFvMesh::~movingBodyTopoFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::movingBodyTopoFvMesh::update()
+{
+    // Store points to recreate mesh motion
+    pointField oldPointsNew = allPoints();
+    pointField newPoints = allPoints();
+
+    autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
+
+    bool localMeshChanged = topoChangeMap->morphing();
+    bool globalMeshChanged = localMeshChanged;
+    reduce(globalMeshChanged, orOp<bool>());
+
+    if (globalMeshChanged)
+    {
+        Pout<< "Topology change. Calculating motion point mask" << endl;
+        motionMask_ = calcMotionMask();
+    }
+
+    if (localMeshChanged)
+    {
+//         // Map old points onto the new mesh
+//         pointField mappedOldPointsNew(allPoints().size());
+//         mappedOldPointsNew.map(oldPointsNew, topoChangeMap->pointMap());
+
+//         movePoints(mappedOldPointsNew);
+//         resetMotion();
+//         setV0();
+
+        // Get new points from preMotion
+        newPoints = topoChangeMap().preMotionPoints();
+    }
+//     else
+//     {
+//         // No change, use old points
+//         movePoints(oldPointsNew);
+//         resetMotion();
+//         setV0();
+//     }
+
+    // Calculate new points using a velocity transformation
+    newPoints += motionMask_*
+        transform(SBMFPtr_().velocity(), newPoints)*time().deltaT().value();
+
+    Info << "Executing mesh motion" << endl;
+    movePoints(newPoints);
+
+    return localMeshChanged;
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/movingBodyTopoFvMesh/movingBodyTopoFvMesh.H

@@ -23,22 +23,30 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::movingConeTopoFvMesh
+    Foam::movingBodyTopoFvMesh
 
 Description
-    Sample topoChangerFvMesh that moves an object in x direction
+    Sample topoChangerFvMesh that moves a cell zone in a given direction,
-    and introduces/removes layers.
+    introducing/removing layers around it.
+
+    Cells in the movingCells zone shall be moved given the prescribed velocity
+    and will be bounded in "front" and "back" by other cell zones.
+    Layer addition/removal interfaces are inserted at boundaries between the
+    moving zone and front and back, pointing outside of the moving cell zone
+
+Author and rewrite
+    Hrvoje Jasak, Wikki Ltd.
 
 SourceFiles
-    movingConeTopoFvMesh.C
+    movingBodyTopoFvMesh.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef movingConeTopoFvMesh_H
+#ifndef movingBodyTopoFvMesh_H
-#define movingConeTopoFvMesh_H
+#define movingBodyTopoFvMesh_H
 
 #include "topoChangerFvMesh.H"
-#include "motionSolver.H"
+#include "solidBodyMotionFunction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,35 +54,29 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class movingConeTopoFvMesh Declaration
+                      Class movingBodyTopoFvMesh Declaration
 \*---------------------------------------------------------------------------*/
 
-class movingConeTopoFvMesh
+class movingBodyTopoFvMesh
 :
     public topoChangerFvMesh
 {
     // Private data
 
-        //- Motion dictionary
+        //- Dynamic mesh dictionary
-        dictionary motionDict_;
+        dictionary dict_;
 
-        //- Motion velocity amplitude
+        //- Moving cell zone name
-        vector motionVelAmplitude_;
+        word movingCellsName_;
 
-        //- Motion velocity period
+        //- Front face zone name
-        scalar motionVelPeriod_;
+        word frontFacesName_;
 
-        //- Motion velocity period
+        //- Back face zone name
-        vector curMotionVel_;
+        word backFacesName_;
 
-        //- Left edge
+        //- Motion control function
-        scalar leftEdge_;
+        autoPtr<solidBodyMotionFunction> SBMFPtr_;
-
-        //- Current left obstacle position
-        scalar curLeft_;
-
-        //- Current right obstacle position
-        scalar curRight_;
 
         //- Vertex motion mask
         scalarField motionMask_;
@@ -83,39 +85,34 @@ class movingConeTopoFvMesh
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        movingConeTopoFvMesh(const movingConeTopoFvMesh&);
+        movingBodyTopoFvMesh(const movingBodyTopoFvMesh&);
 
         //- Disallow default bitwise assignment
-        void operator=(const movingConeTopoFvMesh&);
+        void operator=(const movingBodyTopoFvMesh&);
 
 
         //- Add mixer zones and modifiers
         void addZonesAndModifiers();
 
-        //- Markup motion vertices
+        //- Mark motion vertices
-        tmp<scalarField> vertexMarkup
+        tmp<scalarField> calcMotionMask() const;
-        (
-            const pointField& p,
-            const scalar& curLeft,
-            const scalar& curRight
-        ) const;
 
 
 public:
 
     //- Runtime type information
-    TypeName("movingConeTopoFvMesh");
+    TypeName("movingBodyTopoFvMesh");
 
 
     // Constructors
 
         //- Construct from database
-        explicit movingConeTopoFvMesh(const IOobject& io);
+        explicit movingBodyTopoFvMesh(const IOobject& io);
 
 
     // Destructor
 
-        virtual ~movingConeTopoFvMesh();
+        virtual ~movingBodyTopoFvMesh();
 
 
     // Member Functions

src/dynamicMesh/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C

@@ -1,429 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 "movingConeTopoFvMesh.H"
-#include "Time.H"
-#include "mapPolyMesh.H"
-#include "layerAdditionRemoval.H"
-#include "addToRunTimeSelectionTable.H"
-#include "volMesh.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(movingConeTopoFvMesh, 0);
-
-    addToRunTimeSelectionTable
-    (
-        topoChangerFvMesh,
-        movingConeTopoFvMesh,
-        IOobject
-    );
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
-(
-    const pointField& p,
-    const scalar& curLeft,
-    const scalar& curRight
-) const
-{
-    Info<< "Updating vertex markup.  curLeft: "
-        << curLeft << " curRight: " << curRight << endl;
-
-    tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
-    scalarField& vertexMarkup = tvertexMarkup();
-
-    forAll (p, pI)
-    {
-        if (p[pI].x() < curLeft - 1e-10)
-        {
-            vertexMarkup[pI] = -1;
-        }
-        else if (p[pI].x() > curRight + 1e-10)
-        {
-            vertexMarkup[pI] = 1;
-        }
-        else
-        {
-            vertexMarkup[pI] = 0;
-        }
-    }
-
-    return tvertexMarkup;
-}
-
-
-void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
-{
-    // Add zones and modifiers for motion action
-
-    if
-    (
-        pointZones().size() > 0
-     || faceZones().size() > 0
-     || cellZones().size() > 0
-    )
-    {
-        Info<< "void movingConeTopoFvMesh::addZonesAndModifiers() : "
-            << "Zones and modifiers already present.  Skipping."
-            << endl;
-
-        if (topoChanger_.size() == 0)
-        {
-            FatalErrorIn
-            (
-                "void movingConeTopoFvMesh::addZonesAndModifiers()"
-            )   << "Mesh modifiers not read properly"
-                << abort(FatalError);
-        }
-
-        return;
-    }
-
-    Info<< "Time = " << time().timeName() << endl
-        << "Adding zones and modifiers to the mesh" << endl;
-
-    const vectorField& fc = faceCentres();
-    const vectorField& fa = faceAreas();
-
-    labelList zone1(fc.size());
-    boolList flipZone1(fc.size(), false);
-    label nZoneFaces1 = 0;
-
-    labelList zone2(fc.size());
-    boolList flipZone2(fc.size(), false);
-    label nZoneFaces2 = 0;
-
-    forAll (fc, faceI)
-    {
-        if
-        (
-            fc[faceI].x() > -0.003501
-         && fc[faceI].x() < -0.003499
-        )
-        {
-            if ((fa[faceI] & vector(1, 0, 0)) < 0)
-            {
-                flipZone1[nZoneFaces1] = true;
-            }
-
-            zone1[nZoneFaces1] = faceI;
-            Info<< "face " << faceI << " for zone 1.  Flip: "
-                << flipZone1[nZoneFaces1] << endl;
-            nZoneFaces1++;
-        }
-        else if
-        (
-            fc[faceI].x() > -0.00701
-         && fc[faceI].x() < -0.00699
-        )
-        {
-            zone2[nZoneFaces2] = faceI;
-
-            if ((fa[faceI] & vector(1, 0, 0)) > 0)
-            {
-                flipZone2[nZoneFaces2] = true;
-            }
-
-            Info << "face " << faceI << " for zone 2.  Flip: "
-                << flipZone2[nZoneFaces2] << endl;
-            nZoneFaces2++;
-        }
-    }
-
-    zone1.setSize(nZoneFaces1);
-    flipZone1.setSize(nZoneFaces1);
-
-    zone2.setSize(nZoneFaces2);
-    flipZone2.setSize(nZoneFaces2);
-
-    Info << "zone: " << zone1 << endl;
-    Info << "zone: " << zone2 << endl;
-
-    List<pointZone*> pz(0);
-    List<faceZone*> fz(2);
-    List<cellZone*> cz(0);
-
-    label nFz = 0;
-
-    fz[nFz] =
-        new faceZone
-        (
-            "rightExtrusionFaces",
-            zone1,
-            flipZone1,
-            nFz,
-            faceZones()
-        );
-    nFz++;
-
-    fz[nFz] =
-        new faceZone
-        (
-            "leftExtrusionFaces",
-            zone2,
-            flipZone2,
-            nFz,
-            faceZones()
-        );
-    nFz++;
-
-    fz.setSize(nFz);
-
-    Info << "Adding mesh zones." << endl;
-    addZones(pz, fz, cz);
-
-    // Add layer addition/removal interfaces
-
-    topoChanger_.setSize(2);
-    label nMods = 0;
-
-    topoChanger_.set
-    (
-        0,
-        new layerAdditionRemoval
-        (
-            "right",
-            nMods,
-            topoChanger_,
-            "rightExtrusionFaces",
-            readScalar
-            (
-                motionDict_.subDict("right").lookup("minThickness")
-            ),
-            readScalar
-            (
-                motionDict_.subDict("right").lookup("maxThickness")
-            )
-        )
-    );
-    nMods++;
-
-    topoChanger_.set
-    (
-        1,
-        new layerAdditionRemoval
-        (
-            "left",
-            nMods,
-            topoChanger_,
-            "leftExtrusionFaces",
-            readScalar
-            (
-                motionDict_.subDict("left").lookup("minThickness")
-            ),
-            readScalar
-            (
-                motionDict_.subDict("left").lookup("maxThickness")
-            )
-        )
-    );
-    nMods++;
-
-    Info << "Adding " << nMods << " mesh modifiers" << endl;
-
-    // Write mesh and modifiers
-    topoChanger_.write();
-    write();
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-Foam::movingConeTopoFvMesh::movingConeTopoFvMesh(const IOobject& io)
-:
-    topoChangerFvMesh(io),
-    motionDict_
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "dynamicMeshDict",
-                time().constant(),
-                *this,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            )
-        ).subDict(typeName + "Coeffs")
-    ),
-    motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
-    motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
-    curMotionVel_
-    (
-        motionVelAmplitude_*
-        Foam::sin(time().value()*M_PI/motionVelPeriod_)
-    ),
-    leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
-    curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
-    curRight_(readScalar(motionDict_.lookup("rightObstacleEdge"))),
-    motionMask_
-    (
-        vertexMarkup
-        (
-            points(),
-            curLeft_,
-            curRight_
-        )
-    )
-{
-    addZonesAndModifiers();
-
-    curLeft_ = average
-    (
-        faceZones()
-        [
-            faceZones().findZoneID("leftExtrusionFaces")
-        ]().localPoints()
-    ).x();
-
-    curRight_ = average
-    (
-        faceZones()
-        [
-            faceZones().findZoneID("rightExtrusionFaces")
-        ]().localPoints()
-    ).x();
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::movingConeTopoFvMesh::~movingConeTopoFvMesh()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::movingConeTopoFvMesh::update()
-{
-    autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
-
-    // Calculate the new point positions depending on whether the
-    // topological change has happened or not
-    pointField newPoints;
-
-    vector curMotionVel_ =
-        motionVelAmplitude_*
-        Foam::sin(time().value()*M_PI/motionVelPeriod_);
-
-    bool meshChanged = topoChangeMap->morphing();
-
-    if (meshChanged)
-    {
-        Info << "Topology change. Calculating motion points" << endl;
-
-        if (topoChangeMap().hasMotionPoints())
-        {
-            motionMask_ =
-                vertexMarkup
-                (
-                    topoChangeMap().preMotionPoints(),
-                    curLeft_,
-                    curRight_
-                );
-        }
-        else
-        {
-            motionMask_ =
-                vertexMarkup
-                (
-                    points(),
-                    curLeft_,
-                    curRight_
-                );
-        }
-
-        // Create new points by moving old points but using the
-        // pre-motion points at the motion selector for the moving
-        // region
-        newPoints =
-            points()
-          + (
-                pos(0.5 - mag(motionMask_)) // cells above the body
-//               + pos(motionMask_ - 0.5)*      // cells in front of the body
-//                 (
-//                     points().component(vector::X)/curRight
-//                 )
-//               + pos(-motionMask_ - 0.5)*      // cells behind the body
-//                 (
-//                     (
-//                         points().component(vector::X)
-//                       - leftEdge
-//                     )/
-//                     (curLeft_ - leftEdge_)
-//                 )
-            )*curMotionVel_*time().deltaT().value();
-
-        // Correct mesh motion for correct volume continuity
-        movePoints(topoChangeMap().preMotionPoints());
-        resetMotion();
-        setV0();
-    }
-    else
-    {
-        Info << "No topology change" << endl;
-        // Set the mesh motion
-        newPoints =
-            points()
-          + (
-                pos(0.5 - mag(motionMask_)) // cells above the body
-//               + pos(motionMask_ - 0.5)*  // cells in front of the body
-//                 (
-//                     points().component(vector::X)/curRight
-//                 )
-//               + pos(-motionMask_ - 0.5)*  // cells behind the body
-//                 (
-//                     (
-//                         points().component(vector::X)
-//                       - leftEdge
-//                     )/
-//                     (curLeft_ - leftEdge_)
-//                 )
-            )*curMotionVel_*time().deltaT().value();
-    }
-
-    curLeft_ += curMotionVel_.x()*time().deltaT().value();
-    curRight_ += curMotionVel_.x()*time().deltaT().value();
-
-    // The mesh now contains the cells with zero volume
-
-    Info << "Executing mesh motion" << endl;
-    movePoints(newPoints);
-
-    //  The mesh now has got non-zero volume cells
-
-    return meshChanged;
-}
-
-
-// ************************************************************************* //

src/engine/engineTopoChangerMesh/engineValveSliding/setEngineValveSlidingConstraint.H

@@ -9,7 +9,6 @@
 
     forAll(valves_, valveI)
     {
-    
         vector valveVel =
             valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
 
@@ -17,10 +16,9 @@
         {
             valveVel = vector::zero;
         }
-    
+
-    
         {
-    
+
 //            label movingPtsIndex = pZones.findZoneID("movingPointsV"+Foam::name(valveI + 1));
 //            const labelList& movingPointsV = pZones[movingPtsIndex];
 /*
@@ -28,17 +26,17 @@
             {
                 constrainedPoints.append(movingPointsV[mpI]);
                 constrainedVelocity.append(valveVel);
-           
+
                 constraintSize++;
             }
 */
-            
+
             label movingCellsIndex = cZones.findZoneID("movingCellsZoneV"+Foam::name(valveI + 1));
             const labelList& movingCellsV = cZones[movingCellsIndex];
-            
+
             labelList movingPointsV;
             boolList count(newPoints.size(), false);
-            
+
             forAll(movingCellsV, cellI)
             {
                 const labelList& curCellPoints = cp[movingCellsV[cellI]];
@@ -76,17 +74,15 @@
             {
                 constrainedPoints.append(movingPointsV[mpI]);
                 constrainedVelocity.append(valveVel);
-           
+
                 constraintSize++;
             }
-    
         }
 
         {
-    
 //            label staticPtsIndex = pZones.findZoneID("staticPointsV"+Foam::name(valveI + 1));
 //            const labelList& staticPointsV = pZones[staticPtsIndex];
-/*    
+/*
             forAll(staticPointsV, spI)
             {
                 constrainedPoints.append(staticPointsV[spI]);
@@ -101,7 +97,7 @@
             labelList staticPointsV;
 
             boolList count(newPoints.size(), false);
-            
+
             forAll(staticCellsV, cellI)
             {
                 const labelList& curCellPoints = cp[staticCellsV[cellI]];
@@ -134,7 +130,7 @@
                     nCounted++;
                 }
             }
-            
+
             forAll(staticPointsV, spI)
             {
                 constrainedPoints.append(staticPointsV[spI]);
@@ -145,16 +141,16 @@
         }
 
     }
-    
+
     if(piston().patchID().active())
     {
         vector pistonVel =
             piston().cs().axis()*engineTime_.pistonSpeed().value();
 
-/*        
+/*
         label pistonPtsIndex = pZones.findZoneID("movingPistonPoints");
         const labelList& movingPointsP = pZones[pistonPtsIndex];
-    
+
         forAll(movingPointsP, mpI)
         {
             constrainedPoints.append(movingPointsP[mpI]);
@@ -165,11 +161,11 @@
 
         label pistonCellsIndex = cZones.findZoneID("movingPistonCells");
         const labelList& movingCellsP = cZones[pistonCellsIndex];
-        
+
         labelList movingPointsP;
-        
+
         boolList count(newPoints.size(), false);
-        
+
         forAll(movingCellsP, cellI)
         {
             const labelList& curCellPoints = cp[movingCellsP[cellI]];
@@ -202,8 +198,7 @@
                 nCounted++;
             }
         }
-        
+
-        
         forAll(movingPointsP, mpI)
         {
             constrainedPoints.append(movingPointsP[mpI]);
@@ -213,6 +208,4 @@
 
 
     }
-    
-
 }

src/engine/engineTopoChangerMesh/simpleTwoStroke/addSimpleTwoStrokeModifiers.C

@@ -301,8 +301,6 @@ void Foam::simpleTwoStroke::addZonesAndModifiers()
 
         nFaceZones++;
 
-        Info << "cut p" << endl;
-
         pz[nPointZones] = new pointZone
         (
             "cutPointZone",

src/engine/engineTopoChangerMesh/twoStrokeEngine/addPistonLayer.H

@@ -1,6 +1,11 @@
 if (piston().patchID().active())
 {
-    // add faces for piston layering
+    // Add faces for piston layering
+    // Note: because of operation of layer addition/removal
+    // (reduce function in layer addition/removal thickness)
+    // the layering modifier needs to be present on all processors
+    // even if the patch size is zero
+    // HJ, 7/Mar/2011
 
     faceSet pistonFaceSet(*this, piston().pistonFaceSetName());
 
@@ -23,7 +28,7 @@ if (piston().patchID().active())
     Info << "nSet = " << nSet << endl;
     Info << "nFlip = " << nFlip << endl;
 
-    fz.append
+    fz[nFaceZones] =
     (
         new faceZone
         (
@@ -37,27 +42,10 @@ if (piston().patchID().active())
 
     nFaceZones++;
 
-//     {
-//         pointSet movingPistonPoints(*this, piston().pistonPointSetName());
-
-//         pz.append
-//         (
-//             new pointZone
-//             (
-//                 "pistonPoints",
-//                 movingPistonPoints.toc(),
-//                 nPointZones,
-//                 pointZones()
-//             )
-//         );
-
-//         nPointZones++;
-//     }
-
     cellSet movingPistonCells(*this, piston().pistonCellSetName());
 
     Info<< "Adding piston cell set" << endl;
-    cz.append
+    cz[nCellZones] =
     (
         new cellZone
         (

src/engine/engineTopoChangerMesh/twoStrokeEngine/addTwoStrokeEngineModifiers.C

@@ -62,25 +62,24 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
         return;
     }
 
-    Info << "checkAndCalculate()" << endl;
     checkAndCalculate();
 
     Info<< "Time = " << engTime().theta() << endl
         << "Adding zones to the engine mesh" << endl;
 
-
+    // Zones to add
-    //fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
+    // fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
-    //pz = 2: piston points, cutPointZone
+    // pz = 2: piston points, cutPointZone
-    //cz = 1: moving mask
+    // cz = 1: moving mask
 
     label nPorts = scavInCylPatches_.size();
 
 
-    DynamicList<pointZone*> pz(2 + nPorts);
+    List<pointZone*> pz(nPorts + 2);
-    DynamicList<faceZone*> fz(3*nPorts + 1);
+    List<faceZone*> fz(3*nPorts + 1);
 
     // Added piston cells and head cells
-    DynamicList<cellZone*> cz(3);
+    List<cellZone*> cz(3);
 
     label nPointZones = 0;
     label nFaceZones = 0;
@@ -88,24 +87,11 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
     Info << "Adding piston layer faces" << endl;
 
-#   include "addPistonLayerHrv.H"
+#   include "addPistonLayer.H"
 
-//  adding head points (does not move)
+    // Add head points that do not move
 
     {
-//         pointSet headPointSet(*this, headPointsSetName_);
-
-//         pz[nPointZones] =
-//             new pointZone
-//             (
-//                 "headPoints",
-//                 headPointSet.toc(),
-//                 nPointZones,
-//                 pointZones()
-//             );
-
-//         nPointZones++;
-
         cellSet headCellSet(*this, headCellsSetName_);
 
         cz[nCellZones] =
@@ -120,9 +106,9 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
         nCellZones++;
     }
 
-//  Sliding interface for scavenging ports
+    //  Sliding interface for scavenging ports
 
-    if(nPorts > 0)
+    if (nPorts > 0)
     {
         forAll(scavInCylPatches_, patchi)
         {
@@ -133,76 +119,85 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
                     boundaryMesh().findPatchID(scavInCylPatches_[patchi])
                 ];
 
-            labelList isf(innerScav.size());
-
-            forAll (isf, i)
-            {
-                isf[i] = innerScav.start() + i;
-            }
-
-            fz[nFaceZones] = new faceZone
-            (
-                scavInCylPatches_[patchi] + "Zone" + Foam::name(patchi + 1),
-                isf,
-                boolList(innerScav.size(), false),
-                nFaceZones,
-                faceZones()
-            );
-
-            nFaceZones++;
-
             // Outer slider
-
             const polyPatch& outerScav =
                 boundaryMesh()
                 [
                     boundaryMesh().findPatchID(scavInPortPatches_[patchi])
                 ];
 
-            labelList osf(outerScav.size());
+            Pout<< "A: " << innerScav.size() << " " << outerScav.size()
-
+                << endl;
-            forAll (osf, i)
+            // Add zone if both patches has got faces
+            if (!innerScav.empty() && !outerScav.empty())
             {
-                osf[i] = outerScav.start() + i;
+                // Inner
+                labelList isf(innerScav.size());
+
+                forAll (isf, i)
+                {
+                    isf[i] = innerScav.start() + i;
+                }
+
+                fz[nFaceZones] = new faceZone
+                (
+                    scavInCylPatches_[patchi] + "Zone"
+                  + Foam::name(patchi + 1),
+                    isf,
+                    boolList(innerScav.size(), false),
+                    nFaceZones,
+                    faceZones()
+                );
+
+                nFaceZones++;
+
+                // Outer
+                labelList osf(outerScav.size());
+
+                forAll (osf, i)
+                {
+                    osf[i] = outerScav.start() + i;
+                }
+
+                fz[nFaceZones] = new faceZone
+                (
+                    scavInPortPatches_[patchi] + "Zone"
+                  + Foam::name(patchi + 1),
+                    osf,
+                    boolList(outerScav.size(), false),
+                    nFaceZones,
+                    faceZones()
+                );
+
+                nFaceZones++;
+
+                // Cut faces
+                fz[nFaceZones] = new faceZone
+                (
+                    "cutFaceZone" + Foam::name(patchi + 1),
+                    labelList(0),
+                    boolList(0, false),
+                    nFaceZones,
+                    faceZones()
+                );
+
+                nFaceZones++;
+
+                // Cut points
+                pz[nPointZones] = new pointZone
+                (
+                    "cutPointZone" + Foam::name(patchi + 1),
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                );
+
+                nPointZones++;
             }
-
-            fz[nFaceZones] = new faceZone
-            (
-                scavInPortPatches_[patchi] + "Zone" + Foam::name(patchi + 1),
-                osf,
-                boolList(outerScav.size(), false),
-                nFaceZones,
-                faceZones()
-            );
-
-            nFaceZones++;
-
-            fz[nFaceZones] = new faceZone
-            (
-                "cutFaceZone" + Foam::name(patchi + 1),
-                labelList(0),
-                boolList(0, false),
-                nFaceZones,
-                faceZones()
-            );
-
-            nFaceZones++;
-
-            Info << "cut p" << endl;
-
-            pz[nPointZones] = new pointZone
-            (
-                "cutPointZone" + Foam::name(patchi + 1),
-                labelList(0),
-                nPointZones,
-                pointZones()
-            );
-
-            nPointZones++;
         }
     }
 
-    Info << "moving cells" << endl;
+    Info << "Adding moving cells zone" << endl;
 
     {
         cellSet movingCells(*this, movingCellSetName_);
@@ -218,9 +213,7 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
         nCellZones++;
     }
 
-    Info << "moving cells added" << endl;
+    Pout<< "Adding " << nPointZones << " point, "
-
-    Info<< "Adding " << nPointZones << " point, "
         << nFaceZones << " face zones and " << nCellZones
         << " cell zones" << endl;
 
@@ -229,40 +222,64 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
     cz.setSize(nCellZones);
     addZones(pz, fz, cz);
 
-    List<polyMeshModifier*> tm(1 + nPorts);
+    List<polyMeshModifier*> tm(nPorts + 1);
     label nMods = 0;
 
     // Add piston layer addition
 
 #   include "addPistonLayerAdditionRemovalMeshModifier.H"
 
-    if(nPorts > 0)
+    if (nPorts > 0)
     {
-        forAll(scavInPortPatches_, i)
+        forAll (scavInPortPatches_, i)
         {
-            topoChanger_.setSize(topoChanger_.size() + 1);
+            // Check if patches are present on local processor
+            const label sipID =
+                boundaryMesh().findPatchID(scavInPortPatches_[i]);
 
-            topoChanger_.set
+            const label sicID =
-            (
+                boundaryMesh().findPatchID(scavInCylPatches_[i]);
-                nMods,
+
-                new slidingInterface
+            if (sipID > -1 && sicID > -1)
+            {
+                if
                 (
-                    "portCylinderInterface" + Foam::name(i + 1),
+                    boundaryMesh()[sipID].size() > 0
-                    nMods,
+                 && boundaryMesh()[sicID].size() > 0
-                    topoChanger_,
-                    scavInPortPatches_[i] + "Zone" + Foam::name(i + 1),
-                    scavInCylPatches_[i] + "Zone" + Foam::name(i + 1),
-                    "cutPointZone" + Foam::name(i + 1),
-                    "cutFaceZone" + Foam::name(i + 1),
-                    scavInPortPatches_[i],
-                    scavInCylPatches_[i],
-                    slidingInterface::INTEGRAL,
-                    true,                          // Attach-detach action
-                    intersection::VISIBLE         // Projection algorithm
                 )
-            );
+                {
+                    Pout<< "Adding slider for pair " << scavInPortPatches_[i]
+                        << " and " << scavInCylPatches_[i]
+                        << " with sizes "
+                        << boundaryMesh()[sipID].size() << " "
+                        << boundaryMesh()[sicID].size() << endl;
 
-            nMods++;
+                    // Patches present.  Add modifier
+                    topoChanger_.setSize(topoChanger_.size() + 1);
+
+                    topoChanger_.set
+                    (
+                        nMods,
+                        new slidingInterface
+                        (
+                            "portCylinderInterface" + Foam::name(i + 1),
+                            nMods,
+                            topoChanger_,
+                            scavInPortPatches_[i] + "Zone" + Foam::name(i + 1),
+                            scavInCylPatches_[i] + "Zone" + Foam::name(i + 1),
+                            "cutPointZone" + Foam::name(i + 1),
+                            "cutFaceZone" + Foam::name(i + 1),
+                            scavInPortPatches_[i],
+                            scavInCylPatches_[i],
+                            slidingInterface::INTEGRAL,
+                            true,                   // Attach-detach action
+                            intersection::VISIBLE   // Projection algorithm
+                        )
+                    );
+
+                    nMods++;
+                }
+            }
         }
     }
 

src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngine.C

@@ -152,13 +152,13 @@ void Foam::twoStrokeEngine::setBoundaryVelocity(volVectorField& U)
 
     //  On the piston movingWallVelocity is used.
     // There is no need to update the piston velocity
-//    U.boundaryField()[piston().patchID().index()] = pistonVel;
 
-    forAll(scavInPortPatches_, patchi)
+    forAll (scavInPortPatches_, patchi)
     {
-        U.boundaryField()
+        const label curPatchID =
-            [boundaryMesh().findPatchID(scavInPortPatches_[patchi])] ==
+            boundaryMesh().findPatchID(scavInPortPatches_[patchi]);
-            pistonVel;
+
+        U.boundaryField()[curPatchID] == pistonVel;
     }
 }
 

src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineCalculate.C

@@ -61,21 +61,26 @@ void Foam::twoStrokeEngine::calcMovingMasks() const
     const cellList& c = cells();
     const faceList& f = allFaces();
 
-    const labelList& cellAddr =
+    const label movingCellsID = cellZones().findZoneID("movingCells");
-        cellZones()[cellZones().findZoneID("movingCells")];
 
-    forAll (cellAddr, cellI)
+    // If moving cell zone is found, mark the vertices
+    if (movingCellsID > -1)
     {
-        const cell& curCell = c[cellAddr[cellI]];
+        const labelList& cellAddr = cellZones()[movingCellsID];
 
-        forAll (curCell, faceI)
+        forAll (cellAddr, cellI)
         {
-            // Mark all the points as moving
+            const cell& curCell = c[cellAddr[cellI]];
-            const face& curFace = f[curCell[faceI]];
 
-            forAll (curFace, pointI)
+            forAll (curCell, faceI)
             {
-                movingPointsMask[curFace[pointI]] = 1;
+                // Mark all the points as moving
+                const face& curFace = f[curCell[faceI]];
+
+                forAll (curFace, pointI)
+                {
+                    movingPointsMask[curFace[pointI]] = 1;
+                }
             }
         }
     }

src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineInitialize.C

@@ -44,9 +44,8 @@ void Foam::twoStrokeEngine::checkAndCalculate()
     label cylinderHeadIndex = -1;
     bool foundCylinderHead = false;
 
-    forAll(boundary(), i)
+    forAll (boundary(), i)
     {
-        Info << boundary()[i].name() << endl;
         if (boundary()[i].name() == "piston")
         {
             pistonIndex = i;
@@ -111,34 +110,31 @@ void Foam::twoStrokeEngine::checkAndCalculate()
     }
 }
 
+
 void Foam::twoStrokeEngine::setVirtualPistonPosition()
 {
-
     label pistonFaceIndex = faceZones().findZoneID("pistonLayerFaces");
 
-    bool foundPistonFace = (pistonFaceIndex != -1);
+    if(pistonFaceIndex == -1)
-
-    Info << "piston face index = " << pistonFaceIndex << endl;
-
-    if(!foundPistonFace)
     {
         FatalErrorIn("Foam::twoStrokeEngine::setVirtualPistonPosition()")
-            << " : cannot find the pistonLayerFaces"
+            << "Cannot find the pistonLayerFaces"
-            << exit(FatalError);
+            << abort(FatalError);
     }
 
-    const labelList& pistonFaces = faceZones()[pistonFaceIndex];
+    const labelList& pistonPoints =
-    forAll(pistonFaces, i)
+        faceZones()[pistonFaceIndex]().meshPoints();
-    {
-        const face& f = faces()[pistonFaces[i]];
 
-        // should loop over facepoints...
+    const pointField& p = points();
-        forAll(f, j)
+
-        {
+    forAll (pistonPoints, i)
-            virtualPistonPosition() =
+    {
-                Foam::max(virtualPistonPosition(), points()[f[j]].z());
+        virtualPistonPosition() =
-        }
+            Foam::max(virtualPistonPosition_, p[pistonPoints[i]].z());
     }
 
     reduce(virtualPistonPosition(), maxOp<scalar>());
 }
+
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineMove.C

@@ -41,11 +41,11 @@ void Foam::twoStrokeEngine::makeLayersLive()
     // Enable layering
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(layerAdditionRemoval))
+        if (isA<layerAdditionRemoval>(morphs[modI]))
         {
             morphs[modI].enable();
         }
-        else if (typeid(morphs[modI]) == typeid(slidingInterface))
+        else if (isA<slidingInterface>(morphs[modI]))
         {
             morphs[modI].disable();
         }
@@ -59,6 +59,7 @@ void Foam::twoStrokeEngine::makeLayersLive()
     }
 }
 
+
 void Foam::twoStrokeEngine::makeSlidersLive()
 {
     const polyTopoChanger& morphs = topoChanger_;
@@ -66,11 +67,11 @@ void Foam::twoStrokeEngine::makeSlidersLive()
     // Enable sliding interface
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(layerAdditionRemoval))
+        if (isA<layerAdditionRemoval>(morphs[modI]))
         {
             morphs[modI].disable();
         }
-        else if (typeid(morphs[modI]) == typeid(slidingInterface))
+        else if (isA<slidingInterface>(morphs[modI]))
         {
             morphs[modI].enable();
         }
@@ -93,10 +94,9 @@ bool Foam::twoStrokeEngine::attached() const
 
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(slidingInterface))
+        if (isA<slidingInterface>(morphs[modI]))
         {
-            result =
+            result = result
-                result
              || refCast<const slidingInterface>(morphs[modI]).attached();
         }
     }
@@ -104,7 +104,7 @@ bool Foam::twoStrokeEngine::attached() const
     // Check thal all sliders are in sync (debug only)
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(slidingInterface))
+        if (isA<slidingInterface>(morphs[modI]))
         {
             if
             (
@@ -120,6 +120,9 @@ bool Foam::twoStrokeEngine::attached() const
         }
     }
 
+    // Sync across processors
+    reduce(result, orOp<bool>());
+
     return result;
 }
 
@@ -133,13 +136,25 @@ bool Foam::twoStrokeEngine::update()
         makeSlidersLive();
 
         // Changing topology by hand
-        autoPtr<mapPolyMesh> topoChangeMap5 = topoChanger_.changeMesh();
+        autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
 
-        if (topoChangeMap5->hasMotionPoints() && topoChangeMap5->morphing())
+        bool localMorphing1 = topoChangeMap1->morphing();
+
+        // Note: Since we are detaching, global morphing is always true
+        // HJ, 7/Mar/2011
+
+        if (localMorphing1)
         {
             Info << "Topology change; executing pre-motion after "
                 << "sliding detach" << endl;
-            movePoints(topoChangeMap5->preMotionPoints());
+            movePoints(topoChangeMap1->preMotionPoints());
+        }
+        else
+        {
+            pointField newPoints = allPoints();
+
+            // Dummy motion
+            movePoints(newPoints);
         }
 
         Info << "sliding interfaces successfully decoupled!!!" << endl;
@@ -154,25 +169,14 @@ bool Foam::twoStrokeEngine::update()
     // Piston Layering
 
     makeLayersLive();
-    // Changing topology by hand
 
-
+    // Find piston mesh modifier if present on processor
-// /* Tommaso, 23/5/2008
-
-    // Find piston mesh modifier
     const label pistonLayerID =
         topoChanger_.findModifierID("pistonLayer");
 
-    if (pistonLayerID < 0)
-    {
-        FatalErrorIn("void engineFvMesh::moveAndMorph()")
-            << "Piston modifier not found."
-            << abort(FatalError);
-    }
-
     scalar minLayerThickness = piston().minLayer();
     scalar deltaZ = engTime().pistonDisplacement().value();
-    virtualPistonPosition() += deltaZ;
+    virtualPistonPosition_ += deltaZ;
 
     Info << "virtualPistonPosition = " << virtualPistonPosition()
     << ", deckHeight = " << deckHeight()
@@ -181,32 +185,45 @@ bool Foam::twoStrokeEngine::update()
     if (realDeformation())
     {
         // Dectivate piston layer
-        Info << "Mesh deformation mode" << endl;
+        if (pistonLayerID > -1)
-        topoChanger_[pistonLayerID].disable();
+        {
+            Info << "Mesh deformation mode" << endl;
+            topoChanger_[pistonLayerID].disable();
+        }
     }
     else
     {
         // Activate piston layer
-        Info << "Piston layering mode" << endl;
+        if (pistonLayerID > -1)
-        topoChanger_[pistonLayerID].enable();
+        {
+            Info << "Piston layering mode" << endl;
+            topoChanger_[pistonLayerID].enable();
+        }
     }
 
 
     // Changing topology by hand
-    autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
+    autoPtr<mapPolyMesh> topoChangeMap2 = topoChanger_.changeMesh();
+
+    bool localMorphing2 = topoChangeMap2->morphing();
+    bool globalMorphing2 = localMorphing2;
 
     // Work array for new points position.
     pointField newPoints = allPoints();
-    const pointField& refPoints = allPoints();
 
-    if (topoChangeMap->morphing())
+    if (globalMorphing2)
     {
-        if (topoChangeMap->hasMotionPoints())
+        Info<< "Topology change; executing pre-motion after "
+            << "dynamic layering" << endl;
+
+        if (localMorphing2)
         {
-            Info<< "Topology change; executing pre-motion after "
+            movePoints(topoChangeMap2->preMotionPoints());
-                << "dynamic layering" << endl;
+            newPoints = topoChangeMap2->preMotionPoints();
-            movePoints(topoChangeMap->preMotionPoints());
+        }
-            newPoints = topoChangeMap->preMotionPoints();
+        else
+        {
+            movePoints(newPoints);
         }
 
         setV0();
@@ -225,99 +242,92 @@ bool Foam::twoStrokeEngine::update()
     labelList pistonPoints;
     labelList headPoints;
     {
-        label pistonCellIndex = cellZones().findZoneID("pistonCells");
+        // Get cell-point addressing
-
-        if (pistonCellIndex < 0)
-        {
-            FatalErrorIn("bool twoStrokeEngine::update()")
-                << "Cannot find cell zone pistonCells"
-                << abort(FatalError);
-        }
-
-
-        const labelList& pistonCells = cellZones()[pistonCellIndex];
-
-        label headCellIndex = cellZones().findZoneID("headCells");
-
-        if (headCellIndex < 0)
-        {
-            FatalErrorIn("bool twoStrokeEngine::update()")
-                << "Cannot find cell zone headCells"
-                << abort(FatalError);
-        }
-
-        const labelList& headCells = cellZones()[headCellIndex];
-
         const labelListList& cp = cellPoints();
 
         boolList count(newPoints.size(), false);
 
-        forAll (pistonCells, cellI)
+        // Piston points
-        {
+        label pistonCellID = cellZones().findZoneID("pistonCells");
-            const labelList& curCellPoints = cp[pistonCells[cellI]];
 
-            forAll (curCellPoints, i)
+        if (pistonCellID > -1)
+        {
+            const labelList& pistonCells = cellZones()[pistonCellID];
+
+            forAll (pistonCells, cellI)
             {
-                count[curCellPoints[i]] = true;
+                const labelList& curCellPoints = cp[pistonCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
             }
-        }
 
-        // Count the points
+            // Count the points
-        label nCounted = 0;
+            label nCounted = 0;
-        forAll (count, pointI)
+            forAll (count, pointI)
-        {
-            if (count[pointI] == true)
             {
-                nCounted++;
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
             }
-        }
 
-        pistonPoints.setSize(nCounted);
+            pistonPoints.setSize(nCounted);
 
-        // Collect the points
+            // Collect the points
-        nCounted = 0;
+            nCounted = 0;
-        forAll (count, pointI)
+            forAll (count, pointI)
-        {
-            if (count[pointI] == true)
             {
-                pistonPoints[nCounted] = pointI;
+                if (count[pointI] == true)
-                nCounted++;
+                {
+                    pistonPoints[nCounted] = pointI;
+                    nCounted++;
+                }
             }
         }
 
         // Repeat for head points
         count = false;
 
-        forAll (headCells, cellI)
+        const label headCellID = cellZones().findZoneID("headCells");
-        {
-            const labelList& curCellPoints = cp[pistonCells[cellI]];
 
-            forAll (curCellPoints, i)
+        if (headCellID > -1)
+        {
+            const labelList& headCells = cellZones()[headCellID];
+
+            forAll (headCells, cellI)
             {
-                count[curCellPoints[i]] = true;
+                const labelList& curCellPoints = cp[headCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
             }
-        }
 
-        // Count the points
+            // Count the points
-        nCounted = 0;
+            label nCounted = 0;
-        forAll (count, pointI)
+            forAll (count, pointI)
-        {
-            if (count[pointI] == true)
             {
-                nCounted++;
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
             }
-        }
 
-        headPoints.setSize(nCounted);
+            headPoints.setSize(nCounted);
 
-        // Collect the points
+            // Collect the points
-        nCounted = 0;
+            nCounted = 0;
-        forAll (count, pointI)
+            forAll (count, pointI)
-        {
-            if (count[pointI] == true)
             {
-                headPoints[nCounted] = pointI;
+                if (count[pointI] == true)
-                nCounted++;
+                {
+                    headPoints[nCounted] = pointI;
+                    nCounted++;
+                }
             }
         }
     }
@@ -325,12 +335,6 @@ bool Foam::twoStrokeEngine::update()
 
     label nScaled = nPoints();
 
-//     label pistonPtsIndex = pointZones().findZoneID("pistonPoints");
-//     const labelList& pistonPoints = pointZones()[pistonPtsIndex];
-
-//     label headPtsIndex = pointZones().findZoneID("headPoints");
-//     const labelList& headPoints = pointZones()[headPtsIndex];
-
     const scalarField& movingPointsM = movingPointsMask();
 
     forAll(pistonPoints, i)
@@ -377,6 +381,7 @@ bool Foam::twoStrokeEngine::update()
     {
         // Always move piston
         scalar pistonTopZ = -GREAT;
+
         forAll(pistonPoints, i)
         {
             point& p = newPoints[pistonPoints[i]];
@@ -406,15 +411,11 @@ bool Foam::twoStrokeEngine::update()
         }
     }
 
-
-
     movePoints(newPoints);
     deleteDemandDrivenData(movingPointsMaskPtr_);
 
     pistonPosition() += deltaZ;
 
-//*/ //Tommaso, 23/5/2008
-
     {
         // Grab old points to correct the motion
         pointField oldPointsNew = oldAllPoints();
@@ -424,73 +425,98 @@ bool Foam::twoStrokeEngine::update()
         makeSlidersLive();
 
         // Changing topology by hand
-        autoPtr<mapPolyMesh> topoChangeMap4 = topoChanger_.changeMesh();
+        autoPtr<mapPolyMesh> topoChangeMap3 = topoChanger_.changeMesh();
 
-        if (topoChangeMap4->morphing())
+        bool localMorphing3 = topoChangeMap3->morphing();
+        bool globalMorphing3 = localMorphing3;
+
+        reduce(globalMorphing3, orOp<bool>());
+
+        if (globalMorphing3)
         {
-            if (topoChangeMap4->hasMotionPoints())
+            Info<< "Topology change; executing pre-motion after "
+                << "sliding attach" << endl;
+
+            // Grab points
+            newPoints = allPoints();
+
+            if (localMorphing3)
             {
-                Info<< "Topology change; executing pre-motion after "
+                // If there is layering, pick up correct points
-                    << "sliding attach" << endl;
+                if (topoChangeMap3->hasMotionPoints())
+                {
+                    newPoints = topoChangeMap3->preMotionPoints();
+                }
 
-//                 Info<< "topoChangeMap4->preMotionPoints().size() = "
+                // Prepare old points for the move
-//                     << topoChangeMap4->preMotionPoints().size() << nl
-//                     << "allPoints.size() = " << allPoints().size() << nl
-//                     << "points.size() = " << points().size() << endl;
-
-                movePoints(topoChangeMap4->preMotionPoints());
-                newPoints = points();
-            }
-
-            {
                 pointField mappedOldPointsNew(allPoints().size());
 
                 mappedOldPointsNew.map
                 (
-                    oldPointsNew, topoChangeMap4->pointMap()
+                    oldPointsNew, topoChangeMap3->pointMap()
                 );
 
                 forAll(scavInPortPatches_, patchi)
                 {
-                    const labelList& cutPointsAddressing =
+                    // Find cut point zone ID
-                        pointZones()
+                    const label cutPointZoneID = pointZones().findZoneID
-                        [
+                    (
-                            pointZones().findZoneID
+                        "cutPointZone" + Foam::name(patchi + 1)
-                            (
+                    );
-                                "cutPointZone" + Foam::name(patchi + 1)
-                            )
-                        ];
 
-                    forAll(cutPointsAddressing, i)
+                    if (cutPointZoneID > -1)
                     {
-                        mappedOldPointsNew[cutPointsAddressing[i]] =
+                        const labelList& cutPointsAddressing =
-                            newPoints[cutPointsAddressing[i]];
+                            pointZones()[cutPointZoneID];
-                    }
 
-                    forAll(cutPointsAddressing, i)
+                        forAll(cutPointsAddressing, i)
-                    {
-                        if
-                        (
-                            newPoints[cutPointsAddressing[i]].z()
-                          > virtualPistonPosition()
-                        )
                         {
-                            mappedOldPointsNew[cutPointsAddressing[i]].z() =
+                            mappedOldPointsNew[cutPointsAddressing[i]] =
-                                newPoints[cutPointsAddressing[i]].z();
+                                newPoints[cutPointsAddressing[i]];
                         }
-                        else
+
+                        forAll(cutPointsAddressing, i)
                         {
-                            mappedOldPointsNew[cutPointsAddressing[i]].z() =
+                            if
-                                newPoints[cutPointsAddressing[i]].z() - deltaZ;
+                            (
+                                newPoints[cutPointsAddressing[i]].z()
+                              > virtualPistonPosition()
+                            )
+                            {
+                                mappedOldPointsNew
+                                    [cutPointsAddressing[i]].z() =
+                                    newPoints[cutPointsAddressing[i]].z();
+                            }
+                            else
+                            {
+                                mappedOldPointsNew
+                                    [cutPointsAddressing[i]].z() =
+                                    newPoints[cutPointsAddressing[i]].z()
+                                  - deltaZ;
+                            }
                         }
                     }
                 }
-                pointField newPoints = allPoints();
 
+                // Move mesh into correct old configuration
                 movePoints(mappedOldPointsNew);
 
                 resetMotion();
                 setV0();
+
+                // Set new point motion
+                movePoints(newPoints);
+            }
+            else
+            {
+                // No local topological change.  Execute double motion for
+                // sync with topological changes
+                movePoints(oldPointsNew);
+
+                resetMotion();
+                setV0();
+
+                // Set new point motion
                 movePoints(newPoints);
             }
         }

src/engine/include/addPistonLayerAdditionRemovalMeshModifier.H

@@ -1,11 +1,11 @@
     // Add piston layer addition/removal
     if (piston().patchID().active())
     {
+        Info<< "Adding a layer addition/removal mesh modifier to the piston"
+            << endl;
+
+        topoChanger_.setSize(nMods + 1);
 
-        Info << "Adding a layer addition/removal mesh modifier to the piston" << endl;
-    
-        topoChanger_.setSize(nMods+1);
-         
         topoChanger_.set
         (
             nMods,
@@ -17,9 +17,9 @@
                 "pistonLayerFaces",
                 piston().minLayer(),
                 piston().maxLayer()
-            )    
+            )
         );
         nMods++;
-        Info << "pistonLayer" << endl;
+
-        Info << nMods << endl;            
+        Info << "pistonLayer modifier number " << nMods << endl;
     }

src/engine/include/createEngineDynamicMesh.H

@@ -7,7 +7,8 @@
             (
                 engineTopoChangerMesh::defaultRegion,
                 runTime.timeName(),
-                runTime
+                runTime,
+                IOobject::MUST_READ
             )
         )
     );

src/engine/simpleEngineTopoFvMesh/simpleEngineTopoFvMesh.C

@@ -43,7 +43,7 @@ namespace Foam
 
     addToRunTimeSelectionTable
     (
-        topoChangerFvMesh,
+        engineTopoChangerMesh,
         simpleEngineTopoFvMesh,
         IOobject
     );
@@ -211,7 +211,7 @@ bool Foam::simpleEngineTopoFvMesh::attached() const
         {
             if
             (
-                result 
+                result
              != refCast<const slidingInterface>(topoChanger_[modI]).attached()
             )
             {
@@ -458,8 +458,7 @@ Foam::simpleEngineTopoFvMesh::simpleEngineTopoFvMesh
     const IOobject& io
 )
 :
-    topoChangerFvMesh(io),
+    engineTopoChangerMesh(io),
-    engineTime_(refCast<const engineTime>(time())),
     valves_(*this, engineTime_.engineDict().lookup("valves")),
     piston_(*this, engineTime_.engineDict().subDict("piston")),
     msPtr_(motionSolver::New(*this)),

src/engine/simpleEngineTopoFvMesh/simpleEngineTopoFvMesh.H

@@ -27,8 +27,7 @@ License
 #ifndef simpleEngineTopoFvMesh_H
 #define simpleEngineTopoFvMesh_H
 
-#include "topoChangerFvMesh.H"
+#include "engineTopoChangerMesh.H"
-#include "engineTime.H"
 #include "valveBank.H"
 #include "simpleEnginePiston.H"
 #include "motionSolver.H"
@@ -44,13 +43,10 @@ namespace Foam
 
 class simpleEngineTopoFvMesh
 :
-    public topoChangerFvMesh
+    public engineTopoChangerMesh
 {
     // Private data
 
-        //- Engine database
-        const engineTime& engineTime_;
-
         //- Engine valves
         valveBank valves_;
 
@@ -76,6 +72,9 @@ class simpleEngineTopoFvMesh
         void operator=(const simpleEngineTopoFvMesh&);
 
 
+        //- Add valve and piston zones and modifiers
+        void addZonesAndModifiers();
+
         //- Make layering modifiers live
         void makeLayersLive();
 
@@ -140,11 +139,12 @@ public:
              && engineTime_.thetaRevolution() < deformSwitch_;
         }
 
-        //- Add valve and piston zones and modifiers
-        void addZonesAndModifiers();
-
         //- Update the mesh for both mesh motion and topology change
         virtual bool update();
+
+        //- Set boundary velocities
+        virtual void setBoundaryVelocity(volVectorField& U)
+        {}
 };
 
 

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C

@@ -360,7 +360,7 @@ edgeInterpolationScheme<Type>::interpolate
 
             for (label i=0; i<size; i++)
             {
-                const tensorField& curT = 
+                const tensorField& curT =
                     mesh.edgeTransformTensors()[start + i];
 
                 const tensor& Te = curT[0];
@@ -372,13 +372,13 @@ edgeInterpolationScheme<Type>::interpolate
                     (
                         Te.T(),
                         pLambda[i]*transform(TP, pOwnVf[i])
-                      + (1.0 - pLambda[i])*transform(TN, pNgbVf[i])
+                      + (1 - pLambda[i])*transform(TN, pNgbVf[i])
                     );
             }
 
 //             tsf().boundaryField()[pi] =
 //                 pLambda*vf.boundaryField()[pi].patchInternalField()
-//              + (1.0 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
+//              + (1 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
         }
         else
         {

src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C

@@ -118,7 +118,7 @@ void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    // a simpler way of doing this would be nice
+    // A simpler way of doing this would be nice
     scalar avgU = -flowRate_/gSum(patch().magSf());
 
     vectorField n = patch().nf();
@@ -128,7 +128,7 @@ void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs()
 
     if (phi.dimensions() == dimVelocity*dimArea)
     {
-        // volumetric flow-rate
+        // Volumetric flow-rate
         operator==(n*avgU);
     }
     else if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
@@ -136,7 +136,7 @@ void Foam::flowRateInletVelocityFvPatchVectorField::updateCoeffs()
         const fvPatchField<scalar>& rhop =
             patch().lookupPatchField<volScalarField, scalar>(rhoName_);
 
-        // mass flow-rate
+        // Mass flow-rate
         operator==(n*avgU/rhop);
     }
     else

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

@@ -69,10 +69,10 @@ tmp<volScalarField> CoEulerDdtScheme<Type>::CorDeltaT() const
 
     forAll(owner, faceI)
     {
-        corDeltaT[owner[faceI]] = 
+        corDeltaT[owner[faceI]] =
             max(corDeltaT[owner[faceI]], cofrDeltaT[faceI]);
 
-        corDeltaT[neighbour[faceI]] = 
+        corDeltaT[neighbour[faceI]] =
             max(corDeltaT[neighbour[faceI]], cofrDeltaT[faceI]);
     }
 

src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateScheme/multivariateScheme.H

@@ -105,7 +105,7 @@ public:
         public:
 
             // Constructors
-            
+
                 fieldScheme
                 (
                     const GeometricField<Type, fvPatchField, volMesh>& field,

src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/tetPolyPatches/constraint/processor/calcProcessorTetPolyPatchCellDecompPointAddr.C

@@ -59,8 +59,7 @@ void Foam::processorTetPolyPatchCellDecomp::calcMeshPoints() const
             masterFaces[faceI] = pp[faceI].reverseFace();
         }
 
-        mp =
+        mp = primitiveFacePatch
-            primitiveFacePatch
             (
                 masterFaces,
                 pp.points()

tutorials/incompressible/simpleFoam/pitzDaily/system/fvSolution

@@ -63,6 +63,8 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+
+    convergence          1e-4;
 }
 
 relaxationFactors

tutorials/mesh/moveDynamicMesh/circCylinder3d/system/tetFemSolution

@@ -16,7 +16,14 @@ FoamFile
 
 solvers
 {
-    motionU ICCG 1e-9 0.001;
+    motionU
+    {
+        solver           CG;
+        preconditioner   Cholesky;
+
+        tolerance        1e-09;
+        relTol           0.001;
+    }
 }
 
 // ************************************************************************* //