Per-Particle Density

spatialpy/Domain.py

@@ -67,6 +67,7 @@ def __init__(self, numpoints, xlim, ylim, zlim, rho0=1.0, c0=10, P0=10, gravity=
         self.mass = numpy.zeros((numpoints), dtype=float)
         self.type = numpy.zeros((numpoints), dtype=int)
         self.nu = numpy.zeros((numpoints), dtype=float)
+        self.rho = numpy.zeros((numpoints), dtype=float)
         self.fixed = numpy.zeros((numpoints), dtype=bool)
 
         self.rho0 = rho0
@@ -87,7 +88,7 @@ def __str__(self):
         domain_strs.extend(["", "Paritcles", ""])
         for i, vertex in enumerate(self.vertices):
             v_str = f"{pad}{i+1}: {vertex}\n{pad}   Volume:{self.vol[i]}, Mass: {self.mass[i]}, "
-            v_str += f"Type: {self.type[i]}, nu: {self.nu[i]}, Fixed: {self.fixed[i]}"
+            v_str += f"Type: {self.type[i]}, Viscosity: {self.nu[i]}, Density: {self.rho[i]}, Fixed: {self.fixed[i]}"
             domain_strs.append(v_str)
         if self.triangles is not None:
             domain_strs.extend(["", "Triangles", ""])
@@ -103,7 +104,7 @@ def __str__(self):
     def _ipython_display_(self, use_matplotlib=False):
         self.plot_types(width="auto", height="auto", use_matplotlib=use_matplotlib)
 
-    def add_point(self, x, vol, mass, type, nu, fixed):
+    def add_point(self, x, vol, mass, type, nu, fixed, rho=None):
         """ Add a single point particle to the domain space.
 
             :param x: Spatial coordinate vertices of point to be added
@@ -123,12 +124,22 @@ def add_point(self, x, vol, mass, type, nu, fixed):
 
             :param fixed: True if particle is spatially fixed, else False
             :type fixed: bool
+
+            :param rho: Default density of particle to be created
+            :type rho: float
         """
 
+        if vol < 0:
+            raise DomainError("Volume must be a positive value.")
+
+        if rho is None:
+            rho = mass/vol
+
         self.vol = numpy.append(self.vol, vol)
         self.mass = numpy.append(self.mass, mass)
         self.type = numpy.append(self.type, type)
         self.nu = numpy.append(self.nu, nu)
+        self.rho = numpy.append(self.rho, rho)
         self.fixed = numpy.append(self.fixed, fixed)
 
         self.vertices = numpy.append(self.vertices, [x], axis=0)
@@ -458,8 +469,12 @@ def read_xml_mesh(cls, filename):
             obj.tetrahedrons[ int(c.attrib['index']),3] = int(c.attrib['v3'])
         # volume
         obj.calculate_vol()
+        if not numpy.count_nonzero(obj.vol):
+            raise DomainError("Paritcles cannot have 0 volume")
         # set Mass equal to the volume
         obj.mass = obj.vol
+        # Calculate density
+        obj.rho = obj.mass/obj.vol
         # return model ref
         return obj
 
@@ -488,8 +503,12 @@ def import_meshio_object(cls, mesh_obj):
             obj.tetrahedrons = mesh_obj.cells['tetra']
         # volume
         obj.calculate_vol()
+        if not numpy.count_nonzero(obj.vol):
+            raise DomainError("Paritcles cannot have 0 volume")
         # set Mass equal to the volume
         obj.mass = obj.vol
+        # Calculate density
+        obj.rho = obj.mass/obj.vol
         # return model ref
         return obj
 
@@ -507,11 +526,6 @@ def read_msh_file(cls, filename):
             import pygmsh
         except ImportError as e:
             raise DomainError("The python package 'pygmsh' is not installed.")
-       # try:
-       #     _ = pygmsh.get_gmsh_major_version()
-       # except FileNotFoundError as e:
-       #     raise DomainError("The command line program 'gmsh' is not installed or is not found in the current PATH")
-
         try:
             import meshio
         except ImportError as e:
@@ -553,16 +567,18 @@ def read_stochss_domain(cls, filename):
                         rho0=domain['rho_0'], c0=domain['c_0'], P0=domain['p_0'], gravity=domain['gravity'])
 
             for particle in domain['particles']:
+                # StochSS backward compatability check for rho
+                rho = None if "rho" not in particle.keys() else particle['rho']
                 obj.add_point(particle['point'], particle['volume'], particle['mass'],
-                               particle['type'], particle['nu'], particle['fixed'])
+                              particle['type'], particle['nu'], particle['fixed'], rho=rho)
 
             return obj
         except KeyError as e:
             raise DomainError("The file is not a StochSS Domain (.domn) or a StochSS Spatial Model (.smdl).")
 
 
     @classmethod
-    def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=1.0, fixed=False, **kwargs):
+    def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=1.0, rho=None, fixed=False, **kwargs):
         """ Create a filled 3D domain
 
             :param xlim: highest and lowest coordinate in the x dimension
@@ -583,16 +599,19 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=
             :param nz: number of particle spacing in the z dimension
             :type nz: int
 
-            :param type_id: default type ID of particles created to be created. Defaults to 1
+            :param type_id: default type ID of particles to be created. Defaults to 1
             :type type_id: int
 
-            :param mass: default mass of particles created to be created. Defaults to 1.0
+            :param mass: default mass of particles to be created. Defaults to 1.0
             :type mass: float
 
-            :param nu: default viscosity of particles created to be created. Defaults to 1.0
+            :param nu: default viscosity of particles to be created. Defaults to 1.0
             :type nu: float
 
-            :param fixed: spatially fixed flag of particles created to be created. Defaults to false.
+            :param rho: default density of particles to be created.
+            :type rho:
+
+            :param fixed: spatially fixed flag of particles to be created. Defaults to false.
             :type fixed: bool
 
             :param rho0: background density for the system. Defaults to 1.0
@@ -602,7 +621,7 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=
             :type c0: float
 
             :param P0: background pressure for the system. Defaults to 10
-            :type p0: float
+            :type P0: float
 
             :rtype: spatialpy.Domain.Domain
         """
@@ -616,24 +635,30 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=
         z_list = numpy.linspace(zlim[0],zlim[1],nz)
         ndx = 0
         totalvolume = (xlim[1] - xlim[0]) * (ylim[1] - ylim[0]) * (zlim[1] - zlim[0])
+        vol = totalvolume / numberparticles
+        if vol < 0:
+            raise DomainError("Paritcles cannot have 0 volume")
         for x in x_list:
             for y in y_list:
                 for z in z_list:
-                    obj.vol[ndx] = totalvolume / numberparticles
+                    if rho is None:
+                        rho = mass / vol
+                    obj.vol[ndx] = vol
                     obj.vertices[ndx,0] = x
                     obj.vertices[ndx,1] = y
                     obj.vertices[ndx,2] = z
                     obj.type[ndx] = type_id
                     obj.mass[ndx] = mass
                     obj.nu[ndx] = nu
+                    obj.rho[ndx] = rho
                     obj.fixed[ndx] = fixed
                     ndx+=1
 
         # return model ref
         return obj
 
     @classmethod
-    def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, fixed=False, **kwargs):
+    def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, rho=None, fixed=False, **kwargs):
         """ Create a filled 2D domain
 
             :param xlim: highest and lowest coordinate in the x dimension
@@ -648,16 +673,19 @@ def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, fixed
             :param ny: number of particle spacing in the y dimension
             :type ny: int
 
-            :param type_id: default type ID of particles created to be created. Defaults to 1
+            :param type_id: default type ID of particles to be created. Defaults to 1
             :type type_id: int
 
-            :param mass: default mass of particles created to be created. Defaults to 1.0
+            :param mass: default mass of particles to be created. Defaults to 1.0
             :type mass: float
 
-            :param nu: default viscosity of particles created to be created. Defaults to 1.0
+            :param nu: default viscosity of particles to be created. Defaults to 1.0
             :type nu: float
 
-            :param fixed: spatially fixed flag of particles created to be created. Defaults to false.
+            :param rho: default density of particles to be created.
+            :type rho:
+
+            :param fixed: spatially fixed flag of particles to be created. Defaults to false.
             :type fixed: bool
 
             :param rho0: background density for the system. Defaults to 1.0
@@ -667,7 +695,7 @@ def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, fixed
             :type c0: float
 
             :param P0: background pressure for the system. Defaults to 10
-            :type p0: float
+            :type P0: float
 
             :rtype: spatialpy.Domain.Domain
         """
@@ -680,16 +708,21 @@ def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, fixed
         y_list = numpy.linspace(ylim[0],ylim[1],ny)
         ndx = 0
         totalvolume = (xlim[1] - xlim[0]) * (ylim[1] - ylim[0])
-        #print("totalvolume",totalvolume)
+        vol = totalvolume / numberparticles
+        if vol < 0:
+            raise DomainError("Paritcles cannot have 0 volume")
         for x in x_list:
             for y in y_list:
-                obj.vol[ndx] = totalvolume / numberparticles
+                if rho is None:
+                    rho = mass / vol
+                obj.vol[ndx] = vol
                 obj.vertices[ndx,0] = x
                 obj.vertices[ndx,1] = y
                 obj.vertices[ndx,2] = 0.0
                 obj.type[ndx] = type_id
                 obj.mass[ndx] = mass
                 obj.nu[ndx] = nu
+                obj.rho[ndx] = rho
                 obj.fixed[ndx] = fixed
                 ndx+=1
 

spatialpy/Model.py

@@ -237,7 +237,7 @@ def timespan(self, time_span, timestep_size=None):
 
 
 
-    def set_type(self, geometry_ivar, type_id, mass=None, nu=None, fixed=False):
+    def set_type(self, geometry_ivar, type_id, mass=None, nu=None, rho=None, fixed=False):
         """ Add a type definition to the model.  By default, all regions are set to
         type 0. Returns the number of domain points that were tagged with this type_id
 
@@ -248,6 +248,8 @@ def set_type(self, geometry_ivar, type_id, mass=None, nu=None, fixed=False):
             :type type_id: int
             :param mass: The mass of each particle in the type
             :type mass: float
+            :param rho: The density of each particle in the type
+            :type rho: float
             :param nu: The viscosity of each particle in the type
             :type nu: float
             :param fixed: Are the particles in this type immobile
@@ -267,9 +269,11 @@ def set_type(self, geometry_ivar, type_id, mass=None, nu=None, fixed=False):
         for v_ndx in range(self.domain.get_num_voxels()):
             if geometry_ivar.inside( self.domain.coordinates()[v_ndx,:], on_boundary[v_ndx]):
                 self.domain.type[v_ndx] = type_id
-                if (mass is not None):
+                if mass is not None:
                     self.domain.mass[v_ndx] = mass
-                if (nu is not None):
+                if rho is not None:
+                    self.domain.rho[v_ndx] = rho
+                if nu is not None:
                     self.domain.nu[v_ndx] = nu
                 self.domain.fixed[v_ndx] = fixed
                 count +=1

spatialpy/Solver.py

@@ -433,7 +433,7 @@ def __create_propensity_file(self, file_name=None):
             init_particles += "init_create_particle(sys,id++,{0},{1},{2},{3},{4},{5},{6},{7},{8});".format(
                 self.model.domain.coordinates()[i,0],self.model.domain.coordinates()[i,1],self.model.domain.coordinates()[i,2],
                 self.model.domain.type[i],self.model.domain.nu[i],self.model.domain.mass[i],
-                (self.model.domain.mass[i] / self.model.domain.vol[i]),int(self.model.domain.fixed[i]),num_chem_species )+"\n"
+                self.model.domain.rho[i],int(self.model.domain.fixed[i]),num_chem_species )+"\n"
         propfilestr = propfilestr.replace("__INIT_PARTICLES__", init_particles)
 
         # process initial conditions here