some changes

growth_eqns.py

@@ -24,19 +24,27 @@ class Growth(object):
         self.fminus1 = df.Function(self.function_space)           # f at t = -1
         self.f0 = df.Function(self.function_space)                 # f at t =0
 
+        df.plot(self.mesh)
+        plt.show()
+
     def advection(self,c,v,tc):
         return (df.inner((v*c).dx(0),tc))
 
     def reaction_c(self,c,tc):
         return df.inner(self.turnover_c*(c-self.mean_concentration),tc)
 
-    def reaction_rho(self,rho):
-        return self.turnover_rho*rho
+    def reaction_rho(self,rho, trho):
+        return df.inner(self.turnover_rho*(rho - self.mean_density), trho)
 
     def reaction_diffusion_c(self,c,tc):
-        return ( self.diffusion*df.inner(c.dx(0),tc.dx(0))
+        return ( self.diffusion_c*df.inner(c.dx(0),tc.dx(0))
                 + self.reaction_c(c,tc)
                 )
+    def reaction_diffusion_rho(self,c,tc):
+        return ( self.diffusion_rho*df.inner(c.dx(0),tc.dx(0))
+            + self.reaction_c(c,tc)
+            )
+
 
     def setup_initial_conditions(self, u0,rho0,c0):
         
@@ -58,20 +66,21 @@ class Growth(object):
         self.fminus1.assign(self.f0)
 
     def setup_weak_forms(self):
-        uminus1, _, rhominus1, cminus1  = df.split(self.fminus1)
+        uminus1, vminus1, rhominus1, cminus1  = df.split(self.fminus1)
         u0, _ , rho0, c0 = df.split(self.f0)
         u,  v, rho, c = df.split(self.f)
         tu, tv, trho, tc = df.TestFunctions(self.function_space)
 
         uform = (df.inner((u-u0)/self.timestep,tu) 
-                - df.inner(v,tu)
+                -(3/2)*df.inner(v,tu)
+                + (1/2)*df.inner(vminus1,tu)
         )
-
         vform = (df.inner(v,tv)
-                + self.youngs_modulus*df.inner(u.dx(0),tv.dx(0))
-                - self.b*df.inner(rho.dx(0),tv)
+                + (3/2)*self.youngs_modulus*df.inner(u.dx(0),tv.dx(0))
+                - (1/2)*self.youngs_modulus*df.inner(uminus1.dx(0),tv.dx(0))
+                - (3/2)*self.b*df.inner(rho.dx(0),tv)
+                + (1/2)*self.b*df.inner(rho.dx(0),tv)
         )
-
         cform = (df.inner((c-c0)/self.timestep, tc)
                 - 3/2 * self.advection(c0, v, tc)
                 + 1/2 * self.advection(cminus1, v, tc)
@@ -79,13 +88,12 @@ class Growth(object):
                 + 3/8 * self.reaction_diffusion_c(c0, tc)
                 + 1/16 * self.reaction_diffusion_c(cminus1, tc)
             )
-
         rhoform = (df.inner((rho-rho0)/self.timestep,trho)
                 - (3/2)*self.advection(rho0,v,trho)
                 + (1/2)*self.advection(rhominus1,v, trho)
-                + (9/16)*df.inner(self.reaction_rho(rho),trho)
-                + (3/8)*df.inner(self.reaction_rho(rho0),trho)
-                + (1/16)*df.inner(self.reaction_rho(rhominus1),trho)
+                + (9/16)*self.reaction_diffusion_rho(rho, trho)
+                + (3/8)*self.reaction_diffusion_rho(rho0, trho)
+                + (1/16)*self.reaction_diffusion_rho(rhominus1, trho)
         )
         self.form = (uform + vform + rhoform + cform) * df.dx
 
@@ -106,21 +114,25 @@ class Growth(object):
         rho_array[0] = rho.compute_vertex_values(self.mesh)
         c_array[0] = c.compute_vertex_values(self.mesh)      
 
-        for i in progressbar.progressbar(range(0, len(times))):
+        for i in progressbar.progressbar(range(1, len(times))):
             df.solve(self.form == 0,self.f)
             u, v, rho, c = self.f.split(deepcopy=True)
             u_array[i] = u.compute_vertex_values(self.mesh)
             v_array[i] = v.compute_vertex_values(self.mesh)
             rho_array[i] = rho.compute_vertex_values(self.mesh)
             c_array[i] = c.compute_vertex_values(self.mesh)
-            self.fminus1.assign(self.f0)
             self.f0.assign(self.f)
-
+            self.fminus1.assign(self.f0)
+            df.ALE.move(self.mesh,u)
+        df.plot(self.mesh)
+        plt.show()
         return (u_array, v_array, rho_array, c_array)
 
 
 
 if __name__ == '__main__':
+    import dolfin as df
+
     import json
     import matplotlib.pyplot as plt
     from matplotlib.widgets import Slider
@@ -130,6 +142,7 @@ if __name__ == '__main__':
     assert(params['dimension']==1)
 
     g = Growth(params)
+
     u0 = df.Expression('0.5*cos(x[0]/2)', degree=1)  
     rho0 = df.Expression('1+0.2*cos(x[0])', degree=1)
     c0 = df.Expression('cos(2*x[0])+cos(x[0])', degree=1)
@@ -153,7 +166,7 @@ if __name__ == '__main__':
     rhoplot, = axrho.plot(x, rho[0])
     cplot, = axc.plot(x, c[0])
 
-    
+
     def update(value):
         ti = np.abs(times-value).argmin()
         uplot.set_ydata(u[ti])