diff --git a/examples/dem.py b/examples/dem.py
index 550c3902218c6f1aea246d5d89c42b67377a426a..a64139bef97fd8a2353e11da8069b8bf28a9cd10 100644
--- a/examples/dem.py
+++ b/examples/dem.py
@@ -4,8 +4,13 @@ import sys
def update_mass_and_inertia(i):
- mass[i] = (4.0 / 3.0) * pi * radius[i] * radius[i] * radius[i] * densityParticle_SI
- inv_inertia[i] = diagonal_matrix(1.0 / (0.4 * mass[i] * radius[i] * radius[i]))
+ if is_sphere(i):
+ mass[i] = ((4.0 / 3.0) * pi) * radius[i] * radius[i] * radius[i] * densityParticle_SI
+ inv_inertia[i] = inversed(diagonal_matrix(0.4 * mass[i] * radius[i] * radius[i]))
+
+ else:
+ mass[i] = infinity
+ inv_inertia[i] = 0.0
def linear_spring_dashpot(i, j):
@@ -161,7 +166,8 @@ psim.read_particle_data(
pairs.halfspace())
psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI,
- 'pi': math.pi })
+ 'pi': math.pi,
+ 'infinity': math.inf })
psim.build_neighbor_lists(linkedCellWidth + skin)
psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index b31e997c1b4611bd3fe8613d4c7aebb0e67ac222..d38b06134d61acadef76ee93b5c77f1780b707a0 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -1,3 +1,4 @@
+import math
from pairs.ir.assign import Assign
from pairs.ir.atomic import AtomicAdd
from pairs.ir.arrays import Array, ArrayAccess, DeclareStaticArray, RegisterArray, ReallocArray
@@ -56,6 +57,7 @@ class CGen:
if self.target.is_gpu():
self.print("#define PAIRS_TARGET_CUDA")
+ self.print("#include <limits.h>")
self.print("#include <math.h>")
self.print("#include <stdbool.h>")
self.print("#include <stdio.h>")
@@ -806,6 +808,9 @@ class CGen:
assert index is not None, "Index must be set for non-scalar literals."
return ast_node.value[index]
+ if isinstance(ast_node.value, float) and math.isinf(ast_node.value):
+ return "std::numeric_limits<double>::infinity()"
+
return ast_node.value
if isinstance(ast_node, MathFunction):
diff --git a/src/pairs/mapping/keywords.py b/src/pairs/mapping/keywords.py
index cf28e9cff2cb634ca1a31ac5a1b9626315ca9b8a..c5e84597be7691a2addb247f24edc6cb6148749e 100644
--- a/src/pairs/mapping/keywords.py
+++ b/src/pairs/mapping/keywords.py
@@ -5,9 +5,11 @@ from pairs.ir.loops import Continue
from pairs.ir.math import Abs, Cos, Sin, Sqrt
from pairs.ir.matrices import Matrix
from pairs.ir.quaternions import Quaternion
+from pairs.ir.scalars import ScalarOp
from pairs.ir.select import Select
from pairs.ir.types import Types
from pairs.ir.vectors import Vector, ZeroVector
+from pairs.sim.shapes import Shapes
class Keywords:
@@ -27,6 +29,18 @@ class Keywords:
method = self.get_method(f"keyword_{keyword}")
return method is not None
+ def keyword_is_sphere(self, args):
+ assert len(args) == 1, "is_sphere() keyword requires one parameter."
+ particle_id = args[0]
+ assert particle_id.type() == Types.Int32, "Particle ID must be an integer."
+ return ScalarOp.cmp(self.sim.particle_shape[particle_id], Shapes.Sphere)
+
+ def keyword_is_halfspace(self, args):
+ assert len(args) == 1, "is_sphere() keyword requires one parameter."
+ particle_id = args[0]
+ assert particle_id.type() == Types.Int32, "Particle ID must be an integer."
+ return ScalarOp.cmp(self.sim.particle_shape[particle_id], Shapes.Halfspace)
+
def keyword_select(self, args):
assert len(args) == 3, "select() keyword requires three parameters!"
return Select(self.sim, args[0], args[1], args[2])
@@ -84,23 +98,42 @@ class Keywords:
#return vector * inv_length
def keyword_squared_length(self, args):
- assert len(args) == 1, "length() keyword requires one parameter!"
+ assert len(args) == 1, "length() keyword requires one parameter."
vector = args[0]
- assert vector.type() == Types.Vector, "length(): Argument must be a vector!"
+ assert vector.type() == Types.Vector, "length(): Argument must be a vector."
return sum([vector[d] * vector[d] for d in range(self.sim.ndims())])
def keyword_zero_vector(self, args):
- assert len(args) == 0, "zero_vector() keyword requires no parameter!"
+ assert len(args) == 0, "zero_vector() keyword requires no parameter."
return ZeroVector(self.sim)
def keyword_transposed(self, args):
- assert len(args) == 1, "transposed() keyword requires one parameter!"
+ assert len(args) == 1, "transposed() keyword requires one parameter."
matrix = args[0]
- assert matrix.type() == Types.Matrix, "tranposed(): Argument must be a matrix!"
+ assert matrix.type() == Types.Matrix, "tranposed(): Argument must be a matrix."
return Matrix(self.sim, [ matrix[0], matrix[3], matrix[6],
matrix[1], matrix[4], matrix[7],
matrix[2], matrix[5], matrix[8] ])
+ def keyword_inversed(self, args):
+ assert len(args) == 1, "inversed() keyword requires one parameter."
+ matrix = args[0]
+ assert matrix.type() == Types.Matrix, "inversed(): Argument must be a matrix."
+ det = matrix[0] * ((matrix[4] * matrix[8]) - (matrix[7] * matrix[5])) + \
+ matrix[1] * ((matrix[5] * matrix[6]) - (matrix[8] * matrix[3])) + \
+ matrix[2] * ((matrix[3] * matrix[7]) - (matrix[4] * matrix[6]))
+
+ inv_det = 1.0 / det
+ return Matrix(self.sim, [ inv_det * ((matrix[4] * matrix[8]) - (matrix[5] * matrix[7])),
+ inv_det * ((matrix[7] * matrix[2]) - (matrix[8] * matrix[1])),
+ inv_det * ((matrix[1] * matrix[5]) - (matrix[2] * matrix[4])),
+ inv_det * ((matrix[5] * matrix[6]) - (matrix[3] * matrix[8])),
+ inv_det * ((matrix[8] * matrix[0]) - (matrix[6] * matrix[2])),
+ inv_det * ((matrix[2] * matrix[3]) - (matrix[0] * matrix[5])),
+ inv_det * ((matrix[3] * matrix[7]) - (matrix[4] * matrix[6])),
+ inv_det * ((matrix[6] * matrix[1]) - (matrix[7] * matrix[0])),
+ inv_det * ((matrix[0] * matrix[4]) - (matrix[1] * matrix[3])) ])
+
def keyword_diagonal_matrix(self, args):
assert len(args) == 1, "diagonal_matrix() keyword requires one parameter!"
value = args[0]