Material point method¶

Theory¶

Development of the material point method in Ratel is in progress.

Initializing Material Properties¶

There are three ways to assign material properties to material points.

Conforming Mesh¶

The simplest approach simply uses the mesh regions of the background mesh in the initial configuration to assign material properties to points. As for FEM, define a material for each mesh region. Then, before the first timestep, each material point will have its local material properties set according to the mesh cell in which it is contained.

This initialization method is the easiest to use, but requires that a conformal hexahedral mesh of the desired geometry can be defined. Additionally, since the background mesh is used for both initialization and subsequent solves, any poorly shaped elements in the mesh due to the conforming geometry will negatively impact the numerical performance and accuracy of the solver.

Command-Line Options¶

For completeness, MPM Initialization Options shows the command-line options for initializing material points in a conformal mesh. In addition to the MPM specific options, materials and mesh options should be specified in the same fashion as for FEM, see Command Line Options.

Table 15 MPM Initialization Options¶
Option	Description	Default Value
`-method mpm`	Required to enable MPM
`-mpm_point_location_type [uniform \| gauss \| random]`	Initial per-cell point distribution within each background mesh cell.	`uniform`
`-mpm_num_points_per_cell [int]`	Number of material points to place in each cell. For `uniform` or `gauss` point location type, this will be rounded to the nearest cube. Required.

Example¶

Material Mesh¶

Often, one may have a mesh with well-defined material regions for use with FEM which is not suitable for use as an MPM background mesh. For example, a tetrahedral mesh or a mesh with highly distorted elements. In this case, material properties can be assigned to material points based on a secondary material mesh. The material mesh can be defined using any format supported by the PETSc unstructured mesh DMPlex interface. This includes:

Simplex/tetrahedral meshes
Mixed topology meshes, such as hexahedral, pyramid, and simplex hybrid meshes
High-order geometry
Arbitrary number of volumetric mesh regions

Warning

Performance Implications of Using a Material Mesh

The material mesh initialization interface is very flexible, but can be extremely slow in certain circumstances. The implementation reads the material mesh, sets it as the background mesh for the DMSwarm, and then migrates the material points in order to find their location in the material mesh. If the background mesh and material mesh are distributed across processes differently, a large communication cost may be incurred. Then, each point has its material properties assigned using a loop on the CPU, which may also be expensive for large numbers of points. Finally, the background mesh is restored and points are migrated back to their original processes, incurring a similar communication overhead.

While these costs may seem quite dire, this procedure is performed only once per simulation. Amortized over the entire simulation, the initialization cost may be relatively minor.

Command-Line Interface¶

In order to use a material mesh, the option -mpm_use_material_mesh must be provided. Additional options, including specification of material regions, are given in MPM Material Mesh Options.

Table 16 MPM Material Mesh Options¶
Option	Description	Default Value
`-method mpm`	Required to enable MPM
`-mpm_point_location_type [uniform \| gauss \| random]`	Initial per-cell point distribution within each background mesh cell.	`uniform`
`-mpm_num_points_per_cell [int]`	Number of material points to place in each cell. For `uniform` or `gauss` point location type, this will be rounded to the nearest cube. Required.
`-mpm_use_material_mesh`	Required to enable the use of a material mesh.
`-mpm_material_mesh_dm_plex_filename [path to material mesh]`	Path to material mesh file.
`-mpm_material_mesh_dm_plex_hash_location`	Enable spatial hashing for point location. This is highly recommended to improve performance.
`-mpm_material_mesh_[dm options]`	Any valid option for the PETSc `DM` or `DMPlex` objects can be used to specify the material mesh.
`-material [string list]`	List of names to use as labels for each material. If not provided, a single material is assumed and can be specified with the prefix `-mpm_[material options]`.
`-mpm_{material name}_model`	Material to use (`elasticity-mpm-neo-hookean-current`, `elasticity-mpm-neo-hookean-damage-current`) for a given material. This model must be the same as the model used for the background mesh and all other materials. Required.
`-mpm_{material name}_label_name`	Domain label in the material mesh specifying the type of volume to use for specifying materials.	`"Cell Sets"`
`-mpm_{material name}_label_value [int]`	Domain value in the material mesh specifying the volume to use for a given material.	1
`-mpm_{material name}_{material options}`	Material options for the specified model. See Material Properties for more details.

The background mesh is specified via the typical -dm_* options. It is often favorable to choose a background mesh which is either regular or close to regular to improve solver performance and accuracy.

In addition to defining the materials for the material mesh, a “dummy” material must also be defined for the background mesh. The dummy material requires that the following is set:

The -model for the material should be the same as the model used for the material mesh.
Any global options (mainly -use_AT1 and -use_offdiagonal for the elasticity-mpm-neo-hookean-damage-current model) should be set on the background mesh. Values provided to the material mesh materials will be ignored.

Additional material properties on the background mesh material will be ignored, and input validation on these properties can be skipped via the -skip_checks flag.

See the below example for more details.

Example¶

Voxel Data¶

The voxel interface provides a compromise between the simplicity of the conformal mesh and the flexibility of the material mesh. This interface allows material label values to be assigned to points based on a regular grid of cubes. Each cube has the same side length \(s\), specified with -mpm_voxel_pixel_size, and the entire voxel data is assumed to cover the domain \([0,n_x \cdot s] \times [0,n_y \cdot s] \times [0,n_z \cdot s]\).

The voxel data is specified via -mpm_voxel_filename [path], which should be a file in the format:

DIM N_X N_Y N_Z

VALUE_X1_Y1_Z1 VALUE_X1_Y1_Z2 ... VALUE_X1_Y1_Z(N_Z)
VALUE_X1_Y2_Z1 VALUE_X1_Y2_Z2 ... VALUE_X1_Y2_Z(N_Z)
...
VALUE_X1_Y(N_Y)_Z1 VALUE_X1_Y(N_Y)_Z2 ... VALUE_X1_Y(N_Y)_Z(N_Z)

VALUE_X2_Y1_Z1 VALUE_X2_Y1_Z2 ... VALUE_X2_Y1_Z(N_Z)
VALUE_X2_Y2_Z1 VALUE_X2_Y2_Z2 ... VALUE_X2_Y2_Z(N_Z)
...
VALUE_X2_Y(N_Y)_Z1 VALUE_X2_Y(N_Y)_Z2 ... VALUE_X2_Y(N_Y)_Z(N_Z)

...

VALUE_X(N_X)_Y1_Z1 VALUE_X(N_X)_Y1_Z2 ... VALUE_X(N_X)_Y1_Z(N_Z)
VALUE_X(N_X)_Y2_Z1 VALUE_X(N_X)_Y2_Z2 ... VALUE_X(N_X)_Y2_Z(N_Z)
...
VALUE_X(N_X)_Y(N_Y)_Z1 VALUE_X(N_X)_Y(N_Y)_Z2 ... VALUE_X(N_X)_Y(N_Y)_Z(N_Z)

Extra whitespace in the voxel data is ignored to allow for human-readable formatting of small voxel data files. Each VALUE_X(i)_Y(j)_Z(k) is an integer value corresponding to the -mpm_[material name]_label_value of the material for that voxel. If a value does not have a corresponding material, it is automatically removed from the simulation. This allows for the embedding of more complex structures which are surrounded by void.

The primary use of the voxel data interface is to read voxelized and segmented data which has been assigned a numerical index based on its material. For a grayscale TIFF stack, the corresponding Ratel voxel data can be generated using Script to convert a TIFF stack into voxel data which can be read by Ratel and the numpy and tifffile packages.

Listing 1 Script to convert a TIFF stack into voxel data which can be read by Ratel¶

import numpy as np
from tifffile import TiffFile, TiffSequence
from pathlib import Path

def load_tiff_stack(tiffstack: Path) -> np.ndarray:
    """
    Load a stack of TIFF files into a 3D numpy array.
    """
    # Load the voxel data
    if tiffstack.is_dir():
        # If the input is a directory, load all TIFF files in the directory
        with TiffSequence(f'{Path(tiffstack) / "*.tif"}') as tif:
            # Read all images
            voxel_data = tif.asarray()
    else:
        with TiffFile(tiffstack) as tif:
            voxel_data = tif.asarray()

    # Move the axes to (x, y, z) format
    voxel_data = np.moveaxis(voxel_data, [0, 1, 2], [2, 1, 0])
    return voxel_data

def dump(dst: Path, voxel_data: np.ndarray, keep_grain_ids: bool = False) -> None:
    """
    Dump the tiff stack array to a file.
    """
    if not keep_grain_ids:
        # Remove the grain IDs from the voxel data
        voxel_data[voxel_data >= 2] = 2
    # Save the voxel data to a file
    np.savetxt(dst, voxel_data.flatten().reshape((1, voxel_data.size)),
               header=f'{len(voxel_data.shape)} ' + ' '.join(map(str, voxel_data.shape)),
               fmt='%d', delimiter=' ', comments='')

if __name__ == '__main__':
    tiff_file = Path("/path/to/tif/stack")
    output_file = Path("/path/to/output.dat")
    data = load_tiff_stack(tiff_file)
    dump(output_file, data)

Command-Line Options¶

The -mpm_use_voxel option is required to enable voxel data material point initialization. See MPM Voxel Options for a full list of options.

Table 17 MPM Voxel Options¶
Option	Description	Default Value
`-method mpm`	Required to enable MPM
`-mpm_point_location_type [uniform \| gauss \| random]`	Initial per-cell point distribution within each background mesh cell.	`uniform`
`-mpm_num_points_per_cell [int]`	Number of material points to place in each cell. For `uniform` or `gauss` point location type, this will be rounded to the nearest cube. Required
`-mpm_use_voxel`	Required to enable the use of a voxel data.
`-mpm_voxel_pixel_size [double]`	Side length of each voxel cube. Required.
`-mpm_voxel_filename [path]`	Path to voxel data file. Required.
`-material [string list]`	List of names to use as labels for each material. If not provided, a single material is assumed and can be specified with the prefix `-mpm_[material options]`.
`-mpm_{material name}_model`	Material to use (`elasticity-mpm-neo-hookean-current`, `elasticity-mpm-neo-hookean-damage-current`) for a given material. This model must be the same as the model used for the background mesh and all other materials. Required.
`-mpm_{material name}_label_value [int]`	Voxel value for which this material should be used.	1
`-mpm_{material name}_{material options}`	Material options for the specified model. See Material Properties for more details.

Similar to the material mesh interface, the voxel interface requires that a “dummy” material also be defined for the background mesh. The dummy material requires that the following is set:

The -model for the material should be the same as the model used for the material mesh.
Any global options (mainly -use_AT1 and -use_offdiagonal for the elasticity-mpm-neo-hookean-damage-current model) should be set on the background mesh. Values provided to the material mesh materials will be ignored.

Additional material properties on the background mesh material will be ignored, and input validation on these properties can be skipped via the -skip_checks flag.

Example¶

Example - Initialization from Voxel Data

An example of initialization from voxel data can be run with the command

$ ./bin/ratel-quasistatic examples/ymls/ex02-quasistatic-elasticity-mpm-neo-hookean-damage-current-voxel.yml

This example simulates the same problem as the material mesh example above: a multi-material model with a cylindrical rod embedded in a binder. The binder and rod are defined using a voxel data file, with corresponding labels of 1 and 2, respectively. As shown below, the voxel data file consists of \(n_x\) slices of voxel data, each of which has size \(n_y \times n_z\). For brevity, we show only the first slice, which consists of pure binder, and a slice from the center, which shows both materials.

The background mesh is an affine hexahedral mesh with the same dimensions as the material mesh. The “dummy” material only sets the model and the global -use_AT1 option.

# Ratel problem specification
# Dynamic material point material assignment using voxel data
# Neo-Hookean hyperelastic material at finite strain with damage, MPM

order: 1

method: mpm

mpm:
  point_location_type: uniform
  num_points_per_cell: 27
  use_voxel:
  voxel_filename: examples/data/ex02-quasistatic-elasticity-mpm-neo-hookean-damage-current-voxel.dat
  voxel_pixel_size: 0.05 # Each voxel is a cube with side length 0.05
  material: rod,binder
  binder:
    model: elasticity-mpm-neo-hookean-damage-current
    rho: 1.5
    E: 2.8
    nu: 0.3
    fracture_toughness: 1
    characteristic_length: 1.2
    damage_scaling: 1
    residual_stiffness: 1e-4
    label:
      value: 1
  rod:
    model: elasticity-mpm-neo-hookean-damage-current
    rho: 1
    E: 1
    nu: 0.3
    fracture_toughness: 1.6
    characteristic_length: 0.5
    damage_scaling: 1
    residual_stiffness: 1e-4
    label:
      value: 2

dm_plex:
  dim: 3
  box_faces: 3,3,3
  box_upper: 1,1,0.25
  simplex: 0
  hash_location:

material: dummy
dummy:
  model: elasticity-mpm-neo-hookean-damage-current
  use_AT1: false
  skip_checks:

bc:
  slip: 1,2
  slip_1:
    components: 0,1,2
    translate: 0,0,0
  slip_2:
    components: 0,1,2
    translate: 0,0,-0.025

remap:
  direction: z
  scale: 0.9

snes:
  linesearch:
    type: bt

pc:
  type: jacobi

expected_strain_energy: 3.854621021748e-03
diagnostic_components: damage
expected_max_diagnostic: 1.175983078375e-02