Engine updates: parallelisation of topo changes

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/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

@@ -0,0 +1,60 @@
+if (piston().patchID().active())
+{
+    // 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());
+
+    boolList flipPistonFaces(pistonFaceSet.toc().size(), false);
+    label nSet = pistonFaceSet.toc().size();
+    label nFlip = 0;
+
+    forAll (flipPistonFaces, facei)
+    {
+        scalar scalProd =
+            (faceAreas()[pistonFaceSet.toc()[facei]] & vector(0, 0, 1));
+
+        if (scalProd < 0)
+        {
+            flipPistonFaces[facei] = true;
+            nFlip++;
+        }
+    }
+
+    Info << "nSet = " << nSet << endl;
+    Info << "nFlip = " << nFlip << endl;
+
+    fz[nFaceZones] =
+    (
+        new faceZone
+        (
+            "pistonLayerFaces",
+            pistonFaceSet.toc(),
+            flipPistonFaces,
+            nFaceZones,
+            faceZones()
+        )
+    );
+
+    nFaceZones++;
+
+    cellSet movingPistonCells(*this, piston().pistonCellSetName());
+
+    Info<< "Adding piston cell set" << endl;
+    cz[nCellZones] =
+    (
+        new cellZone
+        (
+            "pistonCells",
+            movingPistonCells.toc(),
+            nCellZones,
+            cellZones()
+        )
+    );
+
+    nCellZones++;
+}

src/engine/engineTopoChangerMesh/twoStrokeEngine/addTwoStrokeEngineModifiers.C

@@ -62,25 +62,24 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
         return;
     }
 
-    Info << "checkAndCalculate()" << endl;
     checkAndCalculate();
 
     Info<< "Time = " << engTime().theta() << endl
         << "Adding zones to the engine mesh" << endl;
 
-
+    // Zones to add
-    //fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
+    // fz = 4: virtual piston, outSidePort, insidePort, cutFaceZone
-    //pz = 2: piston points, cutPointZone
+    // pz = 2: piston points, cutPointZone
-    //cz = 1: moving mask
+    // cz = 1: moving mask
 
     label nPorts = scavInCylPatches_.size();
 
 
-    DynamicList<pointZone*> pz(2 + nPorts);
+    List<pointZone*> pz(nPorts + 2);
-    DynamicList<faceZone*> fz(3*nPorts + 1);
+    List<faceZone*> fz(3*nPorts + 1);
 
     // Added piston cells and head cells
-    DynamicList<cellZone*> cz(3);
+    List<cellZone*> cz(3);
 
     label nPointZones = 0;
     label nFaceZones = 0;
@@ -88,24 +87,11 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
 
     Info << "Adding piston layer faces" << endl;
 
-#   include "addPistonLayerHrv.H"
+#   include "addPistonLayer.H"
 
-//  adding head points (does not move)
+    // Add head points that do not move
 
     {
-//         pointSet headPointSet(*this, headPointsSetName_);
-
-//         pz[nPointZones] =
-//             new pointZone
-//             (
-//                 "headPoints",
-//                 headPointSet.toc(),
-//                 nPointZones,
-//                 pointZones()
-//             );
-
-//         nPointZones++;
-
         cellSet headCellSet(*this, headCellsSetName_);
 
         cz[nCellZones] =
@@ -120,9 +106,9 @@ void Foam::twoStrokeEngine::addZonesAndModifiers()
         nCellZones++;
     }
 
-//  Sliding interface for scavenging ports
+    //  Sliding interface for scavenging ports
 
-    if(nPorts > 0)
+    if (nPorts > 0)
     {
         forAll(scavInCylPatches_, patchi)
         {
@@ -133,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
+                // Outer
-
-            const polyPatch& outerScav =
-                boundaryMesh()
-                [
-                    boundaryMesh().findPatchID(scavInPortPatches_[patchi])
-                ];
-
                 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;
+    Pout<< "Adding " << nPointZones << " point, "
-
-    Info<< "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()
+        const label curPatchID =
-            [boundaryMesh().findPatchID(scavInPortPatches_[patchi])] ==
+            boundaryMesh().findPatchID(scavInPortPatches_[patchi]);
-            pistonVel;
+
+        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 =
+    const label movingCellsID = cellZones().findZoneID("movingCells");
-        cellZones()[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);
+    if(pistonFaceIndex == -1)
-
-    Info << "piston face index = " << pistonFaceIndex << endl;
-
-    if(!foundPistonFace)
     {
         FatalErrorIn("Foam::twoStrokeEngine::setVirtualPistonPosition()")
-            << " : cannot find the pistonLayerFaces"
+            << "Cannot find the pistonLayerFaces"
-            << exit(FatalError);
+            << abort(FatalError);
     }
 
-    const labelList& pistonFaces = faceZones()[pistonFaceIndex];
+    const labelList& pistonPoints =
-    forAll(pistonFaces, i)
+        faceZones()[pistonFaceIndex]().meshPoints();
-    {
-        const face& f = faces()[pistonFaces[i]];
 
-        // should loop over facepoints...
+    const pointField& p = points();
-        forAll(f, j)
+
+    forAll (pistonPoints, i)
     {
         virtualPistonPosition() =
-                Foam::max(virtualPistonPosition(), points()[f[j]].z());
+            Foam::max(virtualPistonPosition_, p[pistonPoints[i]].z());
-        }
     }
 
     reduce(virtualPistonPosition(), maxOp<scalar>());
 }
+
+
+// ************************************************************************* //

src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineMove.C

@@ -41,11 +41,11 @@ void Foam::twoStrokeEngine::makeLayersLive()
     // Enable layering
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(layerAdditionRemoval))
+        if (isA<layerAdditionRemoval>(morphs[modI]))
         {
             morphs[modI].enable();
         }
-        else if (typeid(morphs[modI]) == typeid(slidingInterface))
+        else if (isA<slidingInterface>(morphs[modI]))
         {
             morphs[modI].disable();
         }
@@ -59,6 +59,7 @@ void Foam::twoStrokeEngine::makeLayersLive()
     }
 }
 
+
 void Foam::twoStrokeEngine::makeSlidersLive()
 {
     const polyTopoChanger& morphs = topoChanger_;
@@ -66,11 +67,11 @@ void Foam::twoStrokeEngine::makeSlidersLive()
     // Enable sliding interface
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(layerAdditionRemoval))
+        if (isA<layerAdditionRemoval>(morphs[modI]))
         {
             morphs[modI].disable();
         }
-        else if (typeid(morphs[modI]) == typeid(slidingInterface))
+        else if (isA<slidingInterface>(morphs[modI]))
         {
             morphs[modI].enable();
         }
@@ -93,10 +94,9 @@ bool Foam::twoStrokeEngine::attached() const
 
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(slidingInterface))
+        if (isA<slidingInterface>(morphs[modI]))
         {
-            result =
+            result = result
-                result
              || refCast<const slidingInterface>(morphs[modI]).attached();
         }
     }
@@ -104,7 +104,7 @@ bool Foam::twoStrokeEngine::attached() const
     // Check thal all sliders are in sync (debug only)
     forAll (morphs, modI)
     {
-        if (typeid(morphs[modI]) == typeid(slidingInterface))
+        if (isA<slidingInterface>(morphs[modI]))
         {
             if
             (
@@ -120,6 +120,9 @@ bool Foam::twoStrokeEngine::attached() const
         }
     }
 
+    // Sync across processors
+    reduce(result, orOp<bool>());
+
     return result;
 }
 
@@ -133,13 +136,25 @@ bool Foam::twoStrokeEngine::update()
         makeSlidersLive();
 
         // Changing topology by hand
-        autoPtr<mapPolyMesh> topoChangeMap5 = topoChanger_.changeMesh();
+        autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
 
-        if (topoChangeMap5->hasMotionPoints() && topoChangeMap5->morphing())
+        bool localMorphing1 = topoChangeMap1->morphing();
+
+        // Note: Since we are detaching, global morphing is always true
+        // HJ, 7/Mar/2011
+
+        if (localMorphing1)
         {
             Info << "Topology change; executing pre-motion after "
                 << "sliding detach" << endl;
-            movePoints(topoChangeMap5->preMotionPoints());
+            movePoints(topoChangeMap1->preMotionPoints());
+        }
+        else
+        {
+            pointField newPoints = allPoints();
+
+            // Dummy motion
+            movePoints(newPoints);
         }
 
         Info << "sliding interfaces successfully decoupled!!!" << endl;
@@ -154,25 +169,14 @@ bool Foam::twoStrokeEngine::update()
     // Piston Layering
 
     makeLayersLive();
-    // Changing topology by hand
 
-
+    // Find piston mesh modifier if present on processor
-// /* Tommaso, 23/5/2008
-
-    // Find piston mesh modifier
     const label pistonLayerID =
         topoChanger_.findModifierID("pistonLayer");
 
-    if (pistonLayerID < 0)
-    {
-        FatalErrorIn("void engineFvMesh::moveAndMorph()")
-            << "Piston modifier not found."
-            << abort(FatalError);
-    }
-
     scalar minLayerThickness = piston().minLayer();
     scalar deltaZ = engTime().pistonDisplacement().value();
-    virtualPistonPosition() += deltaZ;
+    virtualPistonPosition_ += deltaZ;
 
     Info << "virtualPistonPosition = " << virtualPistonPosition()
     << ", deckHeight = " << deckHeight()
@@ -181,32 +185,45 @@ bool Foam::twoStrokeEngine::update()
     if (realDeformation())
     {
         // Dectivate piston layer
+        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 (globalMorphing2)
-    {
-        if (topoChangeMap->hasMotionPoints())
     {
         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");
+        // Get cell-point addressing
-
-        if (pistonCellIndex < 0)
-        {
-            FatalErrorIn("bool twoStrokeEngine::update()")
-                << "Cannot find cell zone pistonCells"
-                << abort(FatalError);
-        }
-
-
-        const labelList& pistonCells = cellZones()[pistonCellIndex];
-
-        label headCellIndex = cellZones().findZoneID("headCells");
-
-        if (headCellIndex < 0)
-        {
-            FatalErrorIn("bool twoStrokeEngine::update()")
-                << "Cannot find cell zone headCells"
-                << abort(FatalError);
-        }
-
-        const labelList& headCells = cellZones()[headCellIndex];
-
         const labelListList& cp = cellPoints();
 
         boolList count(newPoints.size(), false);
 
+        // 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())
+        bool localMorphing3 = topoChangeMap3->morphing();
-        {
+        bool globalMorphing3 = localMorphing3;
-            if (topoChangeMap4->hasMotionPoints())
+
+        reduce(globalMorphing3, orOp<bool>());
+
+        if (globalMorphing3)
         {
             Info<< "Topology change; executing pre-motion after "
                 << "sliding attach" << endl;
 
-//                 Info<< "topoChangeMap4->preMotionPoints().size() = "
+            // Grab points
-//                     << topoChangeMap4->preMotionPoints().size() << nl
+            newPoints = allPoints();
-//                     << "allPoints.size() = " << allPoints().size() << nl
-//                     << "points.size() = " << points().size() << endl;
 
-                movePoints(topoChangeMap4->preMotionPoints());
+            if (localMorphing3)
-                newPoints = points();
+            {
+                // 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 =
+                    // Find cut point zone ID
-                        pointZones()
+                    const label cutPointZoneID = pointZones().findZoneID
-                        [
-                            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() =
+                                mappedOldPointsNew
-                                newPoints[cutPointsAddressing[i]].z() - deltaZ;
+                                    [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),
+    engineTopoChangerMesh(io),
-    engineTime_(refCast<const engineTime>(time())),
     valves_(*this, engineTime_.engineDict().lookup("valves")),
     piston_(*this, engineTime_.engineDict().subDict("piston")),
     msPtr_(motionSolver::New(*this)),

src/engine/simpleEngineTopoFvMesh/simpleEngineTopoFvMesh.H

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