fixed vform -- no imex for v

no_c_fixed_bdry_fenics.py

@@ -13,8 +13,8 @@ class Growth(object):
 
         # create mesh, mixed element and function space
         self.mesh = df.IntervalMesh(self.resolution, 0.0, 2.0*np.pi*self.system_size_by_2PI)
-        conc_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
-        mixed_element = df.MixedElement([conc_element, conc_element, conc_element]) # u, v, rho
+        scalar_element = df.FiniteElement('P', self.mesh.ufl_cell(), 1)
+        mixed_element = df.MixedElement([scalar_element, scalar_element, scalar_element]) # u, v, rho
 
         # define function space with this mixed element
         self.function_space = df.FunctionSpace(self.mesh, mixed_element)
@@ -32,8 +32,14 @@ class Growth(object):
     def advection(self, conc, vel, tconc):
         return (df.inner((vel * conc).dx(0), tconc))
 
-    def active_stress(self, rho, tv):
-        return self.lamda * df.inner((rho/(rho + self.saturation_rho)).dx(0), tv)
+    def active_stress(self, rho):
+        return self.lamda * rho/(rho + self.saturation_rho)
+    
+    def passive_stress(self, u, v):
+        return (self.elasticity * u.dx(0) + self.viscosity * v.dx(0))
+    
+    def stress(self, u, v, rho):
+        return (self.passive_stress(u, v) + self.active_stress(rho))
 
     def reaction_rho(self, rho, trho):
         return self.turnover_rho * df.inner(rho - self.average_rho, trho)
@@ -46,15 +52,14 @@ class Growth(object):
         u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
         rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
 
-        v0_function_space = df.FunctionSpace(self.mesh, 'P', 1)
+        v0_function_space = self.function_space.sub(1).collapse()
 
         v0  = df.Function(v0_function_space)
         tv0 = df.TestFunction(v0_function_space)
-        v0form = (df.inner(v0,tv0)
-                + self.elasticity * df.inner(u0.dx(0),tv0.dx(0))
-                + self.viscosity * df.inner(v0.dx(0), tv0.dx(0))
-                - self.lamda * df.inner((rho0/(rho0 + self.saturation_rho)).dx(0),tv0)
+        v0form = (self.friction * df.inner(v0, tv0) 
+                    + df.inner(self.stress(u0, v0, rho0), tv0.dx(0))
                 ) * df.dx
+                
         bcv0  = df.DirichletBC(v0_function_space, df.Constant(0.0), 'on_boundary')
         df.solve(v0form == 0, v0, bcv0)
         df.assign(self.f0, [u0, v0, rho0])
@@ -71,18 +76,8 @@ class Growth(object):
                 - 3/8 * df.inner(v0, tu)
                 - 1/16 * df.inner(vminus1, tu)
         )
-        vform = ( df.inner(v, tv)
-                + 9/16 * self.elasticity * df.inner(u.dx(0), tv.dx(0))
-                + 3/8 * self.elasticity * df.inner(u0.dx(0), tv.dx(0))
-                + 1/16 * self.elasticity * df.inner(uminus1.dx(0), tv.dx(0))
-
-                - 9/16 * self.active_stress(rho0, tv)
-                - 3/8 * self.active_stress(rhominus1, tv)
-                - 1/16 * self.active_stress(rhominus1, tv)
-
-                + 9/16 * self.viscosity * df.inner(v.dx(0),tv.dx(0))
-                + 3/8 * self.viscosity * df.inner(v0.dx(0),tv.dx(0))
-                + 1/16 * self.viscosity * df.inner(vminus1.dx(0),tv.dx(0))
+        vform = (self.friction * df.inner(v, tv) 
+                + df.inner(self.stress(u, v, rho), tv.dx(0))
                 )
 
         rhoform = (df.inner((rho - rho0)/self.timestep, trho)
@@ -104,7 +99,7 @@ class Growth(object):
         self.setup_initial_conditions(u0, rho0)
         self.setup_weak_forms()
 
-        u,v,rho = self.f0.split(deepcopy=True)
+        u, v, rho = self.f0.split(deepcopy=True)
         u_array[0] = u.compute_vertex_values(self.mesh)
         v_array[0] = v.compute_vertex_values(self.mesh) 
         rho_array[0] = rho.compute_vertex_values(self.mesh)
@@ -120,23 +115,18 @@ class Growth(object):
         return (u_array, v_array, rho_array)
 
 if __name__ == '__main__':
-    import dolfin as df
-
     import json
     import matplotlib.pyplot as plt
     from matplotlib.widgets import Slider
     
     with open('growth_parameters.json') as jsonFile:
         params = json.load(jsonFile)
-    assert(params['dimension']==1)
 
     g = Growth(params)
 
     u0 = df.Expression('sin(x[0])', degree=1)
-    u0 = df.interpolate(u0, g.function_space.sub(0).collapse())
 
-    rho0 = df.Expression('cos(x[0])+cos(x[0]/2)', degree=1) 
-    rho0 = df.interpolate(rho0, g.function_space.sub(2).collapse())
+    rho0 = df.Expression('rho0 + 0.1 * (cos(x[0])+cos(x[0]/2))', rho0 =g.average_rho, degree=1) 
 
     (u, v, rho) = g.solve(u0, rho0)