Skip to content
Snippets Groups Projects
Commit b026d7ec authored by Rafael Ravedutti's avatar Rafael Ravedutti
Browse files

Small fixes

parent 558978a8
Branches
Tags
1 merge request!1Implement DEM and many other features
...@@ -4,8 +4,13 @@ import sys ...@@ -4,8 +4,13 @@ import sys
def update_mass_and_inertia(i): def update_mass_and_inertia(i):
mass[i] = (4.0 / 3.0) * pi * radius[i] * radius[i] * radius[i] * densityParticle_SI if is_sphere(i):
inv_inertia[i] = diagonal_matrix(1.0 / (0.4 * mass[i] * radius[i] * radius[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): def linear_spring_dashpot(i, j):
...@@ -161,7 +166,8 @@ psim.read_particle_data( ...@@ -161,7 +166,8 @@ psim.read_particle_data(
pairs.halfspace()) pairs.halfspace())
psim.setup(update_mass_and_inertia, {'densityParticle_SI': densityParticle_SI, 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.build_neighbor_lists(linkedCellWidth + skin)
psim.vtk_output(f"output/dem_{target}", frequency=visSpacing) psim.vtk_output(f"output/dem_{target}", frequency=visSpacing)
......
import math
from pairs.ir.assign import Assign from pairs.ir.assign import Assign
from pairs.ir.atomic import AtomicAdd from pairs.ir.atomic import AtomicAdd
from pairs.ir.arrays import Array, ArrayAccess, DeclareStaticArray, RegisterArray, ReallocArray from pairs.ir.arrays import Array, ArrayAccess, DeclareStaticArray, RegisterArray, ReallocArray
...@@ -56,6 +57,7 @@ class CGen: ...@@ -56,6 +57,7 @@ class CGen:
if self.target.is_gpu(): if self.target.is_gpu():
self.print("#define PAIRS_TARGET_CUDA") self.print("#define PAIRS_TARGET_CUDA")
self.print("#include <limits.h>")
self.print("#include <math.h>") self.print("#include <math.h>")
self.print("#include <stdbool.h>") self.print("#include <stdbool.h>")
self.print("#include <stdio.h>") self.print("#include <stdio.h>")
...@@ -806,6 +808,9 @@ class CGen: ...@@ -806,6 +808,9 @@ class CGen:
assert index is not None, "Index must be set for non-scalar literals." assert index is not None, "Index must be set for non-scalar literals."
return ast_node.value[index] 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 return ast_node.value
if isinstance(ast_node, MathFunction): if isinstance(ast_node, MathFunction):
......
...@@ -5,9 +5,11 @@ from pairs.ir.loops import Continue ...@@ -5,9 +5,11 @@ from pairs.ir.loops import Continue
from pairs.ir.math import Abs, Cos, Sin, Sqrt from pairs.ir.math import Abs, Cos, Sin, Sqrt
from pairs.ir.matrices import Matrix from pairs.ir.matrices import Matrix
from pairs.ir.quaternions import Quaternion from pairs.ir.quaternions import Quaternion
from pairs.ir.scalars import ScalarOp
from pairs.ir.select import Select from pairs.ir.select import Select
from pairs.ir.types import Types from pairs.ir.types import Types
from pairs.ir.vectors import Vector, ZeroVector from pairs.ir.vectors import Vector, ZeroVector
from pairs.sim.shapes import Shapes
class Keywords: class Keywords:
...@@ -27,6 +29,18 @@ class Keywords: ...@@ -27,6 +29,18 @@ class Keywords:
method = self.get_method(f"keyword_{keyword}") method = self.get_method(f"keyword_{keyword}")
return method is not None 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): def keyword_select(self, args):
assert len(args) == 3, "select() keyword requires three parameters!" assert len(args) == 3, "select() keyword requires three parameters!"
return Select(self.sim, args[0], args[1], args[2]) return Select(self.sim, args[0], args[1], args[2])
...@@ -84,23 +98,42 @@ class Keywords: ...@@ -84,23 +98,42 @@ class Keywords:
#return vector * inv_length #return vector * inv_length
def keyword_squared_length(self, args): 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] 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())]) return sum([vector[d] * vector[d] for d in range(self.sim.ndims())])
def keyword_zero_vector(self, args): 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) return ZeroVector(self.sim)
def keyword_transposed(self, args): 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] 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], return Matrix(self.sim, [ matrix[0], matrix[3], matrix[6],
matrix[1], matrix[4], matrix[7], matrix[1], matrix[4], matrix[7],
matrix[2], matrix[5], matrix[8] ]) 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): def keyword_diagonal_matrix(self, args):
assert len(args) == 1, "diagonal_matrix() keyword requires one parameter!" assert len(args) == 1, "diagonal_matrix() keyword requires one parameter!"
value = args[0] value = args[0]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment