Skip to content

HOTFIX: Lattice Divide By Zero Error #333

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
seanebum merged 3 commits into from
Dec 22, 2022
Merged

HOTFIX: Lattice Divide By Zero Error #333

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
34eae14
Fixed divide by zero error in spherical lattices.
BryanRumsey Dec 15, 2022
69d80b0
Fixed divide by zero error in cylindrical lattices.
BryanRumsey Dec 15, 2022
740b49e
Fixed import error. Pylint changes.
BryanRumsey Dec 15, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 32 additions & 16 deletions spatialpy/core/lattice.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,15 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
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.
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -309,7 +310,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)))

Expand All @@ -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)
Expand All @@ -339,6 +340,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
Expand Down Expand Up @@ -437,7 +444,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)))

Expand All @@ -448,7 +455,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)
Expand All @@ -464,6 +470,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
Expand Down Expand Up @@ -501,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)

Expand Down Expand Up @@ -587,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)
Expand Down Expand Up @@ -637,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)

Expand Down Expand Up @@ -702,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:
Expand All @@ -719,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)
Expand All @@ -735,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:
Expand Down Expand Up @@ -774,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)

Expand Down Expand Up @@ -814,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('-', '')
Expand Down