Examples
Make 3D topography model from mesh files
This example generates material point method model from external mesh.obj files which defines layer boundaries of
ground of mountain.
import sys
sys.path.append('/work2/08264/baagee/frontera/cbgeopy')
import utils
from mpm import MPMConfig
import demo_utils
from functools import partial
import numpy as np
import vis_utils
import os
import trimesh
save_dir = './'
# Replace the data_dir to the directory containing the obj surface meshes
data_dir = '/work2/08264/baagee/frontera/fundao/fundao_3d-8/'
# Create surface
mesh_base = utils.obj2mesh(f'{data_dir}/mesh_0_v2.obj') # base
mesh_dike = utils.obj2mesh(f'{data_dir}/mesh_1_v2.obj') # dike
mesh_sand1 = utils.obj2mesh(f'{data_dir}/mesh_2_v2.obj') # sand (or slime)
mesh_slime1 = utils.obj2mesh(f'{data_dir}/mesh_2_slime.obj') # slime-1
mesh_sand_right = utils.obj2mesh(f'{data_dir}/mesh_2_right_sand.obj') # sand right
mesh_sand2 = utils.obj2mesh(f'{data_dir}/mesh_7.obj') # sand
mesh_sand3 = utils.obj2mesh(f'{data_dir}/mesh_7_sand.obj') # sand below slime-2
mesh_slime2 = utils.obj2mesh(f'{data_dir}/mesh_7_slime.obj') # slime-2
mesh_top = utils.obj2mesh(f'{data_dir}/mesh_8.obj') # sand
mesh_top_ref = utils.obj2mesh(f'{data_dir}/mesh_top.obj') # reference top surface for stress initialization
# Make it to trimesh object
floor = demo_utils.gen_surface_mesh(
[200, 1130], [-850, 450], 10, partial(demo_utils.generate_slope, slope_angle=0, z_min=740, amplitude=0),
# [180, 1150], [-300, 450], 10, partial(demo_utils.generate_slope, slope_angle=0, z_min=740, amplitude=0),
plot=False)
mesh_base = trimesh.Trimesh(*mesh_base)
mesh_dike = trimesh.Trimesh(*mesh_dike)
mesh_sand1 = trimesh.Trimesh(*mesh_sand1)
mesh_slime1 = trimesh.Trimesh(*mesh_slime1)
mesh_sand_right = trimesh.Trimesh(*mesh_sand_right)
mesh_sand2 = trimesh.Trimesh(*mesh_sand2)
mesh_sand3 = trimesh.Trimesh(*mesh_sand3)
mesh_slime2 = trimesh.Trimesh(*mesh_slime2)
mesh_top = trimesh.Trimesh(*mesh_top)
mesh_top_ref = trimesh.Trimesh(*mesh_top_ref)
# # See surfaces
# vis_utils.plot_surfaces(floor, mesh_base, mesh_sand_up, mesh_top, save_path='./surface.html')
# Set config
lx, ly, lz = 1130 - 200, 450 - (-850), 950 - 740
mpm = MPMConfig(domain_origin=[200, -850, 740], domain_length=[lx, ly, lz])
# lx, ly, lz = 1150 - 180, 450 - (-300), 950 - 780
# mpm = MPMConfig(domain_origin=[180, -300, 780], domain_length=[lx, ly, lz])
# Mesh
cell_size = 10
mpm.add_mesh(
n_cells_per_dim=[int(lx/cell_size), int(ly/cell_size), int(lz/cell_size)])
# Add materials
mpm.add_materials(
[
{
"id": 0,
"name": "liq_sand",
"density": 2200,
"youngs_modulus": 4e6,
"poisson_ratio": 0.30,
"friction": 0.0,
"dilation": 0.0,
"cohesion": 10e3,
"tension_cutoff": 10,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 45.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb3D"
},
{
"id": 1,
"density": 2200,
"name": "slime",
"youngs_modulus": 2e6,
"poisson_ratio": 0.30,
"friction": 0.0,
"dilation": 0.0,
"cohesion": 50e3,
"tension_cutoff": 10,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 45.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb3D"
},
{
"id": 2,
"name": "bedrock",
"type": "LinearElastic3D",
"density": 2600,
"youngs_modulus": 8e6,
"poisson_ratio": 0.3
},
{
"id": 3,
"name": "sand",
"density": 2200,
"youngs_modulus": 4e6,
"poisson_ratio": 0.30,
"friction": 30.0,
"dilation": 0.0,
"cohesion": 1000,
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 45.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb3D"
},
{
"id": 10,
"name": "geostatic_le",
"type": "LinearElastic3D",
"density": 2200,
"youngs_modulus": 8e6,
"poisson_ratio": 0.3
}
]
)
n_particle_per_cell = 2
z_find_method = 'linear'
z_fill_method = 'round'
overlap_tolerance = 1
mpm.add_particles_from_topography(
lower_topography=floor,
upper_topography=mesh_base,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=0,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_base,
upper_topography=mesh_dike,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=1,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_dike,
upper_topography=mesh_sand1,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=2,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_sand1,
upper_topography=mesh_slime1,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=3,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_slime1,
upper_topography=mesh_sand_right,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=4,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_sand_right,
upper_topography=mesh_sand2,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=5,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_sand2,
upper_topography=mesh_sand3,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=6,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_sand3,
upper_topography=mesh_slime2,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=7,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
mpm.add_particles_from_topography(
lower_topography=mesh_slime2,
upper_topography=mesh_top,
n_particle_per_cell=n_particle_per_cell,
material_id=10,
particle_group_id=8,
z_find_method=z_find_method,
base_find_method='simple',
z_fill_method=z_fill_method,
overlap_tolerance=overlap_tolerance
)
# mpm.add_particles_from_topography(
# lower_topography=mesh_top,
# upper_topography=mesh_top2,
# n_particle_per_cell=n_particle_per_cell,
# material_id=10,
# particle_group_id=8,
# z_find_method=z_find_method,
# base_find_method='simple',
# z_fill_method='round'
# )
mpm.define_particle_entity()
# Stress initialization
mpm.add_initial_stress(
option='k0',
top_surface=mesh_top_ref,
density=2200,
k0=0.5
)
# Boundary constraints
mpm.define_boundary_entity()
mpm.add_velocity_constraints(
[
{"axis": "x", "bound_loc": "start", "velocity": 0.0},
{"axis": "x", "bound_loc": "end", "velocity": 0.0},
{"axis": "y", "bound_loc": "start", "velocity": 0.0},
{"axis": "y", "bound_loc": "end", "velocity": 0.0},
{"axis": "z", "bound_loc": "start", "velocity": 0.0},
{"axis": "z", "bound_loc": "end", "velocity": 0.0}
]
)
mpm.add_friction_constrains(
[
{"axis": "x", "bound_loc": "start", "sign_n": -1, "friction": 0.3},
{"axis": "x", "bound_loc": "end", "sign_n": 1, "friction": 0.3},
{"axis": "y", "bound_loc": "start", "sign_n": -1, "friction": 0.3},
{"axis": "y", "bound_loc": "end", "sign_n": 1, "friction": 0.3},
{"axis": "z", "bound_loc": "start", "sign_n": -1, "friction": 0.3},
{"axis": "z", "bound_loc": "end", "sign_n": 1, "friction": 0.3}
]
)
# External loading conditions
mpm.add_external_loadings(
{"gravity": [0, 0, -9.81]}
)
# Analysis settings
mpm.analysis({
"mpm_scheme": "usl",
"locate_particles": False,
"dt": 1e-04,
"damping": {
"type": "Cundall",
"damping_factor": 0.05
},
"resume": {
"resume": False,
"step": 0,
"uuid": "sand3d"
},
"velocity_update": False,
"nsteps": int(1e6),
"type": "MPMExplicit3D",
"uuid": "sand3d"
})
# Post-processing
mpm.post_processing({
"path": "results/",
"output_steps": 2000,
"vtk": [
"displacements",
"stresses"
]
})
mpm.write(save_dir=save_dir)
# mpm.visualize_mesh(save_path=f'{save_dir}/mesh_config.html', node_indices=True)
# mpm.visualize_particles(save_path=f'{save_dir}/particle_config.html')
vis_utils.plot_cross_section(mpm.particle_groups, 'yz', 800, tolerance=3, grid_spacing=cell_size)
# vis_utils.plot_surfaces(floor, mesh_base, mesh_sand_up, mesh_top,
# points=mpm.particle_groups,
# resolution=1,
# save_path=f'{save_dir}/surfaces.html')
vis_utils.save_as_vtk(
meshes=(floor, mesh_base, mesh_dike, mesh_sand1, mesh_slime1, mesh_sand2, mesh_sand3, mesh_slime2, mesh_top),
points=mpm.particle_groups
)
# Save the current script
# Get the path of the currently running script (main.py)
current_script_path = os.path.abspath(__file__)
utils.save_script(
current_script_path,
save_path=f'{save_dir}/input_script.py')
Make MP model based on 2D lines
This example generates material point method model from user-defined points defining the boundary lines that distinguishes the layers.
import random
import utils
from mpm import MPMConfig
import vis_utils
import os
from tools import random_generator
save_dir = '../layers2d_from_lines/'
# Set config
lx, ly = 1250.0, 150.0
origin_x, origin_y = 0, 0
mpm = MPMConfig(domain_origin=[origin_x, origin_y], domain_length=[lx, ly])
# Mesh
cell_size = 4
mpm.add_mesh(
n_cells_per_dim=[round(lx/cell_size), round(ly/cell_size)])
# Materials
mpm.add_materials([
{
"id": 0,
"name": "bedrock",
"type" : "LinearElastic2D",
"density" : 1800,
"youngs_modulus" : 50000000,
"poisson_ratio" : 0.3
},
{
"id": 1,
"name": "slime",
"density": 1800,
"youngs_modulus": 20000000,
"poisson_ratio": 0.3,
"friction": 7,
"dilation": 0.0,
"cohesion": 1000,
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.005,
"residual_friction": 18.0,
"residual_dilation": 0.0,
"residual_cohesion": 50.0,
"residual_pdstrain": 0.05,
"type": "MohrCoulomb2D"
},
{
"id": 2,
"name": "sand",
"density": 1800,
"youngs_modulus": 20000000,
"poisson_ratio": 0.3,
"friction": 10.0,
"dilation": 0.0,
"cohesion": 1000,
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.005,
"residual_friction": 18.0,
"residual_dilation": 0.0,
"residual_cohesion": 50.0,
"residual_pdstrain": 0.05,
"type": "MohrCoulomb2D"
},
{
"id": 3,
"name": "comp_sand",
"density": 2200,
"youngs_modulus": 20000000,
"poisson_ratio": 0.3,
"friction": 13.0,
"dilation": 0.0,
"cohesion": 1000,
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.005,
"residual_friction": 18.0,
"residual_dilation": 0.0,
"residual_cohesion": 50.0,
"residual_pdstrain": 0.05,
"type": "MohrCoulomb2D"
},
{
"id": 100,
"type": "LinearElastic2D",
"density": 1800,
"youngs_modulus": 5e7,
"poisson_ratio": 0.3
}
]
)
# Particle
points = [
[0, 98],
[37, 97],
[97, 53],
[151, 47],
[423, 40],
[423, 53],
[468, 55],
[491, 40],
[746, 36],
[1250, 12],
[62, 78],
[410, 78],
[0, 120],
[290, 123],
[354, 102],
[541, 78],
[574, 78],
[655, 56],
[673, 37],
[331, 123],
[380, 104],
[486, 90],
[573, 89]
]
mpm.add_particles_from_lines(
layer_info=[
{
"line_points": [points[i] for i in [0, 1, 10, 3, 4, 5, 6, 7, 18, 8, 9]],
"material_id": 0,
"particle_group_id": 0
},
{
"line_points": [points[i] for i in [0, 1, 10, 11, 5, 6, 7, 18, 8, 9]],
"material_id": 1,
"particle_group_id": 1
},
{
"line_points": [points[i] for i in [12, 13, 14, 15, 16, 17, 18, 8, 9]],
"material_id": 2,
"particle_group_id": 2
},
{
"line_points": [points[i] for i in [12, 13, 19, 20, 21, 22, 8, 9]],
"material_id": 3,
"particle_group_id": 3
}
],
n_particle_per_cell=2
)
mpm.remove_overlapping_particles(overlap_tolerance=0.001)
mpm.define_particle_entity()
# Boundary constraints
mpm.define_boundary_entity()
mpm.add_velocity_constraints(
[
{"axis": "x", "bound_loc": "start", "velocity": 0.0},
{"axis": "x", "bound_loc": "end", "velocity": 0.0},
{"axis": "y", "bound_loc": "start", "velocity": 0.0},
{"axis": "y", "bound_loc": "end", "velocity": 0.0},
{"axis": "z", "bound_loc": "start", "velocity": 0.0},
{"axis": "z", "bound_loc": "end", "velocity": 0.0}
]
)
mpm.add_friction_constrains(
[
{"axis": "x", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "x", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "y", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "y", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "z", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "z", "bound_loc": "end", "sign_n": 1, "friction": 0.38}
]
)
# External loading conditions
mpm.add_external_loadings(
{"gravity": [0, -9.81]}
)
# Analysis settings
mpm.analysis({
"mpm_scheme": "usf",
"locate_particles": False,
"dt": 1e-04,
"damping": {
"type": "Cundall",
"damping_factor": 0.05
},
"resume": {
"resume": False,
"step": 0,
"uuid": "sand2d"
},
"velocity_update": False,
"nsteps": int(1.8e5),
"type": "MPMExplicit2D",
"uuid": "sand2d"
})
# Post-processing
mpm.post_processing({
"path": "results/",
"output_steps": 375,
"vtk": [
"displacements"
]
})
mpm.write(save_dir=save_dir)
vis_utils.plot_scatter(mpm.particle_groups, mpm.domain_ranges, f'{save_dir}/particle_config.png')
# Save the current script
# Get the path of the currently running script (main.py)
current_script_path = os.path.abspath(__file__)
utils.save_script(
current_script_path,
save_path=f'{save_dir}/input_script.py')
Make MP model based on polygon
This example generates material point method model from user-defined polygons.
import sys
sys.path.append('/work/08264/baagee/frontera/cbgeopy/')
import random
import utils
from mpm import MPMConfig
import vis_utils
import os
from tools import random_generator
from numpy.random import uniform as randf
from random import randint
# Random parameters
sim_id_range = range(900, 901)
for i in sim_id_range:
print(f"Generate mpm inputs for simulation {i}...")
save_dir = f'./sim-{i}'
# Set config
lx, ly = 1000.0, 152.0
origin_x, origin_y = 0, 0
mpm = MPMConfig(domain_origin=[origin_x, origin_y], domain_length=[lx, ly])
# Mesh
cell_size = 4
mpm.add_mesh(
n_cells_per_dim=[round(lx/cell_size), round(ly/cell_size)])
# Materials
cohesion_options = [10e3, 20e3, 30e3, 40e3, 50e3, 70e3, 100e3, 130e3]
friction_options = [0, 0, 0, 0, 10, 17.5, 25, 32.5, 40, 43]
mpm.add_materials([
{
"id": 0,
"type": "LinearElastic2D",
"density": 1800,
"youngs_modulus": 50000000.0,
"poisson_ratio": 0.3
},
{
"id": 1,
"density": 1800,
"youngs_modulus": 40e6,
"poisson_ratio": 0.3,
"friction": random.choice(friction_options),
"dilation": 0.0,
"cohesion": random.choice(cohesion_options),
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 30.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb2D"
},
{
"id": 2,
"density": 1800,
"youngs_modulus": 40e6,
"poisson_ratio": 0.3,
"friction": random.choice(friction_options),
"dilation": 0.0,
"cohesion": random.choice(cohesion_options),
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 30.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb2D"
},
{
"id": 3,
"density": 1800,
"youngs_modulus": 40e6,
"poisson_ratio": 0.3,
"friction": random.choice(friction_options),
"dilation": 0.0,
"cohesion": random.choice(cohesion_options),
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 30.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb2D"
},
{
"id": 4,
"density": 1800,
"youngs_modulus": 40e6,
"poisson_ratio": 0.3,
"friction": random.choice(friction_options),
"dilation": 0.0,
"cohesion": random.choice(cohesion_options),
"tension_cutoff": 100,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 30.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb2D"
}
])
# Define geometry points
p1 = [randint(250, 350), 4]
p2 = [randint(360, 640), 4]
p3 = [randint(650, 750), 4]
p5 = [randint(464, 486), randint(12, 70)]
p7 = [randint(410, 460), randint(74, 85)]
p8 = [randint(490, 540), randint(74, 95)]
# Linear interpolation between p7 and p1
t1 = random.uniform(0.2, 0.8)
t2 = random.uniform(0.2, 0.8)
t3 = random.uniform(0.3, 0.7)
p4 = [p7[i] + t1 * (p1[i] - p7[i]) for i in range(len(p7))]
p6 = [p8[i] + t2 * (p3[i] - p8[i]) for i in range(len(p8))]
p78 = [p7[i] + t3 * (p8[i] - p7[i]) for i in range(len(p7))]
# Particle
mpm.add_particles_from_lines(
layer_info=[
{
"line_points": [[0, 4], [1000, 4]],
"material_id": 0,
"particle_group_id": 0,
"randomness": 0.0,
}
],
n_particle_per_cell=2
)
mpm.add_particles_from_polygon(
polygon_info=[
{
"polygon_points": [p1, p2, p5, p4],
"material_id": 1,
"particle_group_id": 1
},
{
"polygon_points": [p2, p3, p6, p5],
"material_id": 2,
"particle_group_id": 2
},
{
"polygon_points": [p4, p5, p78, p7],
"material_id": 3,
"particle_group_id": 3
},
{
"polygon_points": [p5, p6, p8, p78],
"material_id": 4,
"particle_group_id": 4
}
],
n_particle_per_cell=2,
randomness=0.5
)
mpm.remove_overlapping_particles(overlap_tolerance=0.001)
mpm.define_particle_entity()
# Boundary constraints
mpm.define_boundary_entity()
mpm.add_velocity_constraints(
[
{"axis": "x", "bound_loc": "start", "velocity": 0.0},
{"axis": "x", "bound_loc": "end", "velocity": 0.0},
{"axis": "y", "bound_loc": "start", "velocity": 0.0},
{"axis": "y", "bound_loc": "end", "velocity": 0.0},
{"axis": "z", "bound_loc": "start", "velocity": 0.0},
{"axis": "z", "bound_loc": "end", "velocity": 0.0}
]
)
mpm.add_friction_constrains(
[
{"axis": "x", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "x", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "y", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "y", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "z", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "z", "bound_loc": "end", "sign_n": 1, "friction": 0.38}
]
)
# External loading conditions
mpm.add_external_loadings(
{"gravity": [0, -9.81]}
)
# Analysis settings
mpm.analysis({
"mpm_scheme": "usf",
"locate_particles": False,
"dt": 1e-04,
"damping": {
"type": "Cundall",
"damping_factor": 0.05
},
"resume": {
"resume": False,
"step": 0,
"uuid": "sand2d"
},
"velocity_update": False,
"nsteps": 225000,
"type": "MPMExplicit2D",
"uuid": "sand2d"
})
# Post-processing
mpm.post_processing({
"path": "results/",
"output_steps": 375,
"vtk": [
"displacements",
"stresses"
]
})
# Save mpm json that will be used after stress initialization
mpm.mpm_json["mesh"]["particles_stresses"] = "particles-stresses.txt"
mpm.write(save_dir=save_dir, file_name='mpm-resume.json')
# mpm json for stress initialization
mpm.mpm_json["mesh"].pop("particles_stresses")
mpm.mpm_json["analysis"]["uuid"] = "sand2d-le"
for particle in mpm.mpm_json["particles"]:
# Set all materials to bedrock
particle["generator"]["material_id"] = 0
mpm.mpm_json["analysis"]["resume"]["resume"] = False
mpm.mpm_json["analysis"]["resume"]["uuid"] = "sand2d-le"
mpm.mpm_json["analysis"]["nsteps"] = 70001
# Overwrite
mpm.write(save_dir=save_dir, file_name='mpm-le.json')
vis_utils.plot_scatter(mpm.particle_groups, mpm.domain_ranges, f'{save_dir}/particle_config.png')
# Save the current script
# Get the path of the currently running script (main.py)
current_script_path = os.path.abspath(__file__)
utils.save_script(
current_script_path,
save_path=f'{save_dir}/input_script.py')
Make MP model for sand cube collision
This example generates material point method model from two colliding cubes. This saves
import trimesh
import utils
from mpm import MPMConfig
import demo_utils
from functools import partial
import numpy as np
import vis_utils
import os
import trimesh
save_dir = './'
# Set config
lx, ly = 200.0, 200.0
mpm = MPMConfig(domain_origin=[0, 0], domain_length=[lx, ly])
# Mesh
cell_size = 5
mpm.add_mesh(
n_cells_per_dim=[int(lx/cell_size), int(ly/cell_size)])
# Add materials
mpm.add_materials(
[
{
"id": 0,
"name": "soil-0",
"density": 1800,
"youngs_modulus": 20000000.0,
"poisson_ratio": 0.3,
"friction": 40.0,
"dilation": 0.0,
"cohesion": 100,
"tension_cutoff": 50,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 30.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb2D"
},
{
"id": 1,
"name": "soil-2",
"density": 1800,
"youngs_modulus": 20000000.0,
"poisson_ratio": 0.3,
"friction": 15.0,
"dilation": 0.0,
"cohesion": 100,
"tension_cutoff": 50,
"softening": False,
"peak_pdstrain": 0.0,
"residual_friction": 30.0,
"residual_dilation": 0.0,
"residual_cohesion": 0.0,
"residual_pdstrain": 0.0,
"type": "MohrCoulomb2D"
},
{
"id": 10,
"type": "LinearElastic2D",
"density": 1800,
"youngs_modulus": 50000000.0,
"poisson_ratio": 0.3
}
]
)
# Particles
mpm.add_particles_from_lines(
layer_info=[
{
"line_points": [[0, 50], [80, 10], [200, 30]],
"material_id": 10,
"particle_group_id": 0
}
],
n_particle_per_cell=2
)
mpm.add_particles_cube(
cube_origin=[40, 40],
cube_length=[60, 60],
material_id=0,
n_particle_per_cell=2,
particle_group_id=1
)
mpm.add_particles_cube(
cube_origin=[130, 80],
cube_length=[60, 60],
material_id=1,
n_particle_per_cell=2,
particle_group_id=2
)
mpm.remove_overlapping_particles(overlap_tolerance=0.001)
mpm.define_particle_entity()
vis_utils.plot_scatter(mpm.particle_groups, mpm.domain_ranges, 'particle_config.png')
# Boundary constraints
mpm.define_boundary_entity()
mpm.add_velocity_constraints(
[
{"axis": "x", "bound_loc": "start", "velocity": 0.0},
{"axis": "x", "bound_loc": "end", "velocity": 0.0},
{"axis": "y", "bound_loc": "start", "velocity": 0.0},
{"axis": "y", "bound_loc": "end", "velocity": 0.0},
{"axis": "z", "bound_loc": "start", "velocity": 0.0},
{"axis": "z", "bound_loc": "end", "velocity": 0.0}
]
)
mpm.add_friction_constrains(
[
{"axis": "x", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "x", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "y", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "y", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "z", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "z", "bound_loc": "end", "sign_n": 1, "friction": 0.38}
]
)
mpm.add_particle_constraints(
[
{
"pset_id": 0,
"axis": 'x',
"velocity": 0.0
},
{
"pset_id": 0,
"axis": 'y',
"velocity": 0.0
},
{
"pset_id": 1,
"axis": 'x',
"velocity": 10.0
},
{
"pset_id": 1,
"axis": 'y',
"velocity": 0.0
},
{
"pset_id": 2,
"axis": 'x',
"velocity": -10.0
},
{
"pset_id": 2,
"axis": 'y',
"velocity": 0.0
}
]
)
# External loading conditions
mpm.add_external_loadings(
{"gravity": [0, -9.81]}
)
# Analysis settings
mpm.analysis({
"mpm_scheme": "usf",
"locate_particles": False,
"dt": 1e-04,
"damping": {
"type": "Cundall",
"damping_factor": 0.05
},
"resume": {
"resume": False,
"step": 0,
"uuid": "sand2d"
},
"velocity_update": False,
"nsteps": int(1.5e5),
"type": "MPMExplicit2D",
"uuid": "sand2d"
})
# Post-processing
mpm.post_processing({
"path": "results/",
"output_steps": 375,
"vtk": [
"displacements",
"stresses"
]
})
mpm.write(save_dir=save_dir)
# # Update input json for resuming without particle velocity constraints
# current_particle_constraints = mpm.mpm_json['mesh']['boundary_conditions']['particles_velocity_constraints']
# # Only fix the velocities for the bedrock
# new_constraints = [constraint for constraint in current_particle_constraints if constraint['pset_id'] == 0]
# mpm.mpm_json['mesh']['boundary_conditions']['particles_velocity_constraints'] = new_constraints
mpm.mpm_json['mesh']['boundary_conditions'].pop('particles_velocity_constraints')
mpm.mpm_json['analysis']['resume']['resume'] = True
# Overwrite
mpm.write(save_dir=save_dir, file_name='mpm-resume.json')
vis_utils.save_as_vtk(
meshes=None,
points=mpm.particle_groups
)
# Save the current script
# Get the path of the currently running script (main.py)
current_script_path = os.path.abspath(__file__)
utils.save_script(
current_script_path,
save_path=f'{save_dir}/input_script.py')
Make MP model for sand random field
This example generates material point method model from user-defined polygons and random field parameters.
import sys
import trimesh
import utils
from mpm import MPMConfig
import demo_utils
from functools import partial
import numpy as np
import vis_utils
import os
import trimesh
save_dir = './sand_random_field'
# Set config
lx, ly = 200.0, 60.0
mpm = MPMConfig(domain_origin=[0, 0], domain_length=[lx, ly])
# Mesh
cell_size = 4
mpm.add_mesh(
n_cells_per_dim=[int(lx/cell_size), int(ly/cell_size)])
# Add materials
bedrock = [
{
"id": 0,
"type": "LinearElastic2D",
"density": 1800,
"youngs_modulus": 50000000.0,
"poisson_ratio": 0.3
}
]
# Particles
# Bedrock
mpm.add_particles_from_lines(
layer_info=[
{
"line_points": [[0, 4], [300, 4]],
"material_id": 0,
"particle_group_id": 0
}
],
n_particle_per_cell=2
)
# Slope
mpm.add_particles_from_random_field(
polygons_params=[
{
"polygon_points": [[0, 4], [60, 4], [60, 20], [0, 20]],
"random_params": {"mean": 12.0, "std": 2.0, "len_scale": 10.0},
},
{
"polygon_points": [[0, 20], [60, 20], [0, 40]],
"random_params": {"mean": 20.0, "std": 2.0, "len_scale": 2.0}
}
],
n_particle_per_cell=2
)
mpm.remove_overlapping_particles(overlap_tolerance=0.001)
mpm.define_particle_entity()
# Material
mpm.add_materials(materials=bedrock)
mpm.add_materials(option="random_field", material_type="MohrCoulomb2D")
# Boundary constraints
mpm.define_boundary_entity()
mpm.add_velocity_constraints(
[
{"axis": "x", "bound_loc": "start", "velocity": 0.0},
{"axis": "x", "bound_loc": "end", "velocity": 0.0},
{"axis": "y", "bound_loc": "start", "velocity": 0.0},
{"axis": "y", "bound_loc": "end", "velocity": 0.0},
{"axis": "z", "bound_loc": "start", "velocity": 0.0},
{"axis": "z", "bound_loc": "end", "velocity": 0.0}
]
)
mpm.add_friction_constrains(
[
{"axis": "x", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "x", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "y", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "y", "bound_loc": "end", "sign_n": 1, "friction": 0.38},
{"axis": "z", "bound_loc": "start", "sign_n": -1, "friction": 0.38},
{"axis": "z", "bound_loc": "end", "sign_n": 1, "friction": 0.38}
]
)
# External loading conditions
mpm.add_external_loadings(
{"gravity": [0, -9.81]}
)
# Analysis settings
mpm.analysis({
"mpm_scheme": "usf",
"locate_particles": False,
"dt": 1e-04,
"damping": {
"type": "Cundall",
"damping_factor": 0.05
},
"resume": {
"resume": False,
"step": 0,
"uuid": "sand2d"
},
"velocity_update": False,
"nsteps": 225000,
"type": "MPMExplicit2D",
"uuid": "sand2d"
})
# Post-processing
mpm.post_processing({
"path": "results/",
"output_steps": 375,
"vtk": [
"displacements",
"stresses"
]
})
# Save mpm json that will be used after stress initialization
mpm.mpm_json["mesh"]["particles_stresses"] = "particles-stresses.txt"
mpm.write(save_dir=save_dir, file_name='mpm-resume.json')
# mpm json for stress initialization
mpm.mpm_json["mesh"].pop("particles_stresses")
mpm.mpm_json["analysis"]["uuid"] = "sand2d-le"
for particle in mpm.mpm_json["particles"]:
# Set all materials to bedrock
particle["generator"]["material_id"] = 0
mpm.mpm_json["analysis"]["resume"]["resume"] = False
mpm.mpm_json["analysis"]["resume"]["uuid"] = "sand2d-le"
mpm.mpm_json["analysis"]["nsteps"] = 70001
# Overwrite
mpm.write(save_dir=save_dir, file_name='mpm-le.json')
# vis_utils.plot_scatter(mpm.particle_groups, mpm.domain_ranges, f'{save_dir}/particle_config.png')
vis_utils.plot_random_field_data(
mpm.particle_groups, mpm.domain_ranges, f'{save_dir}/random_field_data.png')
# Save the current script
# Get the path of the currently running script (main.py)
current_script_path = os.path.abspath(__file__)
utils.save_script(
current_script_path,
save_path=f'{save_dir}/input_script.py')