diff --git a/benchmarks/2d/uniaxial_particle_traction/mpm-particle-traction.toml b/benchmarks/2d/uniaxial_particle_traction/mpm-particle-traction.toml new file mode 100644 index 0000000..f5a45f0 --- /dev/null +++ b/benchmarks/2d/uniaxial_particle_traction/mpm-particle-traction.toml @@ -0,0 +1,68 @@ +# The `meta` group contains top level attributes that govern the +# behaviour of the MPM Solver. +# +# Attributes: +# title: The title of the experiment. This is just for the user's +# reference. +# type: The type of simulation to be used. Allowed values are +# {"MPMExplicit"} +# scheme: The MPM Scheme used for simulation. Allowed values are +# {"usl", "usf"} +# dt: Timestep used in the simulation. +# nsteps: Number of steps to run the simulation for. +[meta] +title = "uniaxial-particle-traction" +type = "MPMExplicit" +dimension = 2 +scheme = "usf" +dt = 0.001 +nsteps = 1000 +velocity_update = true + +[output] +format = "npz" +folder = "results/" +step_frequency = 10 + +[mesh] +# type = "file" +# file = "mesh-1d.txt" +# boundary_nodes = "boundary-1d.txt" +# particle_element_ids = "particles-elements.txt" +type = "generator" +nelements = [3, 1] +element_length = [0.1, 0.1] +particle_element_ids = [0] +element = "Quadrilateral4Node" + +[[mesh.constraints]] +node_ids = [0, 4] +dir = 0 +velocity = 0.0 + +[[materials]] +id = 0 +density = 1000 +poisson_ratio = 0 +youngs_modulus = 1000000 +type = "LinearElastic" + +[[particles]] +file = "particles-2d-particle-traction.json" +material_id = 0 +init_velocity = 0.0 + +[external_loading] +gravity = [0, 0] + +[[external_loading.particle_surface_traction]] +pset = [0] +pids = [[5, 11]] +math_function_id = 0 +dir = 0 +traction = 1 + +[[math_functions]] +type = "Linear" +xvalues = [0.0, 0.5, 1.0] +fxvalues = [0.0, 1.0, 1.0] diff --git a/benchmarks/2d/uniaxial_particle_traction/particles-2d-particle-traction.json b/benchmarks/2d/uniaxial_particle_traction/particles-2d-particle-traction.json new file mode 100644 index 0000000..f0143b8 --- /dev/null +++ b/benchmarks/2d/uniaxial_particle_traction/particles-2d-particle-traction.json @@ -0,0 +1,74 @@ +[ + [ + [ + 0.025, + 0.025 + ] + ], + [ + [ + 0.075, + 0.025 + ] + ], + [ + [ + 0.125, + 0.025 + ] + ], + [ + [ + 0.175, + 0.025 + ] + ], + [ + [ + 0.225, + 0.025 + ] + ], + [ + [ + 0.275, + 0.025 + ] + ], + [ + [ + 0.025, + 0.075 + ] + ], + [ + [ + 0.075, + 0.075 + ] + ], + [ + [ + 0.125, + 0.075 + ] + ], + [ + [ + 0.175, + 0.075 + ] + ], + [ + [ + 0.225, + 0.075 + ] + ], + [ + [ + 0.275, + 0.075 + ] + ] +] \ No newline at end of file diff --git a/benchmarks/2d/uniaxial_particle_traction/test_benchmark.py b/benchmarks/2d/uniaxial_particle_traction/test_benchmark.py new file mode 100644 index 0000000..c63d167 --- /dev/null +++ b/benchmarks/2d/uniaxial_particle_traction/test_benchmark.py @@ -0,0 +1,29 @@ +import os +from pathlib import Path +import jax.numpy as jnp +from diffmpm.solver import MPM + +curr_filepath = Path(__file__).absolute() +curr_dir = curr_filepath.parent +os.chdir(curr_dir) +mpm = MPM("mpm-particle-traction.toml") +mpm.solve() +result = jnp.load("results/uniaxial-particle-traction/particles_0300.npz") +## Step 300 +assert jnp.round(result["stress"][0, :, 0].min() - 0.4450086768966724, 5) == 0.0 +assert jnp.round(result["stress"][0, :, 0].max() - 0.5966527842046769, 5) == 0.0 + +## Step 510 +result = jnp.load("results/uniaxial-particle-traction/particles_0510.npz") +assert jnp.round(result["stress"][0, :, 0].min() - 0.7528092313640623, 5) == 0.0 +assert jnp.round(result["stress"][0, :, 0].max() - 1.0109599915279937, 5) == 0.0 + +## Step 750 +result = jnp.load("results/uniaxial-particle-traction/particles_0750.npz") +assert jnp.round(result["stress"][0, :, 0].min() - 0.7500090055681591, 5) == 0.0 +assert jnp.round(result["stress"][0, :, 0].max() - 1.0000224746314728, 5) == 0.0 + +## Step 990 +result = jnp.load("results/uniaxial-particle-traction/particles_0990.npz") +assert jnp.round(result["stress"][0, :, 0].min() - 0.750002924022295, 5) == 0.0 +assert jnp.round(result["stress"][0, :, 0].max() - 0.9999997782938734, 5) == 0.0 diff --git a/diffmpm/element.py b/diffmpm/element.py index d2bce61..f0b9d4f 100644 --- a/diffmpm/element.py +++ b/diffmpm/element.py @@ -239,6 +239,17 @@ def apply_concentrated_nodal_forces(self, particles, curr_time): factor * cnf.force ) + def apply_particle_traction_forces(self, particles): + def _step(pid, args): + f_ext, ptraction, mapped_pos, el_nodes = args + f_ext = f_ext.at[el_nodes[pid]].add(mapped_pos[pid] @ ptraction[pid]) + return f_ext, ptraction, mapped_pos, el_nodes + + mapped_positions = self.shapefn(particles.reference_loc) + mapped_nodes = vmap(self.id_to_node_ids)(particles.element_ids).squeeze(-1) + args = (self.nodes.f_ext, particles.traction, mapped_positions, mapped_nodes) + self.nodes.f_ext, _, _, _ = lax.fori_loop(0, len(particles), _step, args) + def update_nodal_acceleration_velocity(self, particles, dt: float, *args): """Update the nodal momentum based on total force on nodes.""" total_force = self.nodes.get_total_force() diff --git a/diffmpm/forces.py b/diffmpm/forces.py index 59a31c9..eb6d27f 100644 --- a/diffmpm/forces.py +++ b/diffmpm/forces.py @@ -3,7 +3,7 @@ NodalForce = namedtuple("NodalForce", ("node_ids", "function", "dir", "force")) ParticleTraction = namedtuple( - "ParticleTraction", ("pset", "function", "dir", "traction") + "ParticleTraction", ("pset", "pids", "function", "dir", "traction") ) register_pytree_node( NodalForce, diff --git a/diffmpm/io.py b/diffmpm/io.py index 31e9cab..78f56b5 100644 --- a/diffmpm/io.py +++ b/diffmpm/io.py @@ -80,7 +80,7 @@ def _parse_external_loading(self, config): external_loading = {} external_loading["gravity"] = jnp.array(config["external_loading"]["gravity"]) external_loading["concentrated_nodal_forces"] = [] - external_loading["particle_surface_traction"] = [] + particle_surface_traction = [] if "concentrated_nodal_forces" in config["external_loading"]: cnf_list = [] for cnfconfig in config["external_loading"]["concentrated_nodal_forces"]: @@ -102,17 +102,23 @@ def _parse_external_loading(self, config): if "particle_surface_traction" in config["external_loading"]: pst_list = [] for pstconfig in config["external_loading"]["particle_surface_traction"]: - pst = ParticleTraction( - pset=jnp.array(pstconfig["pset"]), - function=self.parsed_config["math_functions"][ + if "math_function_id" in pstconfig: + fn = self.parsed_config["math_functions"][ pstconfig["math_function_id"] - ], + ] + else: + fn = Unit(-1) + pst = ParticleTraction( + pset=pstconfig["pset"], + pids=jnp.array(pstconfig["pids"]), + function=fn, dir=pstconfig["dir"], traction=pstconfig["traction"], ) pst_list.append(pst) - external_loading["particle_surface_traction"] = pst_list + particle_surface_traction.extend(pst_list) self.parsed_config["external_loading"] = external_loading + self.parsed_config["particle_surface_traction"] = particle_surface_traction def _parse_mesh(self, config): element_cls = getattr(mpel, config["mesh"]["element"]) diff --git a/diffmpm/mesh.py b/diffmpm/mesh.py index 9d967f1..1a31c60 100644 --- a/diffmpm/mesh.py +++ b/diffmpm/mesh.py @@ -21,6 +21,12 @@ def __init__(self, config: dict): """Initialize mesh using configuration.""" self.particles: Iterable[Particles, ...] = config["particles"] self.elements: _Element = config["elements"] + self.particle_tractions = config["particle_surface_traction"] + + @property + @abc.abstractmethod + def ndim(self): + ... # TODO: Convert to using jax directives for loop def apply_on_elements(self, function, args=()): @@ -34,15 +40,33 @@ def apply_on_particles(self, function, args=()): f = getattr(particle_set, function) f(self.elements, *args) + def apply_traction_on_particles(self, curr_time): + self.apply_on_particles("zero_traction") + for ptraction in self.particle_tractions: + factor = ptraction.function.value(curr_time) + traction_val = factor * ptraction.traction + for i, pset_id in enumerate(ptraction.pset): + self.particles[pset_id].assign_traction( + ptraction.pids[i], ptraction.dir, traction_val + ) + + # breakpoint() + self.apply_on_elements("apply_particle_traction_forces") + def tree_flatten(self): children = (self.particles, self.elements) - aux_data = tuple() + aux_data = self.particle_tractions return (children, aux_data) @classmethod def tree_unflatten(cls, aux_data, children): - del aux_data - return cls({"particles": children[0], "elements": children[1]}) + return cls( + { + "particles": children[0], + "elements": children[1], + "particle_surface_traction": aux_data, + } + ) @register_pytree_node_class @@ -61,6 +85,10 @@ def __init__(self, config: dict): """ super().__init__(config) + @property + def ndim(self): + return 1 + @register_pytree_node_class class Mesh2D(_MeshBase): @@ -78,6 +106,10 @@ def __init__(self, config: dict): """ super().__init__(config) + @property + def ndim(self): + return 2 + if __name__ == "__main__": from diffmpm.element import Linear1D diff --git a/diffmpm/particle.py b/diffmpm/particle.py index 119ff13..06396ff 100644 --- a/diffmpm/particle.py +++ b/diffmpm/particle.py @@ -55,6 +55,7 @@ def __init__( jnp.ones_like(self.mass) * self.material.properties["density"] ) self.volume = jnp.ones_like(self.mass) + self.size = jnp.zeros_like(self.loc) self.velocity = jnp.zeros_like(self.loc) self.acceleration = jnp.zeros_like(self.loc) self.momentum = jnp.zeros_like(self.loc) @@ -63,6 +64,7 @@ def __init__( self.strain_rate = jnp.zeros((self.loc.shape[0], 6, 1)) self.dstrain = jnp.zeros((self.loc.shape[0], 6, 1)) self.f_ext = jnp.zeros_like(self.loc) + self.traction = jnp.zeros_like(self.loc) self.reference_loc = jnp.zeros_like(self.loc) self.volumetric_strain_centroid = jnp.zeros((self.loc.shape[0], 1)) else: @@ -70,6 +72,7 @@ def __init__( self.mass, self.density, self.volume, + self.size, self.velocity, self.acceleration, self.momentum, @@ -78,6 +81,7 @@ def __init__( self.strain_rate, self.dstrain, self.f_ext, + self.traction, self.reference_loc, self.volumetric_strain_centroid, ) = data @@ -91,6 +95,7 @@ def tree_flatten(self): self.mass, self.density, self.volume, + self.size, self.velocity, self.acceleration, self.momentum, @@ -99,6 +104,7 @@ def tree_flatten(self): self.strain_rate, self.dstrain, self.f_ext, + self.traction, self.reference_loc, self.volumetric_strain_centroid, ) @@ -152,6 +158,7 @@ def compute_volume(self, elements, total_elements): / particles_per_element[self.element_ids] ) self.volume = self.volume.at[:, 0, 0].set(vol) + self.size = self.size.at[:].set(self.volume ** (1 / self.size.shape[-1])) self.mass = self.mass.at[:, 0, 0].set(vol * self.density.squeeze()) def update_natural_coords(self, elements: _Element): @@ -304,6 +311,14 @@ def update_volume(self, *args): self.volume = self.volume.at[:, 0, :].multiply(1 + self.dvolumetric_strain) self.density = self.density.at[:, 0, :].divide(1 + self.dvolumetric_strain) + def assign_traction(self, pids, dir, traction_): + self.traction = self.traction.at[pids, 0, dir].add( + traction_ * self.volume[pids, 0, 0] / self.size[pids, 0, dir] + ) + + def zero_traction(self, *args): + self.traction = self.traction.at[:].set(0) + if __name__ == "__main__": from diffmpm.material import SimpleMaterial diff --git a/diffmpm/scheme.py b/diffmpm/scheme.py index c5de119..83e35ca 100644 --- a/diffmpm/scheme.py +++ b/diffmpm/scheme.py @@ -25,6 +25,7 @@ def compute_stress_strain(self): def compute_forces(self, gravity, step): self.mesh.apply_on_elements("compute_external_force") self.mesh.apply_on_elements("compute_body_force", args=(gravity,)) + self.mesh.apply_traction_on_particles(step * self.dt) self.mesh.apply_on_elements( "apply_concentrated_nodal_forces", args=(step * self.dt,) ) diff --git a/examples/input_2d.toml b/examples/input_2d.toml index 06d8570..d44abdb 100644 --- a/examples/input_2d.toml +++ b/examples/input_2d.toml @@ -20,8 +20,8 @@ nsteps = 1000 velocity_update = true [output] -type = "hdf5" -file = "results/example_2d_out.hdf5" +format = "none" +folder = "results/" step_frequency = 5 [mesh] @@ -30,8 +30,8 @@ step_frequency = 5 # boundary_nodes = "boundary-1d.txt" # particle_element_ids = "particles-elements.txt" type = "generator" -nelements = 1 -element_length = 1.0 +nelements = [3, 1] +element_length = [0.1, 0.1] particle_element_ids = [0] element = "Quadrilateral4Node" @@ -49,11 +49,11 @@ velocity = 0.1 id = 0 density = 1 poisson_ratio = 1 -E = 100 +youngs_modulus = 100 type = "LinearElastic" [[particles]] -file = "examples/particles-2d.json" +file = "examples/particles-2d-nodal-force.json" material_id = 0 init_velocity = 0.1 @@ -67,12 +67,13 @@ dir = 1 force = 10.5 [[external_loading.particle_surface_traction]] -pset = [1] +pset = [0] +pids = [1, 3] dir = 1 math_function_id = 0 -traction = 10.5 +traction = 10 [[math_functions]] type = "Linear" xvalues = [0.0, 0.5, 1.0] -fxvalues = [0.0, 1.0, 1.0] +fxvalues = [1.0, 1.0, 1.0] diff --git a/tests/test_element.py b/tests/test_element.py index 02d095d..2b72460 100644 --- a/tests/test_element.py +++ b/tests/test_element.py @@ -202,3 +202,11 @@ def test_compute_internal_force(self, elements, particles): def test_compute_volume(self, elements): elements.compute_volume() assert jnp.allclose(elements.volume, jnp.array([1]).reshape(1, 1, 1)) + + def test_apply_particle_traction_forces(self, elements, particles): + particles.traction += jnp.array([1, 0]) + elements.apply_particle_traction_forces(particles) + assert jnp.allclose( + elements.nodes.f_ext, + jnp.array([[0.5, 0], [0.5, 0], [0.5, 0], [0.5, 0]]).reshape(4, 1, 2), + ) diff --git a/tests/test_particle.py b/tests/test_particle.py index 6b8d675..9dd9a5e 100644 --- a/tests/test_particle.py +++ b/tests/test_particle.py @@ -47,3 +47,15 @@ def test_compute_strain(self, elements, particles): def test_compute_volume(self, elements, particles): particles.compute_volume(elements, elements.total_elements) assert jnp.allclose(particles.volume, jnp.array([0.5, 0.5]).reshape(2, 1, 1)) + + def test_assign_traction(self, elements, particles): + particles.compute_volume(elements, elements.total_elements) + particles.assign_traction(jnp.array([0]), 1, 10) + assert jnp.allclose( + particles.traction, jnp.array([[0, 7.071068], [0, 0]]).reshape(2, 1, 2) + ) + + def test_zero_traction(self, particles): + particles.traction += 1 + particles.zero_traction() + assert jnp.all(particles.traction == 0)