From 8e7c357baa9ed44d11fb404415bc6f7950dd6abe Mon Sep 17 00:00:00 2001
From: Alvin Zhang <7065326+zhang-alvin@users.noreply.github.com>
Date: Fri, 15 May 2020 12:40:25 -0400
Subject: [PATCH] FIX: recompute geometric info to get correct element/node
 diameters. also included longer history for dt (needed for adaptivity)
 (#1182)

Referring to one of a few problems described in #1179. This PR is simply to get correctness in place. Later PRs will include tests based on a pseudo-adapted mesh.

Update for simplex meshes only in postStep()
---
 proteus/mprans/MoveMesh.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/proteus/mprans/MoveMesh.py b/proteus/mprans/MoveMesh.py
index 01145746dc..e6059626c2 100644
--- a/proteus/mprans/MoveMesh.py
+++ b/proteus/mprans/MoveMesh.py
@@ -10,6 +10,7 @@
 from proteus.Transport import OneLevelTransport, TC_base, NonlinearEquation
 from proteus import Quadrature, FemTools, Comm, Archiver, cfemIntegrals
 from proteus.Profiling import logEvent, memory
+from proteus import cmeshTools
 from . import cArgumentsDict
 
 class Coefficients(proteus.TransportCoefficients.TC_base):
@@ -34,6 +35,7 @@ def __init__(self,
         self.gravityStep = True
         self.meIndex = meIndex
         self.dt_last = None
+        self.dt_last_last = None
         self.solidsList = []
         self.nullSpace = nullSpace
         if initialize:
@@ -137,7 +139,21 @@ def postStep(self, t, firstStep=False):
         else:
             dt = self.dt_last
         self.mesh.nodeVelocityArray /= dt
+        #this is needed for proper restarting
+        if(self.dt_last is not None):
+            self.dt_last_last = self.dt_last
+
         self.dt_last = self.model.timeIntegration.dt
+
+        #update nodal/element diameters:
+        #TODO: unclear if this needs to apply to all mesh types
+        if self.nd == 2:  
+           cmeshTools.computeGeometricInfo_triangle(self.mesh.subdomainMesh.cmesh)
+           self.mesh.buildFromC(self.mesh.cmesh)
+        if self.nd == 3:  
+           cmeshTools.computeGeometricInfo_tetrahedron(self.mesh.subdomainMesh.cmesh)
+           self.mesh.buildFromC(self.mesh.cmesh)
+
         copyInstructions = {'clear_uList': True}
         return copyInstructions