made active stress a saturating function of rho

no_c_fixed_bdry_fenics.py

@@ -19,7 +19,6 @@ class Growth(object):
         # define function space with this mixed element
         # define function space with this mixed element
         self.function_space = df.FunctionSpace(self.mesh, mixed_element)
         self.function_space = df.FunctionSpace(self.mesh, mixed_element)
 
 
-
         # dirichlet boundaries for u, v
         # dirichlet boundaries for u, v
         bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary')   # u,v at boundary = 0
         bc1 = df.DirichletBC(self.function_space.sub(0), df.Constant(0.0), 'on_boundary')   # u,v at boundary = 0
         bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
         bc2 = df.DirichletBC(self.function_space.sub(1), df.Constant(0.0), 'on_boundary')
@@ -28,11 +27,14 @@ class Growth(object):
         # define the reqd functions on this function space
         # define the reqd functions on this function space
         self.f = df.Function(self.function_space)                 # f at current time
         self.f = df.Function(self.function_space)                 # f at current time
         self.fminus1 = df.Function(self.function_space)           # f at t = -1
         self.fminus1 = df.Function(self.function_space)           # f at t = -1
-        self.f0 = df.Function(self.function_space)                 # f at t =0
+        self.f0 = df.Function(self.function_space)                # f at t =0
 
 
     def advection(self, conc, vel, tconc):
     def advection(self, conc, vel, tconc):
         return (df.inner((vel * conc).dx(0), 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 reaction_rho(self, rho, trho):
     def reaction_rho(self, rho, trho):
         return self.turnover_rho * df.inner(rho - self.average_rho, trho)
         return self.turnover_rho * df.inner(rho - self.average_rho, trho)
 
 
@@ -44,7 +46,6 @@ class Growth(object):
         u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
         u0 = df.interpolate(u0, self.function_space.sub(0).collapse())
         rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
         rho0 = df.interpolate(rho0, self.function_space.sub(2).collapse())
 
 
-        # v0_function_space = df.FunctionSpace(self.mesh, 'P', 1, constrained_domain = self.pbc)
         v0_function_space = df.FunctionSpace(self.mesh, 'P', 1)
         v0_function_space = df.FunctionSpace(self.mesh, 'P', 1)
 
 
         v0  = df.Function(v0_function_space)
         v0  = df.Function(v0_function_space)
@@ -52,7 +53,7 @@ class Growth(object):
         v0form = (df.inner(v0,tv0)
         v0form = (df.inner(v0,tv0)
                 + self.elasticity * df.inner(u0.dx(0),tv0.dx(0))
                 + self.elasticity * df.inner(u0.dx(0),tv0.dx(0))
                 + self.viscosity * df.inner(v0.dx(0), tv0.dx(0))
                 + self.viscosity * df.inner(v0.dx(0), tv0.dx(0))
-                - self.lamda * df.inner(rho0.dx(0),tv0)
+                - self.lamda * df.inner((rho0/(rho0 + self.saturation_rho)).dx(0),tv0)
                 ) * df.dx
                 ) * df.dx
         bcv0 = df.DirichletBC(v0_function_space, df.Constant(0.0), 'on_boundary')
         bcv0 = df.DirichletBC(v0_function_space, df.Constant(0.0), 'on_boundary')
         df.solve(v0form == 0, v0, bcv0)
         df.solve(v0form == 0, v0, bcv0)
@@ -75,9 +76,9 @@ class Growth(object):
                 + 3/8 * self.elasticity * df.inner(u0.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))
                 + 1/16 * self.elasticity * df.inner(uminus1.dx(0), tv.dx(0))
 
 
-                - 9/16 * self.lamda * df.inner(rho.dx(0), tv)
-                - 3/8 * self.lamda * df.inner(rho0.dx(0), tv)
-                - 1/16 * self.lamda * df.inner(rhominus1.dx(0), tv)
+                - 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))
                 + 9/16 * self.viscosity * df.inner(v.dx(0),tv.dx(0))
                 + 3/8 * self.viscosity * df.inner(v0.dx(0),tv.dx(0))
                 + 3/8 * self.viscosity * df.inner(v0.dx(0),tv.dx(0))