From cc22953ac4674fce9740915f015d4041bd300f7a Mon Sep 17 00:00:00 2001
From: Hrvoje Jasak <h.jasak@wikki.co.uk>
Date: Wed, 17 Nov 2010 12:50:34 +0000
Subject: [PATCH] Fixed dynamic refinement tutorials: direct topo change does
 not deal with face mapping correctly

---
 .../multiphase/interDyMFoam/correctPhi.H      | 13 ++++++-
 .../multiphase/interDyMFoam/interDyMFoam.C    | 38 +++++++++----------
 2 files changed, 31 insertions(+), 20 deletions(-)

diff --git a/applications/solvers/multiphase/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interDyMFoam/correctPhi.H
index 88472df9b..287afec5c 100644
--- a/applications/solvers/multiphase/interDyMFoam/correctPhi.H
+++ b/applications/solvers/multiphase/interDyMFoam/correctPhi.H
@@ -38,7 +38,14 @@
         pcorrTypes
     );
 
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
+    dimensionedScalar rAUf
+    (
+        "(1|A(U))",
+        dimTime/rho.dimensions(),
+        runTime.deltaT().value()
+    );
+
+    phi = (fvc::interpolate(U) & mesh.Sf());
 
     adjustPhi(phi, U, pcorr);
 
@@ -59,4 +66,8 @@
     }
 
 #   include "continuityErrs.H"
+#   include "CourantNo.H"
+
+    // Recalculate rhoPhi from rho
+    rhoPhi = fvc::interpolate(rho)*phi;
 }
diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
index 285059a8c..1c51f7a8f 100644
--- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
+++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C
@@ -45,17 +45,17 @@ Description
 
 int main(int argc, char *argv[])
 {
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createDynamicFvMesh.H"
-    #include "readGravitationalAcceleration.H"
-    #include "readPISOControls.H"
-    #include "initContinuityErrs.H"
-    #include "createFields.H"
-    #include "readTimeControls.H"
-    #include "correctPhi.H"
-    #include "CourantNo.H"
-    #include "setInitialDeltaT.H"
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+#   include "readGravitationalAcceleration.H"
+#   include "readPISOControls.H"
+#   include "initContinuityErrs.H"
+#   include "createFields.H"
+#   include "readTimeControls.H"
+#   include "correctPhi.H"
+#   include "CourantNo.H"
+#   include "setInitialDeltaT.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -63,13 +63,13 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
-        #include "CourantNo.H"
+#       include "readControls.H"
+#       include "CourantNo.H"
 
         // Make the fluxes absolute
         fvc::makeAbsolute(phi, U);
 
-        #include "setDeltaT.H"
+#       include "setDeltaT.H"
 
         runTime++;
 
@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && correctPhi)
         {
-            #include "correctPhi.H"
+#           include "correctPhi.H"
         }
 
         // Make the fluxes relative to the mesh motion
@@ -100,19 +100,19 @@ int main(int argc, char *argv[])
 
         if (mesh.changing() && checkMeshCourantNo)
         {
-            #include "meshCourantNo.H"
+#           include "meshCourantNo.H"
         }
 
         twoPhaseProperties.correct();
 
-        #include "alphaEqnSubCycle.H"
+#       include "alphaEqnSubCycle.H"
 
-        #include "UEqn.H"
+#       include "UEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-            #include "pEqn.H"
+#           include "pEqn.H"
         }
 
         p = pd + rho*gh;