added boundary condition and more active stress flags

growing_domain/full_model/growing_domain.py

@@ -47,7 +47,12 @@ class Growing_Domain(object):
         self.function_space = df.FunctionSpace(self.mesh, mixed_element)
         self.function_space = df.FunctionSpace(self.mesh, mixed_element)
 
 
         # dirichlet boundaries for rho
         # dirichlet boundaries for rho
-        self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
+        if self.dirichlet_boundary_rho:
+            if self.zero_rho_boundary:
+                self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(0.0), 'on_boundary')
+            elif self.rho0_boundary:
+                self.bc = df.DirichletBC(self.function_space.sub(2), df.Constant(self.average_rho), 'on_boundary')
+
 
 
         # 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
@@ -57,21 +62,20 @@ class Growing_Domain(object):
         return df.inner(df.div(vel * conc), tconc)
         return df.inner(df.div(vel * conc), tconc)
 
 
     def active_stress(self, rho, c, v):
     def active_stress(self, rho, c, v):
-        if self.morphogen_mediated_active_stress:
-            return self.lamda * rho/(rho + self.saturation_rho) * (c-self.average_c) * df.Identity(self.dimension)
-        if self.rho_t_dependent_active_stress:
-            if self.density_dependent_turnover_rho:
-                rho_t_dependence = (self.diffusion_rho*df.div(df.nabla_grad(rho))
-                                    - self.turnover_rho * rho * (rho - self.average_rho)
-                                    - df.div(rho * v))
-            else:
-                rho_t_dependence = (self.diffusion_rho*df.div(df.nabla_grad(rho))
-                    - self.turnover_rho * (rho - self.average_rho)
-                    - df.div(rho * v))
-
-            return self.lamda * rho_t_dependence/(rho + self.saturation_rho) * df.Identity(self.dimension)
-        else:
+        if self.saturating_density_active_stress:
             return self.lamda * rho/(rho + self.saturation_rho)* df.Identity(self.dimension)
             return self.lamda * rho/(rho + self.saturation_rho)* df.Identity(self.dimension)
+        
+        elif self.difference_in_morphogen_active_stress:
+            return self.lamda * (c-self.average_c) * df.Identity(self.dimension)
+        
+        elif self.difference_in_density_active_stress:
+            return self.lamda * (c-self.average_rho) * df.Identity(self.dimension)
+        
+        elif self.saturating_density_and_difference_in_density_active_stress:
+            return self.lamda * (rho*(rho-self.average_rho))/(rho + self.saturation_rho) * df.Identity(self.dimension)
+        
+        elif self.saturating_density_and_difference_in_morphogen_active_stress:
+            return self.lamda * (rho*(c-self.average_c))/(rho + self.saturation_rho) * df.Identity(self.dimension)
 
 
     def epsilon(self, v):
     def epsilon(self, v):
         return df.sym(df.nabla_grad(v))
         return df.sym(df.nabla_grad(v))
@@ -92,22 +96,29 @@ class Growing_Domain(object):
     def stress(self, u, v, rho, c):
     def stress(self, u, v, rho, c):
         return (self.passive_stress(u, v) + self.active_stress(rho, c, v))
         return (self.passive_stress(u, v) + self.active_stress(rho, c, v))
         
         
-    def reaction_rho(self, rho, trho):
+    def reaction_rho(self, rho, trho, c):
+
         if self.morphogen_mediated_turnover_rho:
         if self.morphogen_mediated_turnover_rho:
-            return self.average_c * df.inner(rho - self.average_rho, trho)
+            return self.turnover_rho * df.inner((c - self.average_c/2)*(rho - self.average_rho), trho)
+
         elif self.morphogen_mediated_average_rho:
         elif self.morphogen_mediated_average_rho:
             return self.turnover_rho * df.inner(rho - self.average_c, trho)
             return self.turnover_rho * df.inner(rho - self.average_c, trho)
-        elif self.density_dependent_turnover_rho:        # to model that cells are born or die only where the density is non-zero
-            return self.turnover_rho * df.inner(rho * (rho - self.average_rho), trho)
+
+        elif self.density_dependent_turnover_rho:        
+            # to model that cells are born or die only where the density is non-zero
+
+            # tanh approximation: Theta = 1+ tanh = 1 + rho - rho^3/3! + . . . 
+            return self.turnover_rho * df.inner((1 + rho) * (rho - self.average_rho), trho)  
+
         else:
         else:
             return self.turnover_rho * df.inner(rho - self.average_rho, trho)
             return self.turnover_rho * df.inner(rho - self.average_rho, trho)
 
 
     def reaction_c(self, c, tc):
     def reaction_c(self, c, tc):
         return self.turnover_c * df.inner(c - self.average_c, tc)
         return self.turnover_c * df.inner(c - self.average_c, tc)
     
     
-    def diffusion_reaction_rho(self, rho, trho):
+    def diffusion_reaction_rho(self, rho, trho, c):
         return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho)) 
         return (self.diffusion_rho * df.inner(df.nabla_grad(rho), df.nabla_grad(trho)) 
-                + self.reaction_rho(rho, trho))
+                + self.reaction_rho(rho, trho, c))
 
 
     def diffusion_reaction_c(self, c, tc):
     def diffusion_reaction_c(self, c, tc):
         return (self.diffusion_c * df.inner(df.nabla_grad(c), df.nabla_grad(tc)) 
         return (self.diffusion_c * df.inner(df.nabla_grad(c), df.nabla_grad(tc)) 
@@ -125,14 +136,23 @@ class Growing_Domain(object):
         u0.vector()[:] = u0.vector()[:] + noise_u 
         u0.vector()[:] = u0.vector()[:] + noise_u 
 
 
         # ic for rho
         # ic for rho
+
         r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
         r2 = "+".join(['x['+str(i)+']*x['+str(i)+']' for i in range(self.dimension)])
 
 
         if self.dimension==1:
         if self.dimension==1:
             rho0='1-cos(2*pi*x[0]/R)'
             rho0='1-cos(2*pi*x[0]/R)'
 
 
         else:
         else:
+            # for rho=rho_0 boundary condition
+            # rho0 = "%s" % str(self.average_rho)
+            # rho0 = '(a/2)*(1 + tanh((sqrt(%s) - r0)/w))' %r2
+
+            # grad rho=0 boundary condition
+            rho0 = '1 + 0.1 * a * cos(2.0*pi*%s/R)' %r2 
+
+            # rho = 0 boundary condition
             # circular domain --> circular domain   
             # circular domain --> circular domain   
-            rho0 = '(a/2)*(1-tanh((sqrt(%s) - r0)/w))' %r2
+            # rho0 = '(a/2)*(1 - tanh((sqrt(%s) - r0)/w))' %r2
             # rho0 = 'a*(1-tanh((sqrt(%s) - r0)/w))' %r2   # grows to a different size
             # rho0 = 'a*(1-tanh((sqrt(%s) - r0)/w))' %r2   # grows to a different size
 
 
             # circular domain --> elliptical domain (w=0.19, lambda = -10, k=1)  
             # circular domain --> elliptical domain (w=0.19, lambda = -10, k=1)  
@@ -144,17 +164,18 @@ class Growing_Domain(object):
             # circular --> kidney (w=0.8, r0=0.5, k=0.1, lambda=-30)
             # circular --> kidney (w=0.8, r0=0.5, k=0.1, lambda=-30)
             # rho0 = 'x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2 
             # rho0 = 'x[0]*(1-tanh((sqrt(%s) - r0)/w))' %r2 
 
 
-        rho0 = df.interpolate(df.Expression(rho0, pi=np.pi, R=self.system_size, a=self.average_rho, r0=0.5, w=0.3, degree=1), self.function_space.sub(2).collapse())
+        rho0 = df.interpolate(df.Expression(rho0, pi=np.pi, R=self.system_size, a=self.average_rho, r0=0.5, w=0.19, degree=1), self.function_space.sub(2).collapse())
         # add noise
         # add noise
         noise_rho = self.noise_level * (2 * np.random.random(rho0.vector().size()) - 1)
         noise_rho = self.noise_level * (2 * np.random.random(rho0.vector().size()) - 1)
         rho0.vector()[:] = rho0.vector()[:] + noise_rho
         rho0.vector()[:] = rho0.vector()[:] + noise_rho
 
 
         # ic for c
         # ic for c
         if self.morphogen:
         if self.morphogen:
+            # c0 = "%s" % str(5.0)
             # c0 = '1 + 0.1*a*cos(2.0*pi*sqrt(%s)/R)' %r2  # shrinks and stops  
             # c0 = '1 + 0.1*a*cos(2.0*pi*sqrt(%s)/R)' %r2  # shrinks and stops  
             c0 = '1 + 0.1*a*cos(2.0*pi*%s/R)' %r2          # grows and stops    
             c0 = '1 + 0.1*a*cos(2.0*pi*%s/R)' %r2          # grows and stops    
             # c0 = '1 + 0.1*a*cos(2.0*pi*x[0]/R)'          # grows very slightly, shrinks even more slightly, changes shape 
             # c0 = '1 + 0.1*a*cos(2.0*pi*x[0]/R)'          # grows very slightly, shrinks even more slightly, changes shape 
-
+            
             c0 = df.interpolate(df.Expression(c0, a = self.average_c, R=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse())
             c0 = df.interpolate(df.Expression(c0, a = self.average_c, R=self.system_size, pi=np.pi, degree=1), self.function_space.sub(3).collapse())
             # add noise
             # add noise
             noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
             noise_c = self.noise_level * (2*np.random.random(c0.vector().size())-1)
@@ -195,7 +216,7 @@ class Growing_Domain(object):
 
 
         rhoform = (df.inner((rho - rho0)/self.timestep, trho)
         rhoform = (df.inner((rho - rho0)/self.timestep, trho)
                 + self.advection(rho0, v0, trho)
                 + self.advection(rho0, v0, trho)
-                + self.diffusion_reaction_rho(rho, trho)
+                + self.diffusion_reaction_rho(rho, trho, c)
                 )
                 )
         
         
         if self.morphogen:
         if self.morphogen:
@@ -236,7 +257,12 @@ class Growing_Domain(object):
 
 
         for steps in progressbar.progressbar(range(1, maxsteps + 1)):
         for steps in progressbar.progressbar(range(1, maxsteps + 1)):
             # solve
             # solve
-            df.solve(self.form == 0, self.f, self.bc)
+            if self.dirichlet_boundary_rho:
+                df.solve(self.form == 0, self.f, self.bc)
+            elif self.neumann_boundary_rho:
+                if self.zero_grad_rho_boundary:
+                    df.solve(self.form==0, self.f)
+
             # update
             # update
             self.f0.assign(self.f)
             self.f0.assign(self.f)
             if self.morphogen:
             if self.morphogen:
@@ -257,6 +283,7 @@ class Growing_Domain(object):
             df.ALE.move(self.mesh, dr)
             df.ALE.move(self.mesh, dr)
             # update time
             # update time
             self.time += self.timestep
             self.time += self.timestep
+
         self.uFile.close()
         self.uFile.close()
         self.vFile.close()
         self.vFile.close()
         self.rhoFile.close()    # def save_data_uvrho(self):
         self.rhoFile.close()    # def save_data_uvrho(self):

growing_domain/full_model/parameters.json

@@ -5,32 +5,42 @@
     "resolution" : 100,
     "resolution" : 100,
     "system_size" : 1,
     "system_size" : 1,
     "meshfile" : "disk.xml.gz",
     "meshfile" : "disk.xml.gz",
-    "morphogen_mediated_active_stress": false,
-    "rho_t_dependent_active_stress" : true,
+
+    "dirichlet_boundary_rho": false,
+    "zero_rho_boundary":false,
+    "rho0_boundary": false,
+
+    "neumann_boundary_rho": true,
+    "zero_grad_rho_boundary": true,
+
+    "difference_in_density_active_stress": false,
+    "difference_in_morphogen_active_stress": false,
+    "saturating_density_active_stress": true,
+    "saturating_density_and_difference_in_morphogen_active_stress": false,
+    "saturating_density_and_difference_in_density_active_stress":false,
+
+    "morphogen_mediated_turnover_rho": false,
+    "density_dependent_turnover_rho": true,
+
+    "morphogen_mediated_average_rho": false,
 
 
     "bulk_elasticity" : 1.0,
     "bulk_elasticity" : 1.0,
     "shear_elasticity": 1.0,
     "shear_elasticity": 1.0,
     "shear_viscosity" : 1.0,
     "shear_viscosity" : 1.0,
     "bulk_viscosity" : 1.0,
     "bulk_viscosity" : 1.0,
     "friction" : 1.0,
     "friction" : 1.0,
-    "lamda": -0.1,
-
+    "lamda": -1.0,
     "diffusion_rho" : 1.0,
     "diffusion_rho" : 1.0,
-    "morphogen_mediated_turnover_rho": false,
-    "density_dependent_turnover_rho": false,
-    "morphogen_mediated_average_rho": false,
     "turnover_rho" : 1.0,
     "turnover_rho" : 1.0,
     "average_rho" : 1.0,
     "average_rho" : 1.0,
     "saturation_rho" : 1.0,
     "saturation_rho" : 1.0,
-
     "diffusion_c" : 1.0,
     "diffusion_c" : 1.0,
     "turnover_c" : 1.0,
     "turnover_c" : 1.0,
     "average_c" : 1.0,
     "average_c" : 1.0,
-
     "noise_level": 0.0,
     "noise_level": 0.0,
 
 
     "timestep" : 0.1,
     "timestep" : 0.1,
     "savetime": 0.2,
     "savetime": 0.2,
-    "maxtime" : 40
+    "maxtime" : 100.0
 
 
 }
 }
\ No newline at end of file