Updates on engine motion

applications/utilities/mesh/manipulation/moveDyMEngineMesh/Make/options

@@ -11,4 +11,5 @@ EXE_LIBS = \
     -lmeshTools \
     -ldynamicMesh \
     -lengine \
-    -lmeshTools
+    -lmeshTools \
+    -llduSolvers

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

@@ -56,10 +56,11 @@ int main(int argc, char *argv[])
             "U",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedVector("zero", dimVelocity, vector::zero)
     );
 
     volScalarField motionContErr

src/engine/Make/files

@@ -9,10 +9,14 @@ engineVerticalValve/engineVerticalValve.C
 thoboisSlidingValve/thoboisSlidingValve.C
 dieselEngineValve/dieselEngineValve.C
 thoboisValve/thoboisValve.C
+accordionValve/accordionValve.C
 
 simpleEnginePiston/simpleEnginePiston.C
 enginePiston/enginePiston.C
 
+engineTopoChangerMesh/regionSide/regionSide.C
+engineTopoChangerMesh/attachDetachFunctions/attachDetachFunctions.C
+
 engineMesh/engineMesh/engineMesh.C
 engineMesh/engineMesh/newEngineMesh.C
 engineMesh/staticEngineMesh/staticEngineMesh.C
@@ -28,10 +32,14 @@ engineTopoChangerMesh/engineTopoChangerMesh/engineTopoChangerMesh.C
 engineTopoChangerMesh/engineTopoChangerMesh/newEngineTopoChangerMesh.C
 
 engineTopoChangerMesh/accordionEngineMesh/accordionEngineMesh.C
-engineTopoChangerMesh/accordionEngineMesh/accordionEngineMeshCalculate.C
 engineTopoChangerMesh/accordionEngineMesh/accordionEngineMeshInitialize.C
 engineTopoChangerMesh/accordionEngineMesh/accordionEngineMeshMove.C
-engineTopoChangerMesh/accordionEngineMesh/addAccordionEngineMeshModifiers.C
+engineTopoChangerMesh/accordionEngineMesh/addAccordionEngineMeshZones.C
+
+engineTopoChangerMesh/deformingEngineMesh/deformingEngineMesh.C
+engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshInitialize.C
+engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshMove.C
+engineTopoChangerMesh/deformingEngineMesh/addDeformingEngineMeshZones.C
 
 engineTopoChangerMesh/engineValveSliding/addEngineValveSlidingMeshModifiers.C
 engineTopoChangerMesh/engineValveSliding/engineValveSliding.C

src/engine/accordionValve/accordionValve.C

@@ -0,0 +1,265 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "accordionValve.H"
+#include "engineTime.H"
+#include "polyMesh.H"
+#include "interpolateXY.H"
+#include "IFstream.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::accordionValve::adjustCrankAngle(const scalar theta) const
+{
+    if (theta < liftProfileStart_)
+    {
+        scalar adjustedTheta = theta;
+
+        while (adjustedTheta < liftProfileStart_)
+        {
+            adjustedTheta += liftProfileEnd_ - liftProfileStart_;
+        }
+
+        return adjustedTheta;
+    }
+    else if (theta > liftProfileEnd_)
+    {
+        scalar adjustedTheta = theta;
+
+        while (adjustedTheta > liftProfileEnd_)
+        {
+            adjustedTheta -= liftProfileEnd_ - liftProfileStart_;
+        }
+
+        return adjustedTheta;
+    }
+    else
+    {
+        return theta;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::accordionValve::accordionValve
+(
+    const word& name,
+    const polyMesh& mesh,
+    const autoPtr<coordinateSystem>& valveCS,
+    const word& bottomPatchName,
+    const word& poppetPatchName,
+    const word& sidePatchName,
+    const word& stemPatchName,
+    const word& detachInCylinderPatchName,
+    const word& detachInPortPatchName,
+    const word& detachFacesName,
+    const graph& liftProfile,
+    const scalar minLift,
+    const scalar diameter,
+    const word& staticCellsName,
+    const word& movingCellsName
+)
+:
+    name_(name),
+    mesh_(mesh),
+    engineDB_(refCast<const engineTime>(mesh.time())),
+    csPtr_(valveCS),
+    bottomPatch_(bottomPatchName, mesh.boundaryMesh()),
+    poppetPatch_(poppetPatchName, mesh.boundaryMesh()),
+    sidePatch_(sidePatchName, mesh.boundaryMesh()),
+    stemPatch_(stemPatchName, mesh.boundaryMesh()),
+    detachInCylinderPatch_(detachInCylinderPatchName, mesh.boundaryMesh()),
+    detachInPortPatch_(detachInPortPatchName, mesh.boundaryMesh()),
+    detachFacesName_(detachFacesName),
+    liftProfile_(liftProfile),
+    liftProfileStart_(min(liftProfile_.x())),
+    liftProfileEnd_(max(liftProfile_.x())),
+    minLift_(minLift),
+    diameter_(diameter),
+    staticCellsName_(staticCellsName),
+    movingCellsName_(movingCellsName)
+{}
+
+
+// Construct from dictionary
+Foam::accordionValve::accordionValve
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    name_(name),
+    mesh_(mesh),
+    engineDB_(refCast<const engineTime>(mesh_.time())),
+    csPtr_
+    (
+        coordinateSystem::New
+        (
+            "coordinateSystem",
+            dict.subDict("coordinateSystem")
+        )
+    ),
+    bottomPatch_(dict.lookup("bottomPatch"), mesh.boundaryMesh()),
+    poppetPatch_(dict.lookup("poppetPatch"), mesh.boundaryMesh()),
+    sidePatch_(dict.lookup("sidePatch"), mesh.boundaryMesh()),
+    stemPatch_(dict.lookup("stemPatch"), mesh.boundaryMesh()),
+    detachInCylinderPatch_
+    (
+        dict.lookup("detachInCylinderPatch"),
+        mesh.boundaryMesh()
+    ),
+    detachInPortPatch_
+    (
+        dict.lookup("detachInPortPatch"),
+        mesh.boundaryMesh()
+    ),
+    detachFacesName_(dict.lookup("detachFaces")),
+    liftProfile_
+    (
+        "theta",
+        "lift",
+        name_,
+        IFstream
+        (
+            mesh.time().path()/mesh.time().caseConstant()/
+            word(dict.lookup("liftProfileFile"))
+        )()
+    ),
+    liftProfileStart_(min(liftProfile_.x())),
+    liftProfileEnd_(max(liftProfile_.x())),
+    minLift_(readScalar(dict.lookup("minLift"))),
+    diameter_(readScalar(dict.lookup("diameter"))),
+    staticCellsName_(dict.lookup("staticCells")),
+    movingCellsName_(dict.lookup("movingCells"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::accordionValve::lift(const scalar theta) const
+{
+
+    label curCycle = theta/720.0;
+    scalar curTheta = theta - curCycle * 720.0;
+
+    return interpolateXY
+    (
+        curTheta,
+        liftProfile_.x(),
+        liftProfile_.y()
+    );
+}
+
+
+bool Foam::accordionValve::isOpen() const
+{
+    return lift(engineDB_.theta()) >= minLift_;
+}
+
+
+Foam::scalar Foam::accordionValve::curLift() const
+{
+    return max
+    (
+        lift(engineDB_.theta()),
+        minLift_
+    );
+}
+
+
+Foam::scalar Foam::accordionValve::curVelocity() const
+{
+    return
+       -(
+             curLift()
+           - max
+             (
+                 lift(engineDB_.theta() - engineDB_.deltaTheta()),
+                 minLift_
+             )
+        )/(engineDB_.deltaT().value() + VSMALL);
+}
+
+Foam::labelList Foam::accordionValve::movingPatchIDs() const
+{
+    labelList mpIDs(3);
+    label nMpIDs = 0;
+
+    if (bottomPatch_.active())
+    {
+        mpIDs[nMpIDs] = bottomPatch_.index();
+        nMpIDs++;
+    }
+
+    if (poppetPatch_.active())
+    {
+        mpIDs[nMpIDs] = poppetPatch_.index();
+        nMpIDs++;
+    }
+
+    if (sidePatch_.active())
+    {
+        mpIDs[nMpIDs] = sidePatch_.index();
+        nMpIDs++;
+    }
+
+    mpIDs.setSize(nMpIDs);
+
+    return mpIDs;
+}
+
+
+void Foam::accordionValve::writeDict(Ostream& os) const
+{
+    os  << nl << name() << nl << token::BEGIN_BLOCK;
+
+    cs().writeDict(os);
+
+    os  << "bottomPatch " << bottomPatch_.name() << token::END_STATEMENT << nl
+        << "poppetPatch " << poppetPatch_.name() << token::END_STATEMENT << nl
+        << "sidePatch " << sidePatch_.name() << token::END_STATEMENT << nl
+        << "stemPatch " << stemPatch_.name() << token::END_STATEMENT << nl
+        << "detachInCylinderPatch " << detachInCylinderPatch_.name()
+        << token::END_STATEMENT << nl
+        << "detachInPortPatch " << detachInPortPatch_.name()
+        << token::END_STATEMENT << nl
+        << "detachFaces " << detachFacesName_ << token::END_STATEMENT << nl
+        << "liftProfile " << nl << token::BEGIN_BLOCK
+        << liftProfile_ << token::END_BLOCK << token::END_STATEMENT << nl
+        << "minLift " << minLift_ << token::END_STATEMENT << nl
+        << "diameter " << diameter_ << token::END_STATEMENT << nl
+        << "staticCells " << staticCellsName_ << token::END_STATEMENT << nl
+        << "movingCells " << movingCellsName_ << token::END_STATEMENT << nl
+        << token::END_BLOCK << endl;
+}
+
+// ************************************************************************* //

src/engine/accordionValve/accordionValve.H

@@ -0,0 +1,284 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    accordionValve
+
+Description
+
+SourceFiles
+    accordionValve.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef accordionValve_H
+#define accordionValve_H
+
+#include "word.H"
+#include "coordinateSystem.H"
+#include "polyPatchID.H"
+#include "graph.H"
+
+#include "faceSet.H"
+#include "pointSet.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+class polyMesh;
+class engineTime;
+
+/*---------------------------------------------------------------------------*\
+                          Class accordionValve Declaration
+\*---------------------------------------------------------------------------*/
+
+class accordionValve
+{
+    // Private data
+
+        //- Name of valve
+        word name_;
+
+        //- Reference to engine mesh
+        const polyMesh& mesh_;
+
+        //- Reference to engine database
+        const engineTime& engineDB_;
+
+        //- Coordinate system
+        autoPtr<coordinateSystem> csPtr_;
+
+
+        // Patch and zone names
+
+            //- Valve bottom patch
+            polyPatchID bottomPatch_;
+
+            //- Valve poppet patch
+            polyPatchID poppetPatch_;
+
+            //- Valve poppet patch
+            polyPatchID sidePatch_;
+
+            //- Valve stem patch
+            polyPatchID stemPatch_;
+
+            //- Valve detach in cylinder patch
+            polyPatchID detachInCylinderPatch_;
+
+            //- Valve detach in port patch
+            polyPatchID detachInPortPatch_;
+
+            //- Faces to detach
+            //labelList detachFaces_;
+            word detachFacesName_;
+
+        // Valve lift data
+
+            //- Valve lift profile
+            graph liftProfile_;
+
+            //- Lift curve start angle
+            scalar liftProfileStart_;
+
+            //- Lift curve end angle
+            scalar liftProfileEnd_;
+
+            //- Minimum valve lift.  On this lift the valve is considered closed
+            const scalar minLift_;
+
+            //- Valve diameter
+            const scalar diameter_;
+
+        // Mesh motion data
+
+            //- Name of the cellSet for the static cells
+            word staticCellsName_;
+
+            //- Name of the cellSet for the moving cells
+            word movingCellsName_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        accordionValve(const accordionValve&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const accordionValve&);
+
+
+        //- Adjust crank angle to drop within the limits of the lift profile
+        scalar adjustCrankAngle(const scalar theta) const;
+
+public:
+
+    // Static data members
+
+    // Constructors
+
+        //- Construct from components
+        accordionValve
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const autoPtr<coordinateSystem>& valveCS,
+            const word& bottomPatchName,
+            const word& poppetPatchName,
+            const word& sidePatchName,
+            const word& stemPatchName,
+            const word& detachInCylinderPatchName,
+            const word& detachInPortPatchName,
+            const word& detachFacesName,
+            const graph& liftProfile,
+            const scalar minLift,
+            const scalar diameter,
+            const word& staticCellsName,
+            const word& movingCellsName
+        );
+
+        //- Construct from dictionary
+        accordionValve
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+    // Destructor - default
+
+
+    // Member Functions
+
+        //- Return name
+        const word& name() const
+        {
+            return name_;
+        }
+
+        //- Return coordinate system
+        const coordinateSystem& cs() const
+        {
+            return csPtr_();
+        }
+
+        //- Return lift profile
+        const graph& liftProfile() const
+        {
+            return liftProfile_;
+        }
+
+        //- Return valve diameter
+        scalar diameter() const
+        {
+            return diameter_;
+        }
+
+        // Valve patches
+
+            //- Return ID of bottom patch
+            const polyPatchID& bottomPatchID() const
+            {
+                return bottomPatch_;
+            }
+
+            //- Return ID of poppet patch
+            const polyPatchID& poppetPatchID() const
+            {
+                return poppetPatch_;
+            }
+
+            //- Return ID of side patch
+            const polyPatchID& sidePatchID() const
+            {
+                return sidePatch_;
+            }
+
+            //- Return ID of stem patch
+            const polyPatchID& stemPatchID() const
+            {
+                return stemPatch_;
+            }
+
+            //- Return ID of detach in cylinder patch
+            const polyPatchID& detachInCylinderPatchID() const
+            {
+                return detachInCylinderPatch_;
+            }
+
+            //- Return ID of detach in port patch
+            const polyPatchID& detachInPortPatchID() const
+            {
+                return detachInPortPatch_;
+            }
+
+            //- Return face labels of detach curtain
+            const word& detachFacesName() const
+            {
+                return detachFacesName_;
+            }
+
+            const word& staticCellsName() const
+            {
+                return staticCellsName_;
+            }
+
+            const word& movingCellsName() const
+            {
+                return movingCellsName_;
+            }
+
+        // Valve position and velocity
+
+            //- Return valve lift given crank angle in degrees
+            scalar lift(const scalar theta) const;
+
+            //- Is the valve open?
+            bool isOpen() const;
+
+            //- Return current lift
+            scalar curLift() const;
+
+            //- Return valve velocity for current time-step
+            scalar curVelocity() const;
+
+            //- Return list of active patch labels for the valve head
+            //  (stem is excluded)
+            labelList movingPatchIDs() const;
+
+        //- Write dictionary
+        void writeDict(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/engine/accordionValve/accordionValveBank.H

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    accordionValveBank
+
+Description
+    A list of valves.
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef accordionValveBank_H
+#define accordionValveBank_H
+
+#include "PtrList.H"
+#include "accordionValve.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class accordionValveBank Declaration
+\*---------------------------------------------------------------------------*/
+
+class accordionValveBank
+:
+    public PtrList<accordionValve>
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        accordionValveBank(const accordionValveBank&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const accordionValveBank&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from Istream
+        accordionValveBank
+        (
+            const polyMesh& mesh,
+            Istream& is
+        )
+        {
+            PtrList<entry> valveEntries(is);
+
+            setSize(valveEntries.size());
+
+            forAll (valveEntries, valveI)
+            {
+                set
+                (
+                    valveI,
+                    new accordionValve
+                    (
+                        valveEntries[valveI].keyword(),
+                        mesh,
+                        valveEntries[valveI].dict()
+                    )
+                );
+            }
+        }
+
+    // Destructor - default
+
+    // Member Functions
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/accordionEngineMesh/accordionEngineMesh.C

@@ -24,9 +24,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+
 #include "accordionEngineMesh.H"
-#include "layerAdditionRemoval.H"
-#include "attachDetach.H"
 #include "componentMixedTetPolyPatchVectorField.H"
 #include "mapPolyMesh.H"
 #include "polyTopoChange.H"
@@ -45,7 +44,12 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(accordionEngineMesh, 0);
-    addToRunTimeSelectionTable(engineTopoChangerMesh, accordionEngineMesh, IOobject);
+    addToRunTimeSelectionTable
+    (
+        engineTopoChangerMesh,
+        accordionEngineMesh,
+        IOobject
+    );
 }
 
 
@@ -57,7 +61,11 @@ namespace Foam
 bool Foam::accordionEngineMesh::realDeformation() const
 {
 
-    if (virtualPistonPosition() + engTime().pistonDisplacement().value() > deckHeight_ - SMALL)
+    if
+    (
+        virtualPistonPosition() + engTime().pistonDisplacement().value()
+      > deckHeight_ - SMALL
+    )
     {
         return true;
     }
@@ -73,32 +81,26 @@ bool Foam::accordionEngineMesh::realDeformation() const
 Foam::accordionEngineMesh::accordionEngineMesh
 (
     const IOobject& io
+//    bool addZonesAndMods
 )
 :
     engineTopoChangerMesh(io),
     piston_(*this, engTime().engineDict().subDict("piston")),
     valves_(*this, engTime().engineDict().lookup("accordionEngineMesh")),
-    deformSwitch_(readScalar(engTime().engineDict().lookup("deformAngle"))),
-    delta_(readScalar(engTime().engineDict().lookup("delta"))),
-    offSet_(readScalar(engTime().engineDict().lookup("offSet"))),
     pistonPosition_(-GREAT),
     virtualPistonPosition_(-GREAT),
     deckHeight_(GREAT),
-    msPtr_(motionSolver::New(*this)),
     cylinderHeadName_(engTime().engineDict().lookup("cylinderHeadName")),
     linerName_(engTime().engineDict().lookup("linerName")),
-    pistonAuxPoints_(engTime().engineDict().lookup("pistonAuxPoints")),
-    moveDetach_(engTime().engineDict().lookup("moveDetach"))
+    headCellSetName_(engTime().engineDict().lookup("headCellSetName"))
 
 {
     // Add zones and modifiers if not already there.
-    addZonesAndModifiers();
+    addMeshZones();
+    msPtr_ = motionSolver::New(*this);
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
@@ -118,32 +120,8 @@ void Foam::accordionEngineMesh::setBoundaryVelocity(volVectorField& U)
             U.boundaryField()[valves()[valveI].stemPatchID().index()] ==
                 valveVel;
         }
-
-        // If valve is present in geometry, set the motion
-        if (valves()[valveI].detachInPortPatchID().active())
-        {
-            // Bottom of the valve moves with given velocity
-            U.boundaryField()[valves()[valveI].detachInPortPatchID().index()] ==
-                vector::zero;
-            U.oldTime().boundaryField()[valves()[valveI].detachInPortPatchID().index()] ==
-                vector::zero;
     }
-
-        // If valve is present in geometry, set the motion
-        if (valves()[valveI].detachInCylinderPatchID().active())
-        {
-            // Bottom of the valve moves with given velocity
-            U.boundaryField()[valves()[valveI].detachInCylinderPatchID().index()] ==
-               vector::zero;
-            U.oldTime().boundaryField()[valves()[valveI].detachInCylinderPatchID().index()] ==
-               vector::zero;
-        }
-     
-    }
-
 }
 
 
-
-
 // ************************************************************************* //

src/engine/engineTopoChangerMesh/accordionEngineMesh/accordionEngineMesh.H

@@ -26,13 +26,11 @@ CLASS
     accordionEngineMesh
 
 DESCRIPTION
-    Engine mesh class for canted valves engine. Layers are added/removed on the
-    piston, while mesh is deformed close to the valves. Attach/detach boundary
-    is used. There are no sliding interfaces.
 
 */
 // ------------------------------------------------------------------------- //
 
+
 #ifndef accordionEngineMesh_H
 #define accordionEngineMesh_H
 
@@ -40,7 +38,7 @@ DESCRIPTION
 #include "enginePiston.H"
 #include "motionSolver.H"
 #include "polyPatchID.H"
-#include "thoboisValveBank.H"
+#include "accordionValveBank.H"
 #include "twoDPointCorrector.H"
 
 #include "faceSet.H"
@@ -69,16 +67,7 @@ class accordionEngineMesh
         enginePiston piston_;
 
         //- Lateral valves
-        thoboisValveBank valves_;
-                         
-        //- Layering-to-deformation switch in crank-angle degrees
-        scalar deformSwitch_;
-        
-        //- Tolerance used when the piston faceZone is created
-        scalar delta_;
-
-        //- Distance from the piston patch
-        scalar offSet_;
+        accordionValveBank valves_;
 
         //- Piston Position
         scalar pistonPosition_;
@@ -99,9 +88,11 @@ class accordionEngineMesh
         word linerName_;
 
         //- PistonFaceSet
-        word pistonAuxPoints_;
+//        word pistonAuxPoints_;
+
+        //- PistonFaceSet
+        word headCellSetName_;
 
-	Switch moveDetach_;
 
     // Private Member Functions
 
@@ -111,30 +102,16 @@ class accordionEngineMesh
         //- Disallow default bitwise assignment
         void operator=(const accordionEngineMesh&);
 
-        //- Make layering modifiers live
-        void makeLayersLive();
-
-        //- Make attach-detach live
-        void makeDetachLive();
-       
-        //- Prepare valve attach/detach
-        void valveDetach();
-
-        //- Prepare valve attach/detach
-        void valveAttach();
-
-        //- Prepare valve attach/detach
-        void prepareValveDetach();
-
         //- Check if all patches exist, then calculate virtualPistonPosition,
         //- pistonPosition and deckHeight for the first time
         void checkAndCalculate();
 
-        //- Calculate the virtualPistonPosition, 
+        //- Calculate the virtualPistonPosition
         void setVirtualPositions();
 
         bool realDeformation() const;
 
+
 public:
 
     //- Runtime type information
@@ -147,16 +124,16 @@ public:
         explicit accordionEngineMesh(const IOobject& io);
 
 
-    // Destructor - default
+    // Destructor
+
+        virtual ~accordionEngineMesh()
+        {}
 
 
     // Member Functions
 
         inline const enginePiston& piston() const;
-        inline const thoboisValveBank& valves() const;
-        inline const scalar& deformSwitch() const;
-        inline const scalar& delta() const;
-        inline const scalar& offSet() const;
+        inline const accordionValveBank& valves() const;
         inline const scalar& pistonPosition() const;
         inline scalar& pistonPosition();
         inline const scalar& virtualPistonPosition() const;
@@ -168,9 +145,7 @@ public:
         //- Return true for mesh deformation mode
         bool deformation() const
         {
-            return
-                engTime().thetaRevolution() > -deformSwitch_
-             && engTime().thetaRevolution() < deformSwitch_;
+            return true;
         }
 
         //- Return number of valves
@@ -180,7 +155,7 @@ public:
         }
 
         //- Add valve and piston zones and modifiers
-        void addZonesAndModifiers();
+        void addMeshZones();
 
         //- Move and morph
         virtual bool update();

src/engine/engineTopoChangerMesh/accordionEngineMesh/accordionEngineMeshI.H

@@ -37,26 +37,11 @@ inline const enginePiston& accordionEngineMesh::piston() const
 }
 
 
-inline const thoboisValveBank& accordionEngineMesh::valves() const
+inline const accordionValveBank& accordionEngineMesh::valves() const
 {
     return valves_;
 }
                  
-inline const scalar& accordionEngineMesh::deformSwitch() const
-{
-    return deformSwitch_;
-}
-
-inline const scalar& accordionEngineMesh::delta() const
-{
-    return delta_;
-}
-
-inline const scalar& accordionEngineMesh::offSet() const
-{
-    return offSet_;
-}
-
 
 inline const scalar& accordionEngineMesh::pistonPosition() const
 {

src/engine/engineTopoChangerMesh/accordionEngineMesh/accordionEngineMeshInitialize.C

@@ -25,8 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "accordionEngineMesh.H"
-#include "slidingInterface.H"
-#include "layerAdditionRemoval.H"
 #include "surfaceFields.H"
 #include "regionSplit.H"
 

src/engine/engineTopoChangerMesh/accordionEngineMesh/accordionEngineMeshMove.C

@@ -25,9 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "accordionEngineMesh.H"
-#include "slidingInterface.H"
-#include "layerAdditionRemoval.H"
-#include "attachDetach.H"
 #include "surfaceFields.H"
 #include "regionSplit.H"
 #include "componentMixedTetPolyPatchVectorField.H"
@@ -56,210 +53,6 @@ License
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-void Foam::accordionEngineMesh::makeLayersLive()
-{
-    // Enable layering
-    forAll (topoChanger_, modI)
-    {
-        if (isA<layerAdditionRemoval>(topoChanger_[modI]))
-        {
-            topoChanger_[modI].enable();
-        }
-        else if (isA<slidingInterface>(topoChanger_[modI]))
-        {
-            topoChanger_[modI].disable();
-        }
-        else if (isA<attachDetach>(topoChanger_[modI]))
-        {
-            topoChanger_[modI].enable();
-        }
-        else
-        {
-            FatalErrorIn("void Foam::engineTopoFvMesh::makeLayersLive()")
-                << "Don't know what to do with mesh modifier "
-                << modI << " of type " << topoChanger_[modI].type()
-                << abort(FatalError);
-        }
-    }
-}
-
-void Foam::accordionEngineMesh::makeDetachLive()
-{
-    // Enable sliding interface
-    forAll (topoChanger_, modI)
-    {
-        if (isA<attachDetach>(topoChanger_[modI]))
-        {
-            topoChanger_[modI].enable();
-        }
-        else
-        {
-            FatalErrorIn("void Foam::accordionEngineMesh::makeDetachLive()")
-                << "Don't know what to do with mesh modifier "
-                << modI << " of type " << topoChanger_[modI].type()
-                << abort(FatalError);
-        }
-    }
-}
-
-
-void Foam::accordionEngineMesh::valveDetach()
-{
-    // Enable sliding interface
-    forAll (topoChanger_, modI)
-    {
-        if (isA<attachDetach>(topoChanger_[modI]))
-        {
-            const attachDetach& ad =
-                refCast<const attachDetach>(topoChanger_[modI]);
-
-            const word masterName = ad.masterPatchID().name();
-
-            // Find the valve with that name
-            label valveIndex = -1;
-
-            forAll (valves_, valveI)
-            {
-                if
-                (
-                    valves_[valveI].detachInCylinderPatchID().name()
-                 == masterName
-                )
-                {
-                    valveIndex = valveI;
-                    break;
-                }
-            }
-
-            if (valveIndex < 0)
-            {
-                FatalErrorIn
-                (
-                    "void Foam::engineTopoFvMesh::prepareValveDetach()"
-                )   << "Cannot match patch for attach/detach " << modI
-                    << abort(FatalError);
-            }
-
-            if (debug)
-            {
-                Info<< " valveI: " << valveIndex << " attached: "
-                    << ad.attached()
-                    << " valve open: " << valves_[valveIndex].isOpen()
-                    << endl;
-            }
-
-            ad.setDetach();
-        }
-    }
-}
-
-void Foam::accordionEngineMesh::valveAttach()
-{
-    // Enable sliding interface
-    forAll (topoChanger_, modI)
-    {
-        if (isA<attachDetach>(topoChanger_[modI]))
-        {
-            const attachDetach& ad =
-                refCast<const attachDetach>(topoChanger_[modI]);
-
-            const word masterName = ad.masterPatchID().name();
-
-            // Find the valve with that name
-            label valveIndex = -1;
-
-            forAll (valves_, valveI)
-            {
-                if
-                (
-                    valves_[valveI].detachInCylinderPatchID().name()
-                 == masterName
-                )
-                {
-                    valveIndex = valveI;
-                    break;
-                }
-            }
-
-            if (valveIndex < 0)
-            {
-                FatalErrorIn
-                (
-                    "void Foam::engineTopoFvMesh::prepareValveDetach()"
-                )   << "Cannot match patch for attach/detach " << modI
-                    << abort(FatalError);
-            }
-
-            if (debug)
-            {
-                Info<< " valveI: " << valveIndex << " attached: "
-                    << ad.attached()
-                    << " valve open: " << valves_[valveIndex].isOpen()
-                    << endl;
-            }
-
-            ad.setAttach();
-        }
-    }
-}
-
-
-void Foam::accordionEngineMesh::prepareValveDetach()
-{
-    // Enable sliding interface
-    forAll (topoChanger_, modI)
-    {
-        if (isA<attachDetach>(topoChanger_[modI]))
-        {
-            const attachDetach& ad =
-                refCast<const attachDetach>(topoChanger_[modI]);
-
-            const word masterName = ad.masterPatchID().name();
-
-            // Find the valve with that name
-            label valveIndex = -1;
-
-            forAll (valves_, valveI)
-            {
-                if
-                (
-                    valves_[valveI].detachInCylinderPatchID().name()
-                 == masterName
-                )
-                {
-                    valveIndex = valveI;
-                    break;
-                }
-            }
-
-            if (valveIndex < 0)
-            {
-                FatalErrorIn
-                (
-                    "void Foam::engineTopoFvMesh::prepareValveDetach()"
-                )   << "Cannot match patch for attach/detach " << modI
-                    << abort(FatalError);
-            }
-
-            if (debug)
-            {
-                Info<< " valveI: " << valveIndex << " attached: "
-                    << ad.attached()
-                    << " valve open: " << valves_[valveIndex].isOpen()
-                    << endl;
-            }
-
-            if (valves_[valveIndex].isOpen())
-            {
-                ad.setAttach();
-            }
-            else
-            {
-                ad.setDetach();
-            }
-        }
-    }
-}
 
 
 bool Foam::accordionEngineMesh::update()
@@ -267,10 +60,6 @@ bool Foam::accordionEngineMesh::update()
     tetDecompositionMotionSolver& mSolver =
         refCast<tetDecompositionMotionSolver>(msPtr_());
 
-    tetPointVectorField& motionU = mSolver.motionU();
-
-    Info << "motioU.size() = " << motionU.internalField().size() << endl;
-
     scalar deltaZ = engTime().pistonDisplacement().value();
 
     // deltaZ set to zero, FIXED PISTON POSITION
@@ -278,31 +67,12 @@ bool Foam::accordionEngineMesh::update()
 
     virtualPistonPosition() += deltaZ;
 
-    makeDetachLive();
-
-//    valveDetach();
-
-    autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
-
-    if (topoChangeMap->morphing())
-    {
-        mSolver.updateMesh(topoChangeMap());
-        Info << "mSolver.updateMesh(topoChangeMap())" << endl;
-
-        if (topoChangeMap->hasMotionPoints())
-        {
-            movePoints(topoChangeMap->preMotionPoints());
-            resetMotion();
-            setV0();
-        }
-    }
-
     pointField newPoints = points();
 
     {
 #       include "setValveMotionBoundaryConditionAccordionEngineMesh.H"
 #       include "setPistonMotionBoundaryConditionAccordionEngineMesh.H"
-
+        Info << "piston motion" << endl;
 
         DynamicList<label> constrainedPoints(mSolver.curPoints()().size()/100);
         DynamicList<vector> constrainedVelocity
@@ -312,7 +82,6 @@ bool Foam::accordionEngineMesh::update()
 
 #       include "setAccordionEngineMeshConstraints.H"
 
-
         labelList constrainedPointsList(constrainedPoints.shrink());
         vectorField constrainedVelocityField(constrainedVelocity.shrink());
 
@@ -329,88 +98,25 @@ bool Foam::accordionEngineMesh::update()
 
         //set to zero the motion U along the x and y directions
 
+
         newPoints = mSolver.curPoints();
         movePoints(newPoints);
         setVirtualPositions();
+
         mSolver.clearConstraints();
+
     }
 
 //    pointField oldPointsNew = oldPoints();
     pointField oldPointsNew = oldAllPoints();
     newPoints = points();
-
+    movePoints(newPoints);
     Info << "max mesh phi before = " << max((phi())) << endl;
     Info << "min mesh phi before = " << min((phi())) << endl;
 
-//    makeDetachLive();
-    prepareValveDetach();
-
-    autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
-
-    if (topoChangeMap1->morphing())
-    {
-        mSolver.updateMesh(topoChangeMap1());
-
-        if (topoChangeMap1->hasMotionPoints())
-        {
-//              movePoints(topoChangeMap1->preMotionPoints());
-//              resetMotion();
-//              setV0();
-
-            // correct the motion after attaching the sliding interface
-            pointField mappedOldPointsNew(allPoints().size());
-
-            mappedOldPointsNew.map(oldPointsNew, topoChangeMap1->pointMap());
-            pointField newPoints = allPoints();
-
-            movePoints(mappedOldPointsNew);
-            resetMotion();
-            setV0();
-            movePoints(newPoints);
-
-        }
-    }
 
     Info << "max mesh phi after = " << max((phi())) << endl;
     Info << "min mesh phi after = " << min((phi())) << endl;
 
-/*         if(correctPointsMotion_)
-        {
-            // correct the motion after attaching the sliding interface
-            pointField mappedOldPointsNew(allPoints().size());
-
-            mappedOldPointsNew.map(oldPointsNew, topoChangeMap3->pointMap());
-            pointField newPoints = allPoints();
-
-            movePoints(mappedOldPointsNew);
-            resetMotion();
-            setV0();
-            movePoints(newPoints);
-        }
-
-//    Info << motionU << endl;
-
-    if (meshChanged1)
-    {
-            Info << "meshChanged1" << endl;
-
-        mSolver.updateMesh(topoChangeMap1());
-
-       {
-           pointField mappedOldPointsNew(allPoints().size());
-
-           mappedOldPointsNew.map(oldPointsNew, topoChangeMap1->pointMap());
-           pointField newPoints = allPoints();
-
-           movePoints(mappedOldPointsNew);
-           resetMotion();
-           setV0();
-           movePoints(newPoints);
-       }
-    }
-*/
-    Info<< "Total cylinder volume at CA " << engTime().timeName() << " = "
-        << sum(V()) << endl;
-
     return false;
 }

src/engine/engineTopoChangerMesh/accordionEngineMesh/addAccordionEngineMeshZones.C

@@ -25,25 +25,33 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "accordionEngineMesh.H"
-#include "slidingInterface.H"
-#include "layerAdditionRemoval.H"
 #include "surfaceFields.H"
 #include "regionSplit.H"
-#include "attachDetach.H"
+#include "regionSide.H"
+#include "attachDetachFunctions.H"
+#include "directTopoChange.H"
+#include "zeroGradientTetPolyPatchFields.H"
+#include "tetDecompositionMotionSolver.H"
+
+#include "fixedValueTetPolyPatchFields.H"
+#include "mixedTetPolyPatchFields.H"
+#include "slipTetPolyPatchFields.H"
+#include "zeroGradientTetPolyPatchFields.H"
+#include "meshTools.H"
+#include "syncTools.H"
+#include "HashSet.H"
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::accordionEngineMesh::addZonesAndModifiers()
+void Foam::accordionEngineMesh::addMeshZones()
 {
     // Add the zones and mesh modifiers to operate piston motion
+    Switch addBoundary(true);
 
     if
     (
         pointZones().size() > 0
-/*
-     && faceZones().size() > 0
-     && cellZones().size() > 0
-     && topoChanger_.size() > 0
-*/
     )
     {
         Info<< "Time = " << engTime().theta() << endl;
@@ -51,15 +59,12 @@ void Foam::accordionEngineMesh::addZonesAndModifiers()
             << "Zones and modifiers already present.  Skipping."
             << endl;
 
-     //   setVirtualPositions();
         checkAndCalculate();
 
         Info << "Point zones found = " << pointZones().size() << endl;
         Info << "Face zones found = " << faceZones().size() << endl;
         Info << "Cell zones found = " << cellZones().size() << endl;
 
-        return;
-
     }
     else
     {
@@ -67,7 +72,6 @@ void Foam::accordionEngineMesh::addZonesAndModifiers()
         faceZones().setSize(0);
         cellZones().setSize(0);
         topoChanger_.setSize(0);
-    }
 
         Info << "checkAndCalculate()" << endl;
         checkAndCalculate();
@@ -75,54 +79,29 @@ void Foam::accordionEngineMesh::addZonesAndModifiers()
         Info<< "Time = " << engTime().theta() << endl
             << "Adding zones to the engine mesh" << endl;
 
-/*
-    Point zones
-    1) Piston points
-*/
-
         DynamicList<pointZone*> pz;
 
-/*
-    Face zones
-    1) Piston layer faces
-
-*/
         DynamicList<faceZone*> fz;
 
-/*
-*/
-
         DynamicList<cellZone*> cz;
 
         label nPointZones = 0;
         label nFaceZones = 0;
         label nCellZones = 0;
 
-/*
-
-    Adding the following faces zones:
-    1:  pistonLayerFaces
-    nv: movingFaces
-    nv: staticFaces
-    nv: detachFaces
-
-    Adding the following point zones:
-    1:  pistonPoints
-    nv:  movingValvePoints
-    nv:  staticValvePoints
-
-*/
 
 #       include "addPistonFacesPointZonesAccordionEngineMesh.H"
 
-#   include "addAttachDetachFacesAccordionEngineMesh.H"
 #       include "addValveFaceZonesAccordionEngineMesh.H"
 
 #       include "addOutputCellsAccordionEngineMesh.H"
+
+
         Info<< "Adding " << nPointZones << " point, "
             << nFaceZones << " face zones and "
             << nCellZones << " cell zones" << endl;
 
+        Info << boundary().size() << endl;
 
         pz.setSize(nPointZones);
         Info << "setSize pz" << endl;
@@ -132,17 +111,15 @@ void Foam::accordionEngineMesh::addZonesAndModifiers()
         Info << "setSize cz" << endl;
 
         addZones(pz, fz, cz);
+    }
 
-#   include "addMeshModifiersAccordionEngineMesh.H"
 
     // Calculating the virtual positions of piston and valves
 
  //   setVirtualPositions();
 
     // Write mesh and modifiers
-    topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
     write();
-
     Info << "virtualPistonPosition = " << virtualPistonPosition() << endl;
     Info << "piston position = " << pistonPosition() << endl;
 }

src/engine/engineTopoChangerMesh/accordionEngineMesh/addAttachDetachFacesAccordionEngineMesh.H

@@ -1,49 +0,0 @@
-
-    for (label valveI = 0; valveI < nValves(); valveI++)
-    {
-        
-        label nFlipDetach = 0;
-        
-        faceSet detachFaces(*this, valves_[valveI].detachFacesName());
-        
-        boolList flipDetach(detachFaces.size(), false);
-                   
-        const polyPatch& bottomPatch =
-        boundaryMesh()[valves()[valveI].bottomPatchID().index()];
-        
-        vector bottomPatchCenter = sum(mag(bottomPatch.faceAreas()) *
-        bottomPatch.faceCentres())/sum(mag(bottomPatch.faceAreas()));
-                
-        forAll(detachFaces.toc(), i)
-        {
-            vector radius = faceCentres()[detachFaces.toc()[i]] - bottomPatchCenter;
-            
-            radius/= mag(radius);
-            
-            if(((faceAreas()[detachFaces.toc()[i]]) & (radius) / mag(faceAreas()[detachFaces.toc()[i]])) < 0)
-            {
-//                flipDetach[i] = true;
-//                nFlipDetach++;
-            }
-               
-        }
-        
-        Info << "valve n. " << valveI << " found " << nFlipDetach << " flipped faces for detach" << endl;
-                
-        fz.append
-        (
-            new faceZone
-            (
-                "detachFaceZoneV" + Foam::name(valveI + 1),
-                detachFaces.toc(),
-                flipDetach,
-                nFaceZones,
-                faceZones()
-            )
-        );
-        
-        nFaceZones++;     
-    }   
-
-    Info << "Added attach/detach faces zones for valves" << endl;
-    

src/engine/engineTopoChangerMesh/accordionEngineMesh/addMeshModifiersAccordionEngineMesh.H

@@ -1,70 +0,0 @@
-//    List<polyMeshModifier*> tm(3*nValves() + 2);
-//    DynamicList<polyMeshModifier*> tm;
-
-    label nMods = 0;
-
-    for (label valveI = 0; valveI < nValves(); valveI++)
-    {
-        
-        
-        // Add attach-detach for valve
-        if
-        (
-            valves_[valveI].detachInCylinderPatchID().active()
-         && valves_[valveI].detachInPortPatchID().active()
-        )
-        {
-
-            Info << "Adding a attach/detach boundary mesh modifier to the valve " <<
-            valveI + 1 << endl;
-        
-            topoChanger_.setSize(nMods+1);
-            topoChanger_.set
-            (
-                nMods,
-                new attachDetach
-                (
-                    "valveAttachDetach" + Foam::name(valveI + 1),
-                    nMods,
-                    topoChanger_,
-                    "detachFaceZoneV" + Foam::name(valveI + 1),
-                    valves_[valveI].detachInCylinderPatchID().name(),
-                    valves_[valveI].detachInPortPatchID().name(),
-                    scalarField(0),
-                    true                // Manual triggering
-                )    
-            );
-            nMods++;
-            Info << "Attach detach" << endl;
-        }        
-                   
-    }
-
-    // Add piston layer addition
-    if (piston().patchID().active())
-    {
-
-/*
-        Info << "Adding a layer addition/removal mesh modifier to the piston" << endl;
-    
-        topoChanger_.setSize(nMods+1);
-        topoChanger_.set
-        (
-            nMods,
-            new layerAdditionRemoval
-            (
-                "pistonLayer",
-                nMods,
-                topoChanger_,
-                "pistonLayerFaces",
-                piston().minLayer(),
-                piston().maxLayer()
-            )    
-        );
-        nMods++;
-        Info << "pistonLayer" << endl;
-        Info << nMods << endl;            
-*/
-    }
-
-    Info << "Adding " << nMods << " topology modifiers" << endl;

src/engine/engineTopoChangerMesh/accordionEngineMesh/addOutputCellsAccordionEngineMesh.H

@@ -1,21 +1,55 @@
-{
-    word cylinderOutputName
+    {
+        word cylinderOutputName(engTime().engineDict().lookup("cylinderSetName"));
+
+        IOobject outputHeader
         (
-        engTime().engineDict().lookup("cylinderSetName")
+            cylinderOutputName,
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
         );
 
-    cellSet outputCellSet(*this, cylinderOutputName);
+        if(!outputHeader.headerOk())
+        {
+            WarningIn
+            (
+                "accordionEngineMesh::addZonesAndModifiers()"
+            )   << "faceSet called " << cylinderOutputName
+                << " does not exist. Continuing mesh motion without cylinder output cells" << endl;
 
         cz.append
         (
           new cellZone
           (
                 "cylinderOutputCells",
-            outputCellSet.toc(),
+                labelList(0),
                 nCellZones,
                 cellZones()
           )
       );
         
       nCellZones++;        
+  
+
+        }
+        else
+        {
+
+        cellSet outputCellSet(*this, cylinderOutputName);
+
+        labelList outputCellSetList = outputCellSet.toc();
+
+        cz.append
+        (
+          new cellZone
+          (
+                "cylinderOutputCells",
+                outputCellSetList,
+                nCellZones,
+                cellZones()
+          )
+      );
+
+      nCellZones++;        
+    }
 }

src/engine/engineTopoChangerMesh/accordionEngineMesh/addPistonFacesPointZonesAccordionEngineMesh.H

@@ -2,14 +2,121 @@
     // Add the piston points zone
     if (piston().patchID().active())
     {
-        pointSet movingPistonPoints(*this, pistonAuxPoints_);
+        IOobject pistonCellsHeader
+        (
+            piston().pistonCellSetName(),
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+        
+        if(!pistonCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "accordionEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << piston().pistonCellSetName()
+                << " does not exist. Continuing mesh motion without piston points" << endl;
 
             pz.append
             (
                 new pointZone
                 (
                     "movingPistonPoints",
-                movingPistonPoints.toc(),
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                )
+            );
+        
+            nPointZones++;        
+                
+        }
+        else
+        {
+            cellSet movingPistonCells(*this, piston().pistonCellSetName());
+
+            labelList pistonCells = movingPistonCells.toc();
+            
+            labelList movingPistonPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (pistonCells, cellI)
+            {   
+                const labelList& curCellPoints = cp[pistonCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }   
+
+            movingPistonPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    movingPistonPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "movingPistonPoints",
+                    movingPistonPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+        
+            nPointZones++;        
+                        
+        }
+                           
+    }
+
+    {
+
+        IOobject headCellsHeader
+        (
+            headCellSetName_,
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+
+        if(!headCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "accordionEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << headCellSetName_
+                << " does not exist. Continuing mesh motion without head points" << endl;
+            pz.append
+            (
+                new pointZone
+                (
+                    "headPoints",
+                    labelList(0),
                     nPointZones,
                     pointZones()
                 )
@@ -17,3 +124,63 @@
         
             nPointZones++;        
         }
+        else
+        {
+            cellSet fixedHeadCells(*this, headCellSetName_);
+
+            labelList headCells = fixedHeadCells.toc();
+            
+            labelList fixedHeadPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (fixedHeadCells, cellI)
+            {   
+                const labelList& curCellPoints = cp[headCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }   
+
+            fixedHeadPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    fixedHeadPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "headPoints",
+                    fixedHeadPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+        
+            nPointZones++;        
+        }
+    
+    }

src/engine/engineTopoChangerMesh/accordionEngineMesh/addValveFaceZonesAccordionEngineMesh.H

@@ -9,79 +9,28 @@
         //  - staticFaceZone
         //  - staticPointZone
         
-/*
-        faceSet movingFaces(*this, valves_[valveI].movingFacesName());
-        
-        fz.append
+        IOobject movingCellsHeader
         (
-            new faceZone
-            (
-                "movingFacesZoneV" + Foam::name(valveI + 1),
-                movingFaces.toc(),
-                boolList(movingFaces.size(), false),
-                nFaceZones,
-                faceZones()
-            )
+            valves_[valveI].movingCellsName(),
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
         );
         
-        nFaceZones++;
-        
-        faceSet staticFaces(*this, valves_[valveI].staticFacesName());
-        
-        fz.append
+        if(!movingCellsHeader.headerOk())
+        {
+            WarningIn
             (
-            new faceZone
-            (
-                "staticFacesZoneV" + Foam::name(valveI + 1),
-                staticFaces.toc(),
-                boolList(staticFaces.size(), false),
-                nFaceZones,
-                faceZones()
-            )
-        );
-        
-        nFaceZones++;
-*/
-
-        cellSet movingCells(*this, valves_[valveI].movingCellsName());
-        
-        cz.append
-        (
-            new cellZone
-            (
-                "movingCellsZoneV" + Foam::name(valveI + 1),
-                movingCells.toc(),
-                nCellZones,
-                cellZones()
-            )
-        );
-        
-        nCellZones++;
-        
-        cellSet staticCells(*this, valves_[valveI].staticCellsName());
-        
-        cz.append
-        (
-            new cellZone
-            (
-                "staticCellsZoneV" + Foam::name(valveI + 1),
-                staticCells.toc(),
-                nCellZones,
-                cellZones()
-            )
-        );
-        
-        nCellZones++;
-        
-
-        pointSet movingPoints(*this, valves_[valveI].movingPointsName());
+                "accordionEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << valves_[valveI].movingCellsName()
+                << " does not exist. Continuing mesh motion without rigid motion points for valve " << valves_[valveI].name() << endl;
 
             pz.append
             (
                 new pointZone
                 (
                     "movingPointsV" + Foam::name(valveI + 1),
-                movingPoints.toc(),
+                    labelList(0),
                     nPointZones,
                     pointZones()
                 )
@@ -89,29 +38,59 @@
         
             nPointZones++;        
 
-        pointSet staticPoints(*this, valves_[valveI].staticPointsName());
+
+        }
+        else
+        {
+            cellSet movingCellsSet(*this, valves_[valveI].movingCellsName());
+
+            labelList movingCells = movingCellsSet.toc();
+            
+            labelList movingPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (movingCells, cellI)
+            {   
+                const labelList& curCellPoints = cp[movingCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }   
+
+            movingPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    movingPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
 
             pz.append
             (
                 new pointZone
                 (
-                "staticPointsV" + Foam::name(valveI + 1),
-                staticPoints.toc(),
-                nPointZones,
-                pointZones()
-            )
-        );
-        
-        nPointZones++;     
-
-        pointSet movingInternalPoints(*this, valves_[valveI].movingInternalPointsName());
-
-        pz.append
-        (
-            new pointZone
-            (
-                "movingInternalPointsV" + Foam::name(valveI + 1),
-                movingInternalPoints.toc(),
+                    "movingPointsV" + Foam::name(valveI + 1),
+                    movingPoints,
                     nPointZones,
                     pointZones()
                 )
@@ -119,4 +98,94 @@
         
             nPointZones++;        
         }
+
+        IOobject staticCellsHeader
+        (
+            valves_[valveI].staticCellsName(),
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+        
+        if(!staticCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "accordionEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << valves_[valveI].staticCellsName()
+                << " does not exist. Continuing mesh motion without static points for valve " << valves_[valveI].name() << endl;
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "staticPointsV" + Foam::name(valveI + 1),
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                )
+            );
+        
+            nPointZones++;     
+
+        }
+        else
+        {
+            cellSet staticCellSet(*this, valves_[valveI].staticCellsName());
+
+            labelList staticCells = staticCellSet.toc();
+            
+            labelList staticPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (staticCells, cellI)
+            {   
+                const labelList& curCellPoints = cp[staticCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }   
+
+            staticPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    staticPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "staticPointsV" + Foam::name(valveI + 1),
+                    staticPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+        
+            nPointZones++;     
+        }
+    }   
 }

src/engine/engineTopoChangerMesh/accordionEngineMesh/setAccordionEngineMeshConstraints.H

@@ -1,4 +1,5 @@
 {
+    boolList isConstrained(points().size(), false);
 
     const pointZoneMesh& pZones = pointZones();
 
@@ -15,33 +16,16 @@
             const labelList& movingPointsV = pZones[movingPtsIndex];
 
             forAll(movingPointsV, mpI)
+            {
+                if(!isConstrained[movingPointsV[mpI]])
                 {
                     constrainedPoints.append(movingPointsV[mpI]);
                     constrainedVelocity.append(valveVel);
 
                     constraintSize++;
+                    isConstrained[movingPointsV[mpI]] = true;
                 }
             }
-
-        {
-
-            label movingInternalPtsIndex =
-                pZones.findZoneID
-                (
-                    "movingInternalPointsV"
-                  + Foam::name(valveI + 1)
-                );
-
-            const labelList& movingInternalPointsV =
-                pZones[movingInternalPtsIndex];
-
-            forAll(movingInternalPointsV, mipI)
-            {
-                constrainedPoints.append(movingInternalPointsV[mipI]);
-                constrainedVelocity.append(valveVel);
-
-                constraintSize++;
-            }
         }
 
         {
@@ -51,11 +35,14 @@
             const labelList& staticPointsV = pZones[staticPtsIndex];
 
             forAll(staticPointsV, spI)
+            {
+                if(!isConstrained[staticPointsV[spI]])
                 {
                     constrainedPoints.append(staticPointsV[spI]);
                     constrainedVelocity.append(vector::zero);
                     constraintSize++;
-
+                    isConstrained[staticPointsV[spI]] = true;
+                }
             }
 
         }
@@ -69,15 +56,36 @@
             const labelList& movingPointsP = pZones[pistonPtsIndex];
 
             forAll(movingPointsP, mpI)
+            {
+                if(!isConstrained[movingPointsP[mpI]])
                 {
                     constrainedPoints.append(movingPointsP[mpI]);
                     constrainedVelocity.append(pistonVel);
-//                constrainedVelocity.append(vector::zero);
                     constraintSize++;
+                    isConstrained[movingPointsP[mpI]] = true;
+                }
+
+            }
+        }
+
+        {
+
+            label headPtsIndex = pZones.findZoneID("headPoints");
+            const labelList& fixedPointsH = pZones[headPtsIndex];
+
+            forAll(fixedPointsH, fpI)
+            {
+                if(!isConstrained[fixedPointsH[fpI]])
+                {
+                    constrainedPoints.append(fixedPointsH[fpI]);
+                    constrainedVelocity.append(vector::zero);
+                    constraintSize++;
+                    isConstrained[fixedPointsH[fpI]] = true;
+                }
 
             }
         }
     }
 
-    Info << "constraintSize = " << constraintSize << endl;
+    Pout << "constraintSize = " << constraintSize << endl;
 }

src/engine/engineTopoChangerMesh/accordionEngineMesh/setAccordionEngineMeshConstraintsNoDeformation.H

@@ -1,72 +0,0 @@
-{
-
-    const faceZoneMesh& fZones = faceZones();
-    const pointZoneMesh& pZones = pointZones();
-
-    label constraintSize = 0;
-
-    forAll(valves_, valveI)
-    {
-    
-        vector valveVel =
-            valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
-
-        {
-    
-            label movingPtsIndex = pZones.findZoneID("movingPointsV"+Foam::name(valveI + 1));
-            const labelList& movingPointsV = pZones[movingPtsIndex];
-    
-            forAll(movingPointsV, mpI)
-            {
-                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]);
-                constrainedVelocity.append(vector::zero);
-                constraintSize++;
-            }
-
-        }
-
-        {
-    
-            label movingFacesIndex = fZones.findZoneID("movingFacesZoneV"+Foam::name(valveI + 1));
-            const labelList& movingFacesV = pZones[movingFacesIndex];
-    
-            forAll(movingFacesV, mfI)
-            {
-                constrainedPoints.append(movingFacesV[mfI]+points().size());
-                constrainedVelocity.append(valveVel);
-                constraintSize++;
-            }
-
-        }
-
-        {
-    
-            label staticFacesIndex = fZones.findZoneID("staticFacesZoneV"+Foam::name(valveI + 1));
-            const labelList& staticFacesV = pZones[staticFacesIndex];
-    
-            forAll(staticFacesV, sfI)
-            {
-                constrainedPoints.append(staticFacesV[sfI]+points().size());
-                constrainedVelocity.append(vector::zero);
-                constraintSize++;
-            }
-
-        }
-
-    }
-
-}

src/engine/engineTopoChangerMesh/accordionEngineMesh/setValveMotionBoundaryConditionAccordionEngineMesh.H

@@ -1,7 +1,6 @@
 // use tetrahedral decomposition of the engine mesh
 
 {
-    
     Info << "setting valve motion b.c." << endl;
 
     tetPointVectorField& motionU = mSolver.motionU();
@@ -9,8 +8,8 @@
     // Set valve velocity
     forAll (valves_, valveI)
     {
-        
-        Info << "Valve n. " << valveI + 1 << " velocity = " <<  valves_[valveI].curVelocity() << endl;
+        Info<< "Valve n. " << valveI + 1 << " velocity = "
+            <<  valves_[valveI].curVelocity() << endl;
 
         vector valveVel =
             valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
@@ -33,7 +32,6 @@
 
         if (valves_[valveI].poppetPatchID().active())
         {
-                            
             // Top of the valve does not move
             motionU.boundaryField()[valves_[valveI].poppetPatchID().index()] ==
                 valveVel;
@@ -45,82 +43,12 @@
             motionU.boundaryField()[valves_[valveI].sidePatchID().index()] ==
                 valveVel;
         }
-        
-        if
-        (
-	        moveDetach_ 
-	        && 
-            valves_[valveI].detachInPortPatchID().active()
-            && 
-            valves_[valveI].detachInCylinderPatchID().active()
-        )
-        {
-        
-            label detachInPortIndex = valves_[valveI].detachInPortPatchID().index();
-            label detachInCylIndex = valves_[valveI].detachInCylinderPatchID().index();
-
-            const coordinateSystem& cSysV = valves_[valveI].cs();
-                        
-            const pointField& detachInPortGlobal = 
-                    motionU.boundaryField()[detachInPortIndex].patch().localPoints();
-            
-            const pointField& detachInCylGlobal = 
-                    motionU.boundaryField()[detachInCylIndex].patch().localPoints();
-
-            scalarField detachInPortLocal =  cSysV.localPosition(detachInPortGlobal)().component(vector::Z);               
-            scalarField detachInCylLocal =  cSysV.localPosition(detachInCylGlobal)().component(vector::Z);               
-                
-            scalar maxDetachInPort = max(detachInPortLocal);
-            scalar minDetachInPort = min(detachInPortLocal);
-
-            scalar maxDetachInCyl = max(detachInCylLocal);
-            scalar minDetachInCyl = min(detachInCylLocal);
-                
-            Info << "maxDetachInPort = " << maxDetachInPort << endl;
-            Info << "minDetachInPort = " << minDetachInPort << endl;
-            Info << "maxDetachInCyl = " << maxDetachInCyl << endl;
-            Info << "minDetachInCyl = " << minDetachInCyl << endl;
-				
-				
-            motionU.boundaryField()[detachInPortIndex] == valveVel*(maxDetachInPort -
-            detachInPortLocal)/(maxDetachInPort-minDetachInPort);
-            
-            motionU.boundaryField()[detachInCylIndex] == valveVel*(maxDetachInCyl -
-            detachInCylLocal)/(maxDetachInCyl-minDetachInCyl);
-        }
-        else
-	    {
-/*        	if (valves_[valveI].detachInPortPatchID().active())
-        	{
-            		motionU.boundaryField()[valves_[valveI].detachInPortPatchID().index()] == vector::zero;
-        	}
-        	Info << "detach in port" << endl;
-
-        	if (valves_[valveI].detachInCylinderPatchID().active())
-        	{	
-            		motionU.boundaryField()[valves_[valveI].detachInCylinderPatchID().index()] == vector::zero;
-       	 	}
-        	Info << "detach in cyl" << endl;
-*/	        if
-	        (
-	            !moveDetach_
-  	            &&
-                valves_[valveI].detachInPortPatchID().active()
-                && 
-                valves_[valveI].detachInCylinderPatchID().active()
-	        )
-
-	        {
-        	    motionU.boundaryField()[valves_[valveI].detachInPortPatchID().index()] == vector::zero;
-        	    motionU.boundaryField()[valves_[valveI].detachInCylinderPatchID().index()] == vector::zero;
     }
 
-	    }
-    }
-//  Setting the boundary position
-
+    //  Setting the boundary position
     {
-        label cylinderHeadIndex = boundaryMesh().findPatchID(cylinderHeadName_);
+        label cylinderHeadIndex =
+            boundaryMesh().findPatchID(cylinderHeadName_);
 
         // Top of the valve does not move
         motionU.boundaryField()[cylinderHeadIndex] ==
@@ -128,5 +56,4 @@
     }
 
     Info << "valve motion boundary conditions set" << endl;
-    
 }

src/engine/engineTopoChangerMesh/attachDetachFunctions/attachDetachFunctions.C

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "attachDetachFunctions.H"
+#include "HashSet.H"
+
+using namespace Foam;
+
+// Find edge between points v0 and v1.
+label Foam::findEdge
+(
+    const primitiveMesh& mesh,
+    const label v0,
+    const label v1
+)
+{
+    const labelList& pEdges = mesh.pointEdges()[v0];
+
+    forAll(pEdges, pEdgeI)
+    {
+        label edgeI = pEdges[pEdgeI];
+
+        const edge& e = mesh.edges()[edgeI];
+
+        if (e.otherVertex(v0) == v1)
+        {
+            return edgeI;
+        }
+    }
+
+    FatalErrorIn
+    (
+        "findEdge(const primitiveMesh&, const label, const label)"
+    )   << "Cannot find edge between mesh points " << v0 << " and " << v1
+        << abort(FatalError);
+
+    return -1;
+}
+
+// Checks whether patch present
+void checkPatch
+(
+    const polyBoundaryMesh& bMesh,
+    const word& name
+)
+{
+    label patchI = bMesh.findPatchID(name);
+
+    if (patchI == -1)
+    {
+        FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
+            << "Cannot find patch " << name << endl
+            << "It should be present but of zero size" << endl
+            << "Valid patches are " << bMesh.names()
+            << exit(FatalError);
+    }
+
+    if (bMesh[patchI].size() != 0)
+    {
+        FatalErrorIn("checkPatch(const polyBoundaryMesh&, const word&)")
+            << "Patch " << name << " is present but not of zero size"
+            << exit(FatalError);
+    }
+}
+
+void changePatchID
+(
+    const polyMesh& mesh,
+    const label faceID,
+    const label patchID,
+    directTopoChange& meshMod
+)
+{
+    const label zoneID = mesh.faceZones().whichZone(faceID);
+
+    bool zoneFlip = false;
+
+    if (zoneID >= 0)
+    {
+        const faceZone& fZone = mesh.faceZones()[zoneID];
+
+        zoneFlip = fZone.flipMap()[fZone.whichFace(faceID)];
+    }
+
+    meshMod.setAction
+    (
+        polyModifyFace
+        (
+            mesh.faces()[faceID],               // face
+            faceID,                             // face ID
+            mesh.faceOwner()[faceID],           // owner
+            -1,                                 // neighbour
+            false,                              // flip flux
+            patchID,                            // patch ID
+            false,                              // remove from zone
+            zoneID,                             // zone ID
+            zoneFlip                            // zone flip
+        )
+    );
+}

src/engine/engineTopoChangerMesh/attachDetachFunctions/attachDetachFunctions.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    engineTopoChangerMesh
+
+Description
+    Base class for engine mesh with topological changes
+
+SourceFiles
+    engineTopoChangerMesh.C
+    newTopoFvMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef attachDetachFunctions_H
+#define attachDetachFunctions_H
+
+#include "argList.H"
+#include "polyMesh.H"
+#include "Time.H"
+#include "polyTopoChange.H"
+#include "mapPolyMesh.H"
+#include "faceSet.H"
+#include "attachDetach.H"
+#include "regionSide.H"
+#include "directTopoChange.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+    class nearestEqOp
+    {
+
+        public:
+
+        void operator()(vector& x, const vector& y) const
+        {
+            if (magSqr(y) < magSqr(x))
+            {
+                x = y;
+            }
+        }
+    };
+
+
+    // Find edge between points v0 and v1.
+    label findEdge
+    (
+       const primitiveMesh& mesh,
+       const label v0,
+       const label v1
+    );
+
+
+// Checks whether patch present
+    void checkPatch
+    (
+        const polyBoundaryMesh& bMesh,
+        const word& name
+    );
+
+    void changePatchID
+    (
+        const polyMesh& mesh,
+        const label faceID,
+        const label patchID,
+        directTopoChange& meshMod
+    );
+
+
+
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/deformingEngineMesh/addDeformingEngineMeshZones.C

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "deformingEngineMesh.H"
+#include "surfaceFields.H"
+#include "regionSplit.H"
+#include "directTopoChange.H"
+#include "zeroGradientTetPolyPatchFields.H"
+#include "tetDecompositionMotionSolver.H"
+
+#include "fixedValueTetPolyPatchFields.H"
+#include "mixedTetPolyPatchFields.H"
+#include "slipTetPolyPatchFields.H"
+#include "zeroGradientTetPolyPatchFields.H"
+#include "meshTools.H"
+#include "syncTools.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::deformingEngineMesh::addMeshZones()
+{
+    // Add the zones and mesh modifiers to operate piston motion
+    Switch addBoundary(true);
+
+    if
+    (
+        pointZones().size() > 0
+    )
+    {
+        Info<< "Time = " << engTime().theta() << endl;
+        Info<< "void Foam::deformingEngineMesh::addZonesAndModifiers() : "
+            << "Zones and modifiers already present.  Skipping."
+            << endl;
+
+        checkAndCalculate();
+
+        Info << "Point zones found = " << pointZones().size() << endl;
+        Info << "Face zones found = " << faceZones().size() << endl;
+        Info << "Cell zones found = " << cellZones().size() << endl;
+
+    }
+    else
+    {
+        pointZones().setSize(0);
+        faceZones().setSize(0);
+        cellZones().setSize(0);
+        topoChanger_.setSize(0);
+
+        Info << "checkAndCalculate()" << endl;
+        checkAndCalculate();
+
+        Info<< "Time = " << engTime().theta() << endl
+            << "Adding zones to the engine mesh" << endl;
+
+        DynamicList<pointZone*> pz;
+
+        DynamicList<faceZone*> fz;
+
+        DynamicList<cellZone*> cz;
+
+        label nPointZones = 0;
+        label nFaceZones = 0;
+        label nCellZones = 0;
+
+
+#       include "addPistonFacesPointZonesDeformingEngineMesh.H"
+
+#       include "addValveFaceZonesDeformingEngineMesh.H"
+
+#       include "addOutputCellsDeformingEngineMesh.H"
+
+
+        Info<< "Adding " << nPointZones << " point, "
+            << nFaceZones << " face zones and "
+            << nCellZones << " cell zones" << endl;
+
+        Info << boundary().size() << endl;
+
+        pz.setSize(nPointZones);
+        Info << "setSize pz" << endl;
+        fz.setSize(nFaceZones);
+        Info << "setSize fz" << endl;
+        cz.setSize(nCellZones);
+        Info << "setSize cz" << endl;
+
+        addZones(pz, fz, cz);
+    }
+
+
+    // Calculating the virtual positions of piston and valves
+
+ //   setVirtualPositions();
+
+    // Write mesh and modifiers
+    write();
+    Info << "virtualPistonPosition = " << virtualPistonPosition() << endl;
+    Info << "piston position = " << pistonPosition() << endl;
+}
+
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/deformingEngineMesh/addOutputCellsDeformingEngineMesh.H

@@ -0,0 +1,55 @@
+    {
+        word cylinderOutputName(engTime().engineDict().lookup("cylinderSetName"));
+
+        IOobject outputHeader
+        (
+            cylinderOutputName,
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+
+        if(!outputHeader.headerOk())
+        {
+            WarningIn
+            (
+                "deformingEngineMesh::addZonesAndModifiers()"
+            )   << "faceSet called " << cylinderOutputName
+                << " does not exist. Continuing mesh motion without cylinder output cells" << endl;
+
+        cz.append
+        (
+          new cellZone
+          (
+                "cylinderOutputCells",
+                labelList(0),
+                nCellZones,
+                cellZones()
+          )
+      );
+
+      nCellZones++;
+
+
+        }
+        else
+    {
+
+        cellSet outputCellSet(*this, cylinderOutputName);
+
+        labelList outputCellSetList = outputCellSet.toc();
+
+        cz.append
+        (
+          new cellZone
+          (
+                "cylinderOutputCells",
+                outputCellSetList,
+                nCellZones,
+                cellZones()
+          )
+      );
+
+      nCellZones++;
+    }
+}

src/engine/engineTopoChangerMesh/deformingEngineMesh/addPistonFacesPointZonesDeformingEngineMesh.H

@@ -0,0 +1,186 @@
+
+    // Add the piston points zone
+    if (piston().patchID().active())
+    {
+        IOobject pistonCellsHeader
+        (
+            piston().pistonCellSetName(),
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+
+        if(!pistonCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "deformingEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << piston().pistonCellSetName()
+                << " does not exist. Continuing mesh motion without piston points" << endl;
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "movingPistonPoints",
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+
+        }
+        else
+        {
+            cellSet movingPistonCells(*this, piston().pistonCellSetName());
+
+            labelList pistonCells = movingPistonCells.toc();
+
+            labelList movingPistonPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (pistonCells, cellI)
+            {
+                const labelList& curCellPoints = cp[pistonCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }
+
+            movingPistonPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    movingPistonPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "movingPistonPoints",
+                    movingPistonPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+
+        }
+
+    }
+
+    {
+
+        IOobject headCellsHeader
+        (
+            headCellSetName_,
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+
+        if(!headCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "deformingEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << headCellSetName_
+                << " does not exist. Continuing mesh motion without head points" << endl;
+            pz.append
+            (
+                new pointZone
+                (
+                    "headPoints",
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+        }
+        else
+        {
+            cellSet fixedHeadCells(*this, headCellSetName_);
+
+            labelList headCells = fixedHeadCells.toc();
+
+            labelList fixedHeadPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (fixedHeadCells, cellI)
+            {
+                const labelList& curCellPoints = cp[headCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }
+
+            fixedHeadPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    fixedHeadPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "headPoints",
+                    fixedHeadPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+        }
+
+    }

src/engine/engineTopoChangerMesh/deformingEngineMesh/addValveFaceZonesDeformingEngineMesh.H

@@ -0,0 +1,190 @@
+{
+    for (label valveI = 0; valveI < nValves(); valveI++)
+    {
+
+        // for each valve the following zones have to be created:
+        //
+        //  - movingFaceZone
+        //  - movingPointZone
+        //  - staticFaceZone
+        //  - staticPointZone
+
+        IOobject movingCellsHeader
+        (
+            valves_[valveI].movingCellsName(),
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+
+        if(!movingCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "deformingEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << valves_[valveI].movingCellsName()
+                << " does not exist. Continuing mesh motion without rigid motion points for valve " << valves_[valveI].name() << endl;
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "movingPointsV" + Foam::name(valveI + 1),
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+
+        }
+        else
+        {
+            cellSet movingCellsSet(*this, valves_[valveI].movingCellsName());
+
+            labelList movingCells = movingCellsSet.toc();
+
+            labelList movingPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (movingCells, cellI)
+            {
+                const labelList& curCellPoints = cp[movingCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }
+
+            movingPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    movingPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "movingPointsV" + Foam::name(valveI + 1),
+                    movingPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+        }
+
+        IOobject staticCellsHeader
+        (
+            valves_[valveI].staticCellsName(),
+            time().constant()+"/polyMesh/sets/",
+            *this,
+            IOobject::MUST_READ
+        );
+
+        if(!staticCellsHeader.headerOk())
+        {
+            WarningIn
+            (
+                "deformingEngineMesh::addZonesAndModifiers()"
+            )   << "cellSet called " << valves_[valveI].staticCellsName()
+                << " does not exist. Continuing mesh motion without static points for valve " << valves_[valveI].name() << endl;
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "staticPointsV" + Foam::name(valveI + 1),
+                    labelList(0),
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+
+        }
+        else
+        {
+            cellSet staticCellSet(*this, valves_[valveI].staticCellsName());
+
+            labelList staticCells = staticCellSet.toc();
+
+            labelList staticPoints;
+
+            const labelListList& cp = cellPoints();
+
+            boolList count(allPoints().size(), false);
+
+            forAll (staticCells, cellI)
+            {
+                const labelList& curCellPoints = cp[staticCells[cellI]];
+
+                forAll (curCellPoints, i)
+                {
+                    count[curCellPoints[i]] = true;
+                }
+            }
+
+            // Count the points
+            label nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    nCounted++;
+                }
+            }
+
+            staticPoints.setSize(nCounted);
+
+            // Collect the points
+            nCounted = 0;
+            forAll (count, pointI)
+            {
+                if (count[pointI] == true)
+                {
+                    staticPoints[nCounted] = pointI;
+                    nCounted++;
+                }
+            }
+
+            pz.append
+            (
+                new pointZone
+                (
+                    "staticPointsV" + Foam::name(valveI + 1),
+                    staticPoints,
+                    nPointZones,
+                    pointZones()
+                )
+            );
+
+            nPointZones++;
+        }
+    }
+}

src/engine/engineTopoChangerMesh/deformingEngineMesh/deformingEngineMesh.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "deformingEngineMesh.H"
+#include "componentMixedTetPolyPatchVectorField.H"
+#include "mapPolyMesh.H"
+#include "polyTopoChange.H"
+#include "addToRunTimeSelectionTable.H"
+#include "GeometricField.H"
+#include "volMesh.H"
+#include "engineTime.H"
+#include "pointField.H"
+#include "fvPatchField.H"
+#include "Switch.H"
+#include "symmetryFvPatch.H"
+#include "tetDecompositionMotionSolver.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(deformingEngineMesh, 0);
+    addToRunTimeSelectionTable(engineTopoChangerMesh, deformingEngineMesh, IOobject);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::deformingEngineMesh::realDeformation() const
+{
+
+    if (virtualPistonPosition() + engTime().pistonDisplacement().value() > deckHeight_ - SMALL)
+    {
+        return true;
+    }
+    else
+    {
+        return deformation();
+    }
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::deformingEngineMesh::deformingEngineMesh
+(
+    const IOobject& io
+)
+:
+    engineTopoChangerMesh(io),
+    piston_(*this, engTime().engineDict().subDict("piston")),
+    valves_(*this, engTime().engineDict().lookup("deformingEngineMesh")),
+    pistonPosition_(-GREAT),
+    virtualPistonPosition_(-GREAT),
+    deckHeight_(GREAT),
+    cylinderHeadName_(engTime().engineDict().lookup("cylinderHeadName")),
+    linerName_(engTime().engineDict().lookup("linerName")),
+    headCellSetName_(engTime().engineDict().lookup("headCellSetName"))
+
+{
+    // Add zones and modifiers if not already there.
+    addMeshZones();
+    msPtr_ = motionSolver::New(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::deformingEngineMesh::setBoundaryVelocity(volVectorField& U)
+{
+    // Set valve velociaty
+    forAll (valves(), valveI)
+    {
+        vector valveVel =
+            valves()[valveI].curVelocity()*valves()[valveI].cs().axis();
+
+
+        // If valve is present in geometry, set the motion
+        if (valves()[valveI].stemPatchID().active())
+        {
+            // Bottom of the valve moves with given velocity
+            U.boundaryField()[valves()[valveI].stemPatchID().index()] ==
+                valveVel;
+        }
+
+    }
+
+}
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/deformingEngineMesh/deformingEngineMesh.H

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    deformingEngineMesh
+
+DESCRIPTION
+
+*/
+// ------------------------------------------------------------------------- //
+
+#ifndef deformingEngineMesh_H
+#define deformingEngineMesh_H
+
+#include "engineTopoChangerMesh.H"
+#include "enginePiston.H"
+#include "motionSolver.H"
+#include "polyPatchID.H"
+#include "accordionValveBank.H"
+#include "twoDPointCorrector.H"
+
+#include "faceSet.H"
+#include "pointSet.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+
+/*---------------------------------------------------------------------------*\
+                           Class lateralValves Declaration
+\*---------------------------------------------------------------------------*/
+class twoDPointCorrector;
+
+class deformingEngineMesh
+:
+    public engineTopoChangerMesh
+{
+    // Private data
+
+        //- Piston patch
+        enginePiston piston_;
+
+        //- Lateral valves
+        accordionValveBank valves_;
+
+
+        //- Piston Position
+        scalar pistonPosition_;
+
+        //- Virtual piston position (pistonPosition + offSet)
+        scalar virtualPistonPosition_;
+
+        //- deckHeight
+        scalar deckHeight_;
+
+        //- Motion solver
+        autoPtr<motionSolver> msPtr_;
+
+        //- Cylinder head name
+        word cylinderHeadName_;
+
+        //- Cylinder head name
+        word linerName_;
+
+        //- PistonFaceSet
+//        word pistonAuxPoints_;
+
+        //- PistonFaceSet
+        word headCellSetName_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        deformingEngineMesh(const deformingEngineMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const deformingEngineMesh&);
+
+        //- Check if all patches exist, then calculate virtualPistonPosition,
+        //- pistonPosition and deckHeight for the first time
+        void checkAndCalculate();
+
+        //- Calculate the virtualPistonPosition,
+        void setVirtualPositions();
+
+        bool realDeformation() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("deformingEngineMesh");
+
+
+    // Constructors
+
+        //- Construct from database
+        explicit deformingEngineMesh(const IOobject& io);
+
+
+    // Destructor - default
+
+
+    // Member Functions
+
+        inline const enginePiston& piston() const;
+        inline const accordionValveBank& valves() const;
+        inline const scalar& pistonPosition() const;
+        inline scalar& pistonPosition();
+        inline const scalar& virtualPistonPosition() const;
+        inline scalar& virtualPistonPosition();
+        inline const scalar& deckHeight() const;
+        inline scalar& deckHeight();
+
+
+        //- Return true for mesh deformation mode
+        bool deformation() const
+        {
+            return true;
+        }
+
+        //- Return number of valves
+        label nValves() const
+        {
+            return valves_.size();
+        }
+
+        //- Add valve and piston zones and modifiers
+        void addMeshZones();
+
+        //- Move and morph
+        virtual bool update();
+
+        //- Set boundary velocities
+        void setBoundaryVelocity(volVectorField& U);
+
+        //- Generate the mesh from a series of square blocks
+        void generateMesh(){};
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "deformingEngineMeshI.H"
+
+#endif
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshI.H

@@ -22,15 +22,60 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
-#include "accordionEngineMesh.H"
-#include "slidingInterface.H"
-#include "layerAdditionRemoval.H"
-#include "surfaceFields.H"
-#include "regionSplit.H"
-#include "componentMixedTetPolyPatchVectorField.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+namespace Foam
+{
+
+inline const enginePiston& deformingEngineMesh::piston() const
+{
+    return piston_;
+}
 
 
+inline const accordionValveBank& deformingEngineMesh::valves() const
+{
+    return valves_;
+}
+
+
+inline const scalar& deformingEngineMesh::pistonPosition() const
+{
+    return pistonPosition_;
+}
+
+inline scalar& deformingEngineMesh::pistonPosition()
+{
+    return pistonPosition_;
+}
+
+inline const scalar& deformingEngineMesh::virtualPistonPosition() const
+{
+    return virtualPistonPosition_;
+}
+
+inline scalar& deformingEngineMesh::virtualPistonPosition()
+{
+    return virtualPistonPosition_;
+}
+
+
+inline const scalar& deformingEngineMesh::deckHeight() const
+{
+    return deckHeight_;
+}
+
+inline scalar& deformingEngineMesh::deckHeight()
+{
+    return deckHeight_;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshInitialize.C

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "deformingEngineMesh.H"
+#include "surfaceFields.H"
+#include "regionSplit.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::deformingEngineMesh::checkAndCalculate()
+{
+
+    label pistonIndex = -1;
+    bool foundPiston = false;
+
+    label linerIndex = -1;
+    bool foundLiner = false;
+
+    label cylinderHeadIndex = -1;
+    bool foundCylinderHead = false;
+
+
+    forAll(boundary(), i)
+    {
+        if (boundary()[i].name() == piston().patchID().name())
+        {
+            pistonIndex = i;
+            foundPiston = true;
+        }
+        else if (boundary()[i].name() == linerName_)
+        {
+            linerIndex = i;
+            foundLiner = true;
+        }
+        else if (boundary()[i].name() == cylinderHeadName_)
+        {
+            cylinderHeadIndex = i;
+            foundCylinderHead = true;
+        }
+    }
+
+    reduce(foundPiston, orOp<bool>());
+    reduce(foundLiner, orOp<bool>());
+    reduce(foundCylinderHead, orOp<bool>());
+
+    if (!foundPiston)
+    {
+        FatalErrorIn("Foam::deformingEngineMesh::checkAndCalculate()")
+            << " : cannot find piston patch"
+            << abort(FatalError);
+    }
+
+    if (!foundLiner)
+    {
+        FatalErrorIn("Foam::deformingEngineMesh::checkAndCalculate()")
+            << " : cannot find liner patch"
+            << abort(FatalError);
+    }
+
+    if (!foundCylinderHead)
+    {
+        FatalErrorIn("Foam::deformingEngineMesh::checkAndCalculate()")
+            << " : cannot find cylinderHead patch"
+            << exit(FatalError);
+    }
+
+    {
+        if (linerIndex != -1)
+        {
+            pistonPosition() =
+                max(boundary()[pistonIndex].patch().localPoints()).z();
+        }
+        reduce(pistonPosition(), minOp<scalar>());
+
+        if (cylinderHeadIndex != -1)
+        {
+            deckHeight() = min
+            (
+                boundary()[cylinderHeadIndex].patch().localPoints()
+            ).z();
+
+        }
+        reduce(deckHeight(), minOp<scalar>());
+
+        Info<< "deckHeight: " << deckHeight() << nl
+            << "piston position: " << pistonPosition() << endl;
+    }
+
+
+}
+
+void Foam::deformingEngineMesh::setVirtualPositions()
+{
+}

src/engine/engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshMove.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "deformingEngineMesh.H"
+#include "surfaceFields.H"
+#include "regionSplit.H"
+#include "componentMixedTetPolyPatchVectorField.H"
+#include "mapPolyMesh.H"
+
+#include "tetPolyMesh.H"
+#include "tetPointFields.H"
+#include "elementFields.H"
+#include "fixedValueTetPolyPatchFields.H"
+#include "slipTetPolyPatchFields.H"
+#include "symmetryTetPolyPatch.H"
+
+#include "tetFem.H"
+
+#include "symmetryFvPatch.H"
+#include "wedgeFvPatch.H"
+#include "emptyFvPatch.H"
+#include "zeroGradientTetPolyPatchFields.H"
+#include "tetDecompositionMotionSolver.H"
+
+#include "fixedValueTetPolyPatchFields.H"
+#include "mixedTetPolyPatchFields.H"
+#include "slipTetPolyPatchFields.H"
+#include "zeroGradientTetPolyPatchFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::deformingEngineMesh::update()
+{
+    tetDecompositionMotionSolver& mSolver =
+        refCast<tetDecompositionMotionSolver>(msPtr_());
+
+    scalar deltaZ = engTime().pistonDisplacement().value();
+
+    // deltaZ set to zero, FIXED PISTON POSITION
+    deltaZ = 0.0;
+
+    virtualPistonPosition() += deltaZ;
+
+    pointField newPoints = points();
+
+    {
+#       include "setValveMotionBoundaryConditionDeformingEngineMesh.H"
+#       include "setPistonMotionBoundaryConditionDeformingEngineMesh.H"
+        Info << "piston motion" << endl;
+
+        DynamicList<label> constrainedPoints(mSolver.curPoints()().size()/100);
+        DynamicList<vector> constrainedVelocity
+        (
+            mSolver.curPoints()().size()/100
+        );
+
+#       include "setDeformingEngineMeshConstraints.H"
+
+
+        labelList constrainedPointsList(constrainedPoints.shrink());
+        vectorField constrainedVelocityField(constrainedVelocity.shrink());
+
+        forAll (constrainedPointsList, i)
+        {
+            mSolver.setConstraint
+            (
+                constrainedPointsList[i],
+                constrainedVelocityField[i]
+            );
+        }
+
+        mSolver.solve();
+
+        newPoints = mSolver.curPoints();
+        movePoints(newPoints);
+        setVirtualPositions();
+        mSolver.clearConstraints();
+
+    }
+
+    pointField oldPointsNew = oldAllPoints();
+    newPoints = points();
+    movePoints(newPoints);
+    Info << "max mesh phi before = " << max((phi())) << endl;
+    Info << "min mesh phi before = " << min((phi())) << endl;
+
+
+    Info << "max mesh phi after = " << max((phi())) << endl;
+    Info << "min mesh phi after = " << min((phi())) << endl;
+
+    return false;
+}

src/engine/engineTopoChangerMesh/deformingEngineMesh/movePistonPointsDeformingEngineMesh.H

@@ -0,0 +1,96 @@
+    Info << "virtualPistonPosition = " << virtualPistonPosition()
+    << ", deckHeight = " << deckHeight() << endl;
+
+    // Mesh in three parts:
+    // - pistonPoints - move with deltaZ
+    // - headPoints - do not move
+
+    const pointZoneMesh& pZones = pointZones();
+    label headPtsIndex = pZones.findZoneID("headPoints");
+    label pistonPtsIndex = pZones.findZoneID("pistonPoints");
+    const labelList& pistonPoints = pZones[pistonPtsIndex];
+    const labelList& headPoints = pZones[headPtsIndex];
+
+    // Whether point displacement is by scaling
+    boolList scaleDisp(nPoints(), true);
+    label nScaled = nPoints();
+    List<bool> pistonPoint(newPoints.size(), false);
+    List<bool> headPoint(newPoints.size(), false);
+
+    forAll(pistonPoints, i)
+    {
+        label pointI = pistonPoints[i];
+        pistonPoint[pointI] = true;
+        point& p = newPoints[pointI];
+
+        if (p.z() < pistonPosition() - 1.0e-6)
+        {
+            scaleDisp[pointI] = false;
+            nScaled--;
+        }
+    }
+
+    forAll(headPoints, i)
+    {
+        headPoint[headPoints[i]] = true;
+        scaleDisp[headPoints[i]] = false;
+        nScaled--;
+    }
+
+    if (realDeformation)
+    {
+
+        forAll(scaleDisp, pointI)
+        {
+            point& p = newPoints[pointI];
+
+            if (scaleDisp[pointI])
+            {
+                p.z() +=
+                    deltaZ
+                  * (deckHeight() - p.z())/(deckHeight() - pistonPosition());
+            }
+            else
+            {
+                if (!headPoint[pointI])
+                {
+                    p.z() += deltaZ;
+                }
+            }
+        }
+    }
+    else
+    {
+
+        // Always move piston
+        scalar pistonTopZ = -GREAT;
+        forAll(pistonPoints, i)
+        {
+            point& p = newPoints[pistonPoints[i]];
+            p.z() += deltaZ;
+            pistonTopZ = max(pistonTopZ, p.z());
+        }
+
+        if (deltaZ > 0.0)
+        {
+            // check if piston-points have moved beyond the layer above
+            forAll(newPoints, i)
+            {
+                if (!pistonPoint[i])
+                {
+                    if (virtualPistonPosition() > newPoints[i].z())
+                    {
+                        newPoints[i].z() = pistonTopZ + 0.9*minLayerThickness;
+                    }
+                }
+            }
+        }
+    }
+
+    movePoints(newPoints);
+
+    pistonPosition() += deltaZ;
+    scalar pistonSpeed = deltaZ/engTime().deltaT().value();
+
+    Info<< "clearance: " << deckHeight() - pistonPosition() << nl
+        << "Piston speed = " << pistonSpeed << " m/s" << endl;

src/engine/engineTopoChangerMesh/deformingEngineMesh/setDeformingEngineMeshConstraints.H

@@ -0,0 +1,90 @@
+{
+    boolList isConstrained(points().size(), false);
+
+    const pointZoneMesh& pZones = pointZones();
+
+    label constraintSize = 0;
+
+    forAll(valves_, valveI)
+    {
+        vector valveVel =
+            valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
+
+        {
+            label movingPtsIndex =
+                pZones.findZoneID("movingPointsV"+Foam::name(valveI + 1));
+            const labelList& movingPointsV = pZones[movingPtsIndex];
+
+            forAll(movingPointsV, mpI)
+            {
+                if(!isConstrained[movingPointsV[mpI]])
+                {
+                    constrainedPoints.append(movingPointsV[mpI]);
+                    constrainedVelocity.append(valveVel);
+
+                    constraintSize++;
+                    isConstrained[movingPointsV[mpI]] = true;
+                }
+            }
+        }
+
+        {
+
+            label staticPtsIndex =
+                pZones.findZoneID("staticPointsV" + Foam::name(valveI + 1));
+            const labelList& staticPointsV = pZones[staticPtsIndex];
+
+            forAll(staticPointsV, spI)
+            {
+                if(!isConstrained[staticPointsV[spI]])
+                {
+                    constrainedPoints.append(staticPointsV[spI]);
+                    constrainedVelocity.append(vector::zero);
+                    constraintSize++;
+                    isConstrained[staticPointsV[spI]] = true;
+                }
+            }
+
+        }
+
+        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)
+            {
+                if(!isConstrained[movingPointsP[mpI]])
+                {
+                    constrainedPoints.append(movingPointsP[mpI]);
+                    constrainedVelocity.append(pistonVel);
+                    constraintSize++;
+                    isConstrained[movingPointsP[mpI]] = true;
+                }
+
+            }
+        }
+
+        {
+
+            label headPtsIndex = pZones.findZoneID("headPoints");
+            const labelList& fixedPointsH = pZones[headPtsIndex];
+
+            forAll(fixedPointsH, fpI)
+            {
+                if(!isConstrained[fixedPointsH[fpI]])
+                {
+                    constrainedPoints.append(fixedPointsH[fpI]);
+                    constrainedVelocity.append(vector::zero);
+                    constraintSize++;
+                    isConstrained[fixedPointsH[fpI]] = true;
+                }
+
+            }
+        }
+    }
+
+}

src/engine/engineTopoChangerMesh/deformingEngineMesh/setPistonMotionBoundaryConditionDeformingEngineMesh.H

@@ -0,0 +1,22 @@
+
+    vector pistonVel =
+            piston().cs().axis()*engineTime_.pistonSpeed().value();
+
+{
+    tetPointVectorField& motionU = mSolver.motionU();
+
+    Info << "setting the piston velocity" << endl;
+
+    if (piston().patchID().active())
+    {
+
+        if (debug)
+        {
+            Info << "Piston velocity: " << pistonVel;
+        }
+
+        motionU.boundaryField()[piston().patchID().index()] == pistonVel;
+
+    }
+
+}

src/engine/engineTopoChangerMesh/deformingEngineMesh/setValveMotionBoundaryConditionDeformingEngineMesh.H

@@ -0,0 +1,62 @@
+// use tetrahedral decomposition of the engine mesh
+
+{
+
+    Info << "setting valve motion b.c." << endl;
+
+    tetPointVectorField& motionU = mSolver.motionU();
+
+    // Set valve velocity
+    forAll (valves_, valveI)
+    {
+
+        Info << "Valve n. " << valveI + 1 << " velocity = " <<  valves_[valveI].curVelocity() << endl;
+
+        vector valveVel =
+            valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
+
+        // If valve is present in geometry, set the motion
+        if (valves_[valveI].bottomPatchID().active())
+        {
+            // Bottom of the valve moves with given velocity
+            motionU.boundaryField()[valves_[valveI].bottomPatchID().index()] ==
+                valveVel;
+
+            if (debug)
+            {
+                Info<< "Valve " << valveI << " lift: "
+                    << valves_[valveI].curLift()
+                    << " velocity: " << valves_[valveI].curVelocity()
+                    << endl;
+            }
+        }
+
+        if (valves_[valveI].poppetPatchID().active())
+        {
+
+            // Top of the valve does not move
+            motionU.boundaryField()[valves_[valveI].poppetPatchID().index()] ==
+                valveVel;
+        }
+
+        if (valves_[valveI].sidePatchID().active())
+        {
+            // Top of the valve does not move
+            motionU.boundaryField()[valves_[valveI].sidePatchID().index()] ==
+                valveVel;
+        }
+
+    }
+//  Setting the boundary position
+
+    {
+        label cylinderHeadIndex = boundaryMesh().findPatchID(cylinderHeadName_);
+
+        // Top of the valve does not move
+        motionU.boundaryField()[cylinderHeadIndex] ==
+            vector::zero;
+    }
+
+    Info << "valve motion boundary conditions set" << endl;
+
+}

src/engine/engineTopoChangerMesh/engineValveSliding/setEngineValveSlidingConstraint.H

@@ -9,7 +9,6 @@
 
     forAll(valves_, valveI)
     {
-    
         vector valveVel =
             valves_[valveI].curVelocity()*valves_[valveI].cs().axis();
 
@@ -18,7 +17,6 @@
             valveVel = vector::zero;
         }
 
-    
         {
 
 //            label movingPtsIndex = pZones.findZoneID("movingPointsV"+Foam::name(valveI + 1));
@@ -79,11 +77,9 @@
 
                 constraintSize++;
             }
-    
         }
 
         {
-    
 //            label staticPtsIndex = pZones.findZoneID("staticPointsV"+Foam::name(valveI + 1));
 //            const labelList& staticPointsV = pZones[staticPtsIndex];
 /*
@@ -203,7 +199,6 @@
             }
         }
 
-        
         forAll(movingPointsP, mpI)
         {
             constrainedPoints.append(movingPointsP[mpI]);
@@ -213,6 +208,4 @@
 
 
     }
-    
-
 }

src/engine/engineTopoChangerMesh/regionSide/regionSide.C

@@ -0,0 +1,483 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+\*---------------------------------------------------------------------------*/
+
+#include "regionSide.H"
+#include "meshTools.H"
+#include "primitiveMesh.H"
+#include "IndirectList.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(regionSide, 0);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Step across edge onto other face on cell
+Foam::label Foam::regionSide::otherFace
+(
+    const primitiveMesh& mesh,
+    const label cellI,
+    const label faceI,
+    const label edgeI
+)
+{
+    label f0I, f1I;
+
+    meshTools::getEdgeFaces(mesh, cellI, edgeI, f0I, f1I);
+
+    if (f0I == faceI)
+    {
+        return f1I;
+    }
+    else
+    {
+        return f0I;
+    }
+}
+
+
+// Step across point to other edge on face
+Foam::label Foam::regionSide::otherEdge
+(
+    const primitiveMesh& mesh,
+    const label faceI,
+    const label edgeI,
+    const label pointI
+)
+{
+    const edge& e = mesh.edges()[edgeI];
+
+    // Get other point on edge.
+    label freePointI = e.otherVertex(pointI);
+
+    const labelList& fEdges = mesh.faceEdges()[faceI];
+
+    forAll(fEdges, fEdgeI)
+    {
+        label otherEdgeI = fEdges[fEdgeI];
+
+        const edge& otherE = mesh.edges()[otherEdgeI];
+
+        if
+        (
+            (
+                otherE.start() == pointI
+             && otherE.end() != freePointI
+            )
+         || (
+                otherE.end() == pointI
+             && otherE.start() != freePointI
+            )
+        )
+        {
+            // otherE shares one (but not two) points with e.
+            return otherEdgeI;
+        }
+    }
+
+    FatalErrorIn
+    (
+        "regionSide::otherEdge(const primitiveMesh&, const label, const label"
+        ", const label)"
+    )   << "Cannot find other edge on face " << faceI << " that uses point "
+        << pointI << " but not point " << freePointI << endl
+        << "Edges on face:" << fEdges
+        << " verts:" << IndirectList<edge>(mesh.edges(), fEdges)()
+        << " Vertices on face:"
+        << mesh.faces()[faceI]
+        << " Vertices on original edge:" << e << abort(FatalError);
+
+    return -1;
+}
+
+
+// Step from faceI (on side cellI) to connected face & cell without crossing
+// fenceEdges.
+void Foam::regionSide::visitConnectedFaces
+(
+    const primitiveMesh& mesh,
+    const labelHashSet& region,
+    const labelHashSet& fenceEdges,
+    const label cellI,
+    const label faceI,
+    labelHashSet& visitedFace
+)
+{
+    if (!visitedFace.found(faceI))
+    {
+        if (debug)
+        {
+            Info<< "visitConnectedFaces : cellI:" << cellI << " faceI:"
+                << faceI << "  isOwner:" << (cellI == mesh.faceOwner()[faceI])
+                << endl;
+        }
+
+        // Mark as visited
+        visitedFace.insert(faceI);
+
+        // Mark which side of face was visited.
+        if (cellI == mesh.faceOwner()[faceI])
+        {
+            sideOwner_.insert(faceI);
+        }
+
+
+        // Visit all neighbouring faces on faceSet. Stay on this 'side' of
+        // face by doing edge-face-cell walk.
+        const labelList& fEdges = mesh.faceEdges()[faceI];
+
+        forAll(fEdges, fEdgeI)
+        {
+            label edgeI = fEdges[fEdgeI];
+
+            if (!fenceEdges.found(edgeI))
+            {
+                // Step along faces on edge from cell to cell until
+                // we hit face on faceSet.
+
+                // Find face reachable from edge
+                label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
+
+                if (mesh.isInternalFace(otherFaceI))
+                {
+                    label otherCellI = cellI;
+
+                    // Keep on crossing faces/cells until back on face on
+                    // surface
+                    while (!region.found(otherFaceI))
+                    {
+                        visitedFace.insert(otherFaceI);
+
+                        if (debug)
+                        {
+                            Info<< "visitConnectedFaces : cellI:" << cellI
+                                << " found insideEdgeFace:" << otherFaceI
+                                << endl;
+                        }
+
+
+                        // Cross otherFaceI into neighbouring cell
+                        otherCellI =
+                            meshTools::otherCell
+                            (
+                                mesh,
+                                otherCellI,
+                                otherFaceI
+                            );
+
+                        otherFaceI =
+                                otherFace
+                                (
+                                    mesh,
+                                    otherCellI,
+                                    otherFaceI,
+                                    edgeI
+                                );
+                    }
+
+                    visitConnectedFaces
+                    (
+                        mesh,
+                        region,
+                        fenceEdges,
+                        otherCellI,
+                        otherFaceI,
+                        visitedFace
+                    );
+                }
+            }
+        }
+    }
+}
+
+
+// From edge on face connected to point on region (regionPointI) cross
+// to all other edges using this point by walking across faces
+// Does not cross regionEdges so stays on one side
+// of region
+void Foam::regionSide::walkPointConnectedFaces
+(
+    const primitiveMesh& mesh,
+    const labelHashSet& regionEdges,
+    const label regionPointI,
+    const label startFaceI,
+    const label startEdgeI,
+    labelHashSet& visitedEdges
+)
+{
+    // Mark as visited
+    insidePointFaces_.insert(startFaceI);
+
+    if (debug)
+    {
+        Info<< "walkPointConnectedFaces : regionPointI:" << regionPointI
+            << " faceI:" << startFaceI
+            << " edgeI:" << startEdgeI << " verts:"
+            << mesh.edges()[startEdgeI]
+            << endl;
+    }
+
+    // Cross faceI i.e. get edge not startEdgeI which uses regionPointI
+    label edgeI = otherEdge(mesh, startFaceI, startEdgeI, regionPointI);
+
+    if (!regionEdges.found(edgeI))
+    {
+        if (!visitedEdges.found(edgeI))
+        {
+            visitedEdges.insert(edgeI);
+
+            if (debug)
+            {
+                Info<< "Crossed face from "
+                    << " edgeI:" << startEdgeI << " verts:"
+                    << mesh.edges()[startEdgeI]
+                    << " to edge:" << edgeI << " verts:"
+                    << mesh.edges()[edgeI]
+                    << endl;
+            }
+
+            // Cross edge to all faces connected to it.
+
+            const labelList& eFaces = mesh.edgeFaces()[edgeI];
+
+            forAll(eFaces, eFaceI)
+            {
+                label faceI = eFaces[eFaceI];
+
+                walkPointConnectedFaces
+                (
+                    mesh,
+                    regionEdges,
+                    regionPointI,
+                    faceI,
+                    edgeI,
+                    visitedEdges
+                );
+            }
+        }
+    }
+}
+
+
+// Find all faces reachable from all non-fence points and staying on
+// regionFaces side.
+void Foam::regionSide::walkAllPointConnectedFaces
+(
+    const primitiveMesh& mesh,
+    const labelHashSet& regionFaces,
+    const labelHashSet& fencePoints
+)
+{
+    //
+    // Get all (internal and external) edges on region.
+    //
+    labelHashSet regionEdges(4*regionFaces.size());
+
+    for
+    (
+        labelHashSet::const_iterator iter = regionFaces.begin();
+        iter != regionFaces.end();
+        ++iter
+    )
+    {
+        label faceI = iter.key();
+
+        const labelList& fEdges = mesh.faceEdges()[faceI];
+
+        forAll(fEdges, fEdgeI)
+        {
+            regionEdges.insert(fEdges[fEdgeI]);
+        }
+    }
+
+
+    //
+    // Visit all internal points on surface.
+    //
+
+    // Storage for visited points
+    labelHashSet visitedPoint(4*regionFaces.size());
+
+    // Insert fence points so we don't visit them
+    for
+    (
+        labelHashSet::const_iterator iter = fencePoints.begin();
+        iter != fencePoints.end();
+        ++iter
+    )
+    {
+        visitedPoint.insert(iter.key());
+    }
+
+    labelHashSet visitedEdges(2*fencePoints.size());
+
+
+    if (debug)
+    {
+        Info<< "Excluding visit of points:" << visitedPoint << endl;
+    }
+
+    for
+    (
+        labelHashSet::const_iterator iter = regionFaces.begin();
+        iter != regionFaces.end();
+        ++iter
+    )
+    {
+        label faceI = iter.key();
+
+        // Get side of face.
+        label cellI;
+
+        if (sideOwner_.found(faceI))
+        {
+            cellI = mesh.faceOwner()[faceI];
+        }
+        else
+        {
+            cellI = mesh.faceNeighbour()[faceI];
+        }
+
+        // Find starting point and edge on face.
+        const labelList& fEdges = mesh.faceEdges()[faceI];
+
+        forAll(fEdges, fEdgeI)
+        {
+            label edgeI = fEdges[fEdgeI];
+
+            // Get the face 'perpendicular' to faceI on region.
+            label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);
+
+            // Edge 
+            const edge& e = mesh.edges()[edgeI];
+
+            if (!visitedPoint.found(e.start()))
+            {
+                Info<< "Determining visibility from point " << e.start()
+                    << endl;
+
+                visitedPoint.insert(e.start());
+
+                //edgeI = otherEdge(mesh, otherFaceI, edgeI, e.start());
+
+                walkPointConnectedFaces
+                (
+                    mesh,
+                    regionEdges,
+                    e.start(),
+                    otherFaceI,
+                    edgeI,
+                    visitedEdges
+                );
+            }
+            if (!visitedPoint.found(e.end()))
+            {
+                Info<< "Determining visibility from point " << e.end()
+                    << endl;
+
+                visitedPoint.insert(e.end());
+
+                //edgeI = otherEdge(mesh, otherFaceI, edgeI, e.end());
+
+                walkPointConnectedFaces
+                (
+                    mesh,
+                    regionEdges,
+                    e.end(),
+                    otherFaceI,
+                    edgeI,
+                    visitedEdges
+                );
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::regionSide::regionSide
+(
+    const primitiveMesh& mesh,
+    const labelHashSet& region,         // faces of region
+    const labelHashSet& fenceEdges,     // outside edges
+    const label startCellI,
+    const label startFaceI
+)
+:
+    sideOwner_(region.size()),
+    insidePointFaces_(region.size())
+{
+    // Storage for visited faces
+    labelHashSet visitedFace(region.size());
+
+    // Visit all faces on this side of region.
+    // Mark for each face whether owner (or neighbour) has been visited.
+    // Sets sideOwner_
+    visitConnectedFaces
+    (
+        mesh,
+        region,
+        fenceEdges,
+        startCellI,
+        startFaceI,
+        visitedFace
+    );
+
+    //
+    // Visit all internal points on region and mark faces visitable from these.
+    // Sets insidePointFaces_.
+    //
+
+    labelHashSet fencePoints(fenceEdges.size());
+
+    for
+    (
+        labelHashSet::const_iterator iter = fenceEdges.begin();
+        iter != fenceEdges.end();
+        ++iter
+    )
+    {
+        const edge& e = mesh.edges()[iter.key()];
+
+        fencePoints.insert(e.start());
+        fencePoints.insert(e.end());
+    }
+
+    walkAllPointConnectedFaces(mesh, region, fencePoints);
+}
+
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/regionSide/regionSide.H

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::regionSide
+
+Description
+    Determines the 'side' for every face and connected to a 
+    singly-connected (through edges) region of faces. Gets set of faces and
+    a list of mesh edges ('fenceEdges') which should not be crossed.
+    Used in splitting a mesh region.
+
+    Determines:
+    - For every face on the surface: whether the owner was visited
+      from starting face.
+    - List of faces using an internal point of the region visitable by
+      edge-face-edge walking from the correct side of the region.
+
+SourceFiles
+    regionSide.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef regionSide_H
+#define regionSide_H
+
+#include "HashSet.H"
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class primitiveMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class regionSide Declaration
+\*---------------------------------------------------------------------------*/
+
+class regionSide
+{
+    // Private data
+
+        //- For every face on region tells whether the owner is on the
+        //  'regionside'.
+        labelHashSet sideOwner_;
+
+        //- Contains the faces using an internal point and visited face
+        labelHashSet insidePointFaces_;
+
+
+    // Private Member Functions
+
+        //- Step across point to other edge on face
+        static label otherEdge
+        (
+            const primitiveMesh& mesh,
+            const label faceI,
+            const label edgeI,
+            const label pointI
+        );
+
+        //- From faceI, side cellI, cross to other faces/cells by
+        //  face-cell walking and store visited faces and update sideOwner_.
+        void visitConnectedFaces
+        (
+            const primitiveMesh& mesh,
+            const labelHashSet& region,
+            const labelHashSet& fenceEdges,
+            const label cellI,
+            const label faceI,
+            labelHashSet& visitedFace
+        );
+
+        //- From edge on face connected to point on region (regionPointI) cross
+        //  to all other edges using this point by walking across faces
+        //  Does not cross regionEdges so stays on one side of region
+        void walkPointConnectedFaces
+        (
+            const primitiveMesh& mesh,
+            const labelHashSet& regionEdges,
+            const label regionPointI,
+            const label startFaceI,
+            const label startEdgeI,
+            labelHashSet& visitedEdges
+        );
+
+        //- Visits all internal points on region and marks edges reachable
+        //  from sideOwner side (using walkPointConnectedFaces)
+        void walkAllPointConnectedFaces
+        (
+            const primitiveMesh& mesh,
+            const labelHashSet& regionFaces,
+            const labelHashSet& fenceEdges
+        );
+
+public:
+
+    //- Runtime type information
+    ClassName("regionSide");
+
+    // Static Functions
+
+        //- Step across edge onto other face on cell
+        static label otherFace
+        (
+            const primitiveMesh& mesh,
+            const label cellI,
+            const label excludeFaceI,
+            const label edgeI
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        regionSide
+        (
+            const primitiveMesh& mesh,
+            const labelHashSet& region,
+            const labelHashSet& fenceEdges,   // labels of fence edges
+            const label startCell,
+            const label startFace
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            const labelHashSet& sideOwner() const
+            {
+                return sideOwner_;
+            }
+
+            const labelHashSet& insidePointFaces() const
+            {
+                return insidePointFaces_;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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;
 
-
-    //fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
-    //pz = 2: piston points, cutPointZone
-    //cz = 1: moving mask
+    // Zones to add
+    // fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
+    // pz = 2: piston points, cutPointZone
+    // cz = 1: moving mask
 
     label nPorts = scavInCylPatches_.size();
 
 
-    DynamicList<pointZone*> pz(2 + nPorts);
-    DynamicList<faceZone*> fz(3*nPorts + 1);
+    List<pointZone*> pz(nPorts + 2);
+    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,6 +119,19 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
                     boundaryMesh().findPatchID(scavInCylPatches_[patchi])
                 ];
 
+            // Outer slider
+            const polyPatch& outerScav =
+                boundaryMesh()
+                [
+                    boundaryMesh().findPatchID(scavInPortPatches_[patchi])
+                ];
+
+            Pout<< "A: " << innerScav.size() << " " << outerScav.size()
+                << endl;
+            // Add zone if both patches has got faces
+            if (!innerScav.empty() && !outerScav.empty())
+            {
+                // Inner
                 labelList isf(innerScav.size());
 
                 forAll (isf, i)
@@ -142,7 +141,8 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
                 fz[nFaceZones] = new faceZone
                 (
-                scavInCylPatches_[patchi] + "Zone" + Foam::name(patchi + 1),
+                    scavInCylPatches_[patchi] + "Zone"
+                  + Foam::name(patchi + 1),
                     isf,
                     boolList(innerScav.size(), false),
                     nFaceZones,
@@ -151,14 +151,7 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
                 nFaceZones++;
 
-            // Outer slider
-
-            const polyPatch& outerScav =
-                boundaryMesh()
-                [
-                    boundaryMesh().findPatchID(scavInPortPatches_[patchi])
-                ];
-
+                // Outer
                 labelList osf(outerScav.size());
 
                 forAll (osf, i)
@@ -168,7 +161,8 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
                 fz[nFaceZones] = new faceZone
                 (
-                scavInPortPatches_[patchi] + "Zone" + Foam::name(patchi + 1),
+                    scavInPortPatches_[patchi] + "Zone"
+                  + Foam::name(patchi + 1),
                     osf,
                     boolList(outerScav.size(), false),
                     nFaceZones,
@@ -177,6 +171,7 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
                 nFaceZones++;
 
+                // Cut faces
                 fz[nFaceZones] = new faceZone
                 (
                     "cutFaceZone" + Foam::name(patchi + 1),
@@ -188,8 +183,7 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
                 nFaceZones++;
 
-            Info << "cut p" << endl;
-
+                // Cut points
                 pz[nPointZones] = new pointZone
                 (
                     "cutPointZone" + Foam::name(patchi + 1),
@@ -201,8 +195,9 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
                 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;
-
-    Info<< "Adding " << nPointZones << " point, "
+    Pout<< "Adding " << nPointZones << " point, "
         << nFaceZones << " face zones and " << nCellZones
         << " cell zones" << endl;
 
@@ -229,17 +222,39 @@ 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)
         {
+            // Check if patches are present on local processor
+            const label sipID =
+                boundaryMesh().findPatchID(scavInPortPatches_[i]);
+
+            const label sicID =
+                boundaryMesh().findPatchID(scavInCylPatches_[i]);
+
+            if (sipID > -1 && sicID > -1)
+            {
+                if
+                (
+                    boundaryMesh()[sipID].size() > 0
+                 && boundaryMesh()[sicID].size() > 0
+                )
+                {
+                    Pout<< "Adding slider for pair " << scavInPortPatches_[i]
+                        << " and " << scavInCylPatches_[i]
+                        << " with sizes "
+                        << boundaryMesh()[sipID].size() << " "
+                        << boundaryMesh()[sicID].size() << endl;
+
+                    // Patches present.  Add modifier
                     topoChanger_.setSize(topoChanger_.size() + 1);
 
                     topoChanger_.set
@@ -265,6 +280,8 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
                     nMods++;
                 }
             }
+        }
+    }
 
     Info << "Adding " << nMods << " topology modifiers" << endl;
 

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()
-            [boundaryMesh().findPatchID(scavInPortPatches_[patchi])] ==
-            pistonVel;
+        const label curPatchID =
+            boundaryMesh().findPatchID(scavInPortPatches_[patchi]);
+
+        U.boundaryField()[curPatchID] == pistonVel;
     }
 }
 

src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineCalculate.C

@@ -61,8 +61,12 @@ void Foam::twoStrokeEngine::calcMovingMasks() const
     const cellList& c = cells();
     const faceList& f = allFaces();
 
-    const labelList& cellAddr =
-        cellZones()[cellZones().findZoneID("movingCells")];
+    const label movingCellsID = cellZones().findZoneID("movingCells");
+
+    // If moving cell zone is found, mark the vertices
+    if (movingCellsID > -1)
+    {
+        const labelList& cellAddr = cellZones()[movingCellsID];
 
         forAll (cellAddr, cellI)
         {
@@ -79,6 +83,7 @@ void Foam::twoStrokeEngine::calcMovingMasks() const
                 }
             }
         }
+    }
 }
 
 

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);
-
-    Info << "piston face index = " << pistonFaceIndex << endl;
-
-    if(!foundPistonFace)
+    if(pistonFaceIndex == -1)
     {
         FatalErrorIn("Foam::twoStrokeEngine::setVirtualPistonPosition()")
-            << " : cannot find the pistonLayerFaces"
-            << exit(FatalError);
+            << "Cannot find the pistonLayerFaces"
+            << abort(FatalError);
     }
 
-    const labelList& pistonFaces = faceZones()[pistonFaceIndex];
-    forAll(pistonFaces, i)
-    {
-        const face& f = faces()[pistonFaces[i]];
+    const labelList& pistonPoints =
+        faceZones()[pistonFaceIndex]().meshPoints();
 
-        // should loop over facepoints...
-        forAll(f, j)
+    const pointField& p = points();
+
+    forAll (pistonPoints, i)
     {
         virtualPistonPosition() =
-                Foam::max(virtualPistonPosition(), points()[f[j]].z());
-        }
+            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
 
-
-// /* Tommaso, 23/5/2008
-
-    // Find piston mesh modifier
+    // Find piston mesh modifier if present on processor
     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
+        if (pistonLayerID > -1)
+        {
             Info << "Mesh deformation mode" << endl;
             topoChanger_[pistonLayerID].disable();
         }
+    }
     else
     {
         // Activate piston layer
+        if (pistonLayerID > -1)
+        {
             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 (topoChangeMap->hasMotionPoints())
+    if (globalMorphing2)
     {
         Info<< "Topology change; executing pre-motion after "
             << "dynamic layering" << endl;
-            movePoints(topoChangeMap->preMotionPoints());
-            newPoints = topoChangeMap->preMotionPoints();
+
+        if (localMorphing2)
+        {
+            movePoints(topoChangeMap2->preMotionPoints());
+            newPoints = topoChangeMap2->preMotionPoints();
+        }
+        else
+        {
+            movePoints(newPoints);
         }
 
         setV0();
@@ -225,33 +242,18 @@ bool Foam::twoStrokeEngine::update()
     labelList pistonPoints;
     labelList headPoints;
     {
-        label pistonCellIndex = cellZones().findZoneID("pistonCells");
-
-        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];
-
+        // Get cell-point addressing
         const labelListList& cp = cellPoints();
 
         boolList count(newPoints.size(), false);
 
+        // Piston points
+        label pistonCellID = cellZones().findZoneID("pistonCells");
+
+        if (pistonCellID > -1)
+        {
+            const labelList& pistonCells = cellZones()[pistonCellID];
+
             forAll (pistonCells, cellI)
             {
                 const labelList& curCellPoints = cp[pistonCells[cellI]];
@@ -284,13 +286,20 @@ bool Foam::twoStrokeEngine::update()
                     nCounted++;
                 }
             }
+        }
 
         // Repeat for head points
         count = false;
 
+        const label headCellID = cellZones().findZoneID("headCells");
+
+        if (headCellID > -1)
+        {
+            const labelList& headCells = cellZones()[headCellID];
+
             forAll (headCells, cellI)
             {
-            const labelList& curCellPoints = cp[pistonCells[cellI]];
+                const labelList& curCellPoints = cp[headCells[cellI]];
 
                 forAll (curCellPoints, i)
                 {
@@ -299,7 +308,7 @@ bool Foam::twoStrokeEngine::update()
             }
 
             // Count the points
-        nCounted = 0;
+            label nCounted = 0;
             forAll (count, pointI)
             {
                 if (count[pointI] == true)
@@ -321,16 +330,11 @@ 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,42 +425,49 @@ bool Foam::twoStrokeEngine::update()
         makeSlidersLive();
 
         // Changing topology by hand
-        autoPtr<mapPolyMesh> topoChangeMap4 = topoChanger_.changeMesh();
+        autoPtr<mapPolyMesh> topoChangeMap3 = topoChanger_.changeMesh();
 
-        if (topoChangeMap4->morphing())
-        {
-            if (topoChangeMap4->hasMotionPoints())
+        bool localMorphing3 = topoChangeMap3->morphing();
+        bool globalMorphing3 = localMorphing3;
+
+        reduce(globalMorphing3, orOp<bool>());
+
+        if (globalMorphing3)
         {
             Info<< "Topology change; executing pre-motion after "
                 << "sliding attach" << endl;
 
-//                 Info<< "topoChangeMap4->preMotionPoints().size() = "
-//                     << topoChangeMap4->preMotionPoints().size() << nl
-//                     << "allPoints.size() = " << allPoints().size() << nl
-//                     << "points.size() = " << points().size() << endl;
+            // Grab points
+            newPoints = allPoints();
 
-                movePoints(topoChangeMap4->preMotionPoints());
-                newPoints = points();
+            if (localMorphing3)
+            {
+                // If there is layering, pick up correct points
+                if (topoChangeMap3->hasMotionPoints())
+                {
+                    newPoints = topoChangeMap3->preMotionPoints();
                 }
 
-            {
+                // Prepare old points for the move
                 pointField mappedOldPointsNew(allPoints().size());
 
                 mappedOldPointsNew.map
                 (
-                    oldPointsNew, topoChangeMap4->pointMap()
+                    oldPointsNew, topoChangeMap3->pointMap()
                 );
 
                 forAll(scavInPortPatches_, patchi)
                 {
-                    const labelList& cutPointsAddressing =
-                        pointZones()
-                        [
-                            pointZones().findZoneID
+                    // Find cut point zone ID
+                    const label cutPointZoneID = pointZones().findZoneID
                     (
                         "cutPointZone" + Foam::name(patchi + 1)
-                            )
-                        ];
+                    );
+
+                    if (cutPointZoneID > -1)
+                    {
+                        const labelList& cutPointsAddressing =
+                            pointZones()[cutPointZoneID];
 
                         forAll(cutPointsAddressing, i)
                         {
@@ -475,22 +483,40 @@ bool Foam::twoStrokeEngine::update()
                               > virtualPistonPosition()
                             )
                             {
-                            mappedOldPointsNew[cutPointsAddressing[i]].z() =
+                                mappedOldPointsNew
+                                    [cutPointsAddressing[i]].z() =
                                     newPoints[cutPointsAddressing[i]].z();
                             }
                             else
                             {
-                            mappedOldPointsNew[cutPointsAddressing[i]].z() =
-                                newPoints[cutPointsAddressing[i]].z() - deltaZ;
+                                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,10 +1,10 @@
     // Add piston layer addition/removal
     if (piston().patchID().active())
     {
+        Info<< "Adding a layer addition/removal mesh modifier to the piston"
+            << endl;
 
-        Info << "Adding a layer addition/removal mesh modifier to the piston" << endl;
-    
-        topoChanger_.setSize(nMods+1);
+        topoChanger_.setSize(nMods + 1);
 
         topoChanger_.set
         (
@@ -20,6 +20,6 @@
             )
         );
         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
     );
@@ -458,8 +458,7 @@ Foam::simpleEngineTopoFvMesh::simpleEngineTopoFvMesh
     const IOobject& io
 )
 :
-    topoChangerFvMesh(io),
-    engineTime_(refCast<const engineTime>(time())),
+    engineTopoChangerMesh(io),
     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 "engineTime.H"
+#include "engineTopoChangerMesh.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)
+        {}
 };