From 34eae14a5a2dc8569031570b369063421f1d9fe1 Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 15 Dec 2022 09:28:36 -0500 Subject: [PATCH 1/3] Fixed divide by zero error in spherical lattices. --- spatialpy/core/lattice.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 6b70864c..0a91bb04 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -309,7 +309,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): count = 0 radius = self.radius - while radius >= 0: + while radius > 0: # Calculate the approximate number of particle with the radius approx_rc = int(round((4 * radius ** 2) / ((self.deltas / 2) ** 2))) @@ -339,6 +339,12 @@ def apply(self, domain, geometry, transform=None, **kwargs): domain.add_point(point + self.center, **kwargs) count += 1 radius -= self.deltar + if radius == 0 and geometry.inside((0, 0, 0), False): + point = [0, 0, 0] if transform is None else transform([0, 0, 0]) + if not isinstance(point, numpy.ndarray): + point = numpy.array(point) + domain.add_point(point + self.center, **kwargs) + count += 1 self._update_limits(domain) if 'vol' not in kwargs: offset = len(domain.vertices) - count From 69d80b0a795fca40e45e48eb2b3b7171c492864f Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 15 Dec 2022 09:35:40 -0500 Subject: [PATCH 2/3] Fixed divide by zero error in cylindrical lattices. --- spatialpy/core/lattice.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 0a91bb04..1548bd88 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -443,7 +443,7 @@ def apply(self, domain, geometry, transform=None, **kwargs): xmin = -h_len xmax = h_len radius = self.radius - while radius >= 0: + while radius > 0: # Calculate the approximate number of particle with the radius approx_rc = int(round((2 * radius * self.length) / ((self.deltas / 2) ** 2))) @@ -454,7 +454,6 @@ def apply(self, domain, geometry, transform=None, **kwargs): x = xmin while x <= xmax: - for mtheta in range(m_theta): theta = 2 * numpy.pi * (mtheta + 0.5) / m_theta y = radius * numpy.cos(theta) @@ -470,6 +469,16 @@ def apply(self, domain, geometry, transform=None, **kwargs): count += 1 x += self.deltas radius -= self.deltar + if radius == 0: + x = xmin + while x <= xmax: + if geometry.inside((x, 0, 0), False): + point = [x, 0, 0] if transform is None else transform([x, 0, 0]) + if not isinstance(point, numpy.ndarray): + point = numpy.array(point) + domain.add_point(point + self.center, **kwargs) + count += 1 + x += self.deltas self._update_limits(domain) if 'vol' not in kwargs: offset = len(domain.vertices) - count From 740b49eb8d8ee61ceef63f0a90887ebf1f97d3be Mon Sep 17 00:00:00 2001 From: Bryan Rumsey Date: Thu, 15 Dec 2022 09:43:11 -0500 Subject: [PATCH 3/3] Fixed import error. Pylint changes. --- spatialpy/core/lattice.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/spatialpy/core/lattice.py b/spatialpy/core/lattice.py index 1548bd88..d23ce53d 100644 --- a/spatialpy/core/lattice.py +++ b/spatialpy/core/lattice.py @@ -14,14 +14,15 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . import json -import numpy +import string import xml.etree.ElementTree as ET +import numpy + from spatialpy.core.geometry import Geometry, CombinatoryGeometry from spatialpy.core.spatialpyerror import LatticeError - class Lattice: """ Lattice class provides a method for creating parts of the spatial domain. @@ -102,7 +103,7 @@ class CartesianLattice(Lattice): :param zmin: Minimum z value of the lattice (optional, defaults to xmin). :type zmin: float - + :param zmax: Maximum z value of the lattice (optional, defaults to xmax). :type zmax: float @@ -319,11 +320,11 @@ def apply(self, domain, geometry, transform=None, **kwargs): m_phi = int(round(numpy.pi * radius / d_a)) d_phi = numpy.pi / m_phi d_theta = p_area / d_phi - + for mphi in range(m_phi): phi = numpy.pi * (mphi + 0.5) / m_phi m_theta = int(round(2 * numpy.pi * numpy.sin(phi) / d_phi)) - + for mtheta in range(m_theta): theta = 2 * numpy.pi * mtheta / m_theta x = radius * numpy.cos(theta) * numpy.sin(phi) @@ -516,7 +517,7 @@ class XMLMeshLattice(Lattice): """ XML mesh lattice class provides a method for creating parts of the spatial domain with a mesh defined by a FEniCS/dolfin style XML mesh file. - + :param center: The center point of the lattice. :type center: float[3] | float(3) @@ -602,7 +603,7 @@ def apply(self, domain, *args, transform=None, **kwargs): point = [x, y, z] else: point = transform([x, y, z]) - if type_ids is not None and i in type_ids.keys(): + if type_ids is not None and i in type_ids: kwargs['type_id'] = type_ids[i] if not isinstance(point, numpy.ndarray): point = numpy.array(point) @@ -652,7 +653,7 @@ class MeshIOLattice(Lattice): """ meshio lattice class provides a method for creating parts of the spatial domain with a mesh defined by a Gmsh style .msh mesh file. - + :param center: The center point of the lattice. :type center: float[3] | float(3) @@ -717,7 +718,7 @@ def apply(self, domain, *args, transform=None, **kwargs): mesh = meshio.read(self.filename) else: mesh = self.mesh - + if self.subdomain_file is not None: type_ids = self.__get_types() else: @@ -734,7 +735,7 @@ def apply(self, domain, *args, transform=None, **kwargs): point = [x, y, z] else: point = transform([x, y, z]) - if type_ids is not None and i in type_ids.keys(): + if type_ids is not None and i in type_ids: kwargs['type_id'] = type_ids[i] if not isinstance(point, numpy.ndarray): point = numpy.array(point) @@ -750,7 +751,7 @@ def apply(self, domain, *args, transform=None, **kwargs): for triangle in triangles[0].data: triangle = triangle + num_points domain.triangles = numpy.append(domain.triangles, [triangle], axis=0) - + #tetrahedrons tetras = list(filter(lambda cell: cell.type == "tetra", mesh.cells)) if tetras: @@ -789,7 +790,7 @@ class StochSSLattice(Lattice): """ stochss lattice class provides a method for creating parts of the spatial domain with a domain defined by a stochss style .domn domain file or .smdl model file. - + :param center: The center point of the lattice. :type center: float[3] | float(3) @@ -829,7 +830,7 @@ def apply(self, domain, *args, transform=None, **kwargs): domain.c0 = s_domain['c_0'] domain.P0 = s_domain['p_0'] domain.gravity = s_domain['gravity'] - + type_ids = {} for s_type in s_domain['types']: type_ids[s_type['typeID']] = s_type['name'].replace('-', '')