diff --git a/spatialpy/Domain.py b/spatialpy/Domain.py index 1ca08b2f..df1dc473 100644 --- a/spatialpy/Domain.py +++ b/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,8 +567,10 @@ 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: @@ -562,7 +578,7 @@ def read_stochss_domain(cls, filename): @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,16 +635,22 @@ 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 @@ -633,7 +658,7 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu= 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 diff --git a/spatialpy/Model.py b/spatialpy/Model.py index 58fd19e1..6b33bd10 100644 --- a/spatialpy/Model.py +++ b/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 diff --git a/spatialpy/Solver.py b/spatialpy/Solver.py index a3ff2929..80bd2a1f 100644 --- a/spatialpy/Solver.py +++ b/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