From a00cf16bd7e709fc6e2bc2aac3d50e8dd3c76745 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Thu, 8 Dec 2022 17:37:37 +0000
Subject: [PATCH] solvers::isothermalFluid: Updated for mesh-motion with NCC
 and mesh-to-mesh mapping

---
 .../modules/isothermalFluid/isothermalFluid.C        | 10 +++-------
 .../modules/isothermalFluid/isothermalFluid.H        |  8 ++++++++
 .../solvers/modules/isothermalFluid/moveMesh.C       |  5 +++++
 .../GeometricFields/GeometricField/GeometricField.C  | 12 ++++++------
 src/finiteVolume/fvMesh/fvMesh.C                     |  4 ++--
 5 files changed, 24 insertions(+), 15 deletions(-)

diff --git a/applications/solvers/modules/isothermalFluid/isothermalFluid.C b/applications/solvers/modules/isothermalFluid/isothermalFluid.C
index 90fd6f37b1..c797b2d325 100644
--- a/applications/solvers/modules/isothermalFluid/isothermalFluid.C
+++ b/applications/solvers/modules/isothermalFluid/isothermalFluid.C
@@ -302,16 +302,13 @@ void Foam::solvers::isothermalFluid::preSolve()
 
     fvModels().preUpdateMesh();
 
-    //- Pointer to the vol momentum field
-    //  used for mesh-change to set rhoUf for introduced faces
-    autoPtr<volVectorField> rhoU;
-    autoPtr<volVectorField> rhoU0;
-
     // Store momentum to set rhoUf for introduced faces
     if (mesh.topoChanging())
     {
         rhoU = new volVectorField("rhoU", rho*U);
 
+        InfoInFunction << "rhoUf().nOldTimes() " << rhoUf().nOldTimes() << endl;
+
         if (rhoUf().nOldTimes() > 1)
         {
             rhoU0 = new volVectorField("rhoU_0", rho.oldTime()*U.oldTime());
@@ -321,8 +318,7 @@ void Foam::solvers::isothermalFluid::preSolve()
     // Update the mesh for topology change, mesh to mesh mapping
     mesh.update();
 
-    rhoU.clear();
-    rhoU0.clear();
+    InfoInFunction << "rhoUf().nOldTimes() " << rhoUf().nOldTimes() << endl;
 }
 
 
diff --git a/applications/solvers/modules/isothermalFluid/isothermalFluid.H b/applications/solvers/modules/isothermalFluid/isothermalFluid.H
index 3cc89f95e7..47c165f523 100644
--- a/applications/solvers/modules/isothermalFluid/isothermalFluid.H
+++ b/applications/solvers/modules/isothermalFluid/isothermalFluid.H
@@ -151,6 +151,14 @@ protected:
         //  Used for buoyant simulations only
         tmp<volVectorField> netForce;
 
+        //- Pointer to the vol momentum field
+        //  used for mesh-change to set rhoUf for introduced faces
+        autoPtr<volVectorField> rhoU;
+
+        //- Pointer to the old-time vol momentum field
+        //  used for mesh-change to set rhoUf.oldTime() for introduced faces
+        autoPtr<volVectorField> rhoU0;
+
         //- Pointer to the surface momentum field
         //  used to recreate the flux after mesh-change
         autoPtr<surfaceVectorField> rhoUf;
diff --git a/applications/solvers/modules/isothermalFluid/moveMesh.C b/applications/solvers/modules/isothermalFluid/moveMesh.C
index 5e422dc5c5..f85f5ca97a 100644
--- a/applications/solvers/modules/isothermalFluid/moveMesh.C
+++ b/applications/solvers/modules/isothermalFluid/moveMesh.C
@@ -35,6 +35,11 @@ bool Foam::solvers::isothermalFluid::moveMesh()
         // Move the mesh
         mesh.move();
 
+        // The rhoU and rhoU0 fields can be cleared following mesh-motion
+        // now the mesh has been re-stitched as necessary
+        rhoU.clear();
+        rhoU0.clear();
+
         if (mesh.changing())
         {
             if (buoyancy.valid())
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
index 701b7aa465..0e73d60205 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
@@ -418,7 +418,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
             << "Constructing as copy" << endl << this->info() << endl;
     }
 
-    if (gf.field0Ptr_)
+    if (gf.field0Ptr_ && notNull(gf.field0Ptr_))
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (
@@ -448,7 +448,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
             << "Constructing by moving" << endl << this->info() << endl;
     }
 
-    if (gf.field0Ptr_)
+    if (gf.field0Ptr_ && notNull(gf.field0Ptr_))
     {
         field0Ptr_ = gf.field0Ptr_;
         gf.field0Ptr_ = nullptr;
@@ -506,7 +506,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
             << endl << this->info() << endl;
     }
 
-    if (!readIfPresent() && gf.field0Ptr_)
+    if (!readIfPresent() && gf.field0Ptr_ && notNull(gf.field0Ptr_))
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (
@@ -568,7 +568,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
             << endl << this->info() << endl;
     }
 
-    if (!readIfPresent() && gf.field0Ptr_)
+    if (!readIfPresent() && gf.field0Ptr_ && notNull(gf.field0Ptr_))
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (
@@ -631,7 +631,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     boundaryField_ == gf.boundaryField_;
 
-    if (!readIfPresent() && gf.field0Ptr_)
+    if (!readIfPresent() && gf.field0Ptr_ && notNull(gf.field0Ptr_))
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (
@@ -673,7 +673,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 
     boundaryField_ == gf.boundaryField_;
 
-    if (!readIfPresent() && gf.field0Ptr_)
+    if (!readIfPresent() && gf.field0Ptr_ && notNull(gf.field0Ptr_))
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (
diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C
index 9f556ea954..61f996c494 100644
--- a/src/finiteVolume/fvMesh/fvMesh.C
+++ b/src/finiteVolume/fvMesh/fvMesh.C
@@ -645,8 +645,6 @@ bool Foam::fvMesh::dynamic() const
 
 bool Foam::fvMesh::update()
 {
-    if (!conformal()) stitcher_->disconnect(true, true);
-
     if
     (
         stitcher_->stitches()
@@ -657,6 +655,8 @@ bool Foam::fvMesh::update()
         nullOldestTimeFields();
     }
 
+    if (!conformal()) stitcher_->disconnect(true, true);
+
     const bool hasV00 = V00Ptr_;
     deleteDemandDrivenData(V00Ptr_);