Functions
list[dict[str, float]]	read_sites (str\|os.PathLike[str] csv=CSV)

Deck	build_deck (str\|os.PathLike[str] output_path="runs/out/", *str preset="short", float\|None final_time=None, int\|None num_steps=None, int\|None output_interval=None, float omega=OMEGA, str\|os.PathLike[str] mesh_dir=MESH_DIR, str\|os.PathLike[str] csv=CSV)

int	main (list[str]\|None argv=None)

None	validate_sites (list[dict[str, float]] rows)

dict[str, object]	pack_geometry (int ncols=NCOLS, int nrows=NROWS, float r=R, float mesh_size=MESH_SIZE)

Deck	build_deck (str\|os.PathLike[str] output_path="runs/", *int ncols=NCOLS, int nrows=NROWS, float final_time=FINAL_TIME, int num_steps=NUM_STEPS, float mesh_size=MESH_SIZE, int search_interval=SEARCH_INTERVAL, str\|os.PathLike[str]\|None mesh_dir=None, bool write_mesh=True, bool in_process_mesh=False)

Deck	build_deck (str\|os.PathLike[str] output_path="runs/", *float final_time=FINAL_TIME, int num_steps=NUM_STEPS, float mesh_size=MESH_SIZE, float\|None horizon=None, str\|os.PathLike[str]\|None mesh_dir=None, tuple[str, str]\|None mesh_files=None, bool write_mesh=True, bool in_process_mesh=False)

float	_nu (float E_, float K_)

list[list[float]]	notch_void_boxes (float h, float notch_half, float notch_w, float notch_depth, float z_lo, float z_hi)

Deck	build_deck (str\|os.PathLike[str] output_path="runs/out/", *bool dim3=False, bool quick=False, float\|None final_time=None, float\|None dt=None, int\|None num_steps=None, str\|os.PathLike[str]\|None mesh_dir=None, float v_impact=IMPACT_V)

int	seed_prenotch (sim, dict[str, float] info)

np.ndarray	run_with_arrival_times (sim, dict[str, float] info, *float threshold=0.30)

tuple[float, int, float]\|None	crack_speed (sim, np.ndarray arrival, dict[str, float] info)

Deck	build_deck (str\|os.PathLike[str] output_path="out/", *float final_time=FINAL_TIME, int num_steps=NUM_STEPS, float\|None mesh_size=None, float\|None horizon=None, bool damping_on=False, float eps=EPSILON, bool zero_ic=False, str\|os.PathLike[str]\|None mesh_dir=None, bool write_mesh=True, bool in_process_mesh=False)

Deck	build_deck (str\|os.PathLike[str] output_path="runs/", *float final_time=FINAL_TIME, int num_steps=NUM_STEPS, int output_interval=1000, str\|os.PathLike[str]\|None mesh=None, float\|None mesh_size=None, list[str]\|None tags=None)

Deck	build_deck (str\|os.PathLike[str] output_path="runs/", *float final_time=FINAL_TIME, int num_steps=NUM_STEPS, int output_interval=2000, float mesh_size=MESH_SIZE, float horizon=HORIZON, list[str]\|None tags=None)

Variables
	HERE = Path(__file__).resolve().parent

str	MESH_DIR = HERE / "meshes"

str	CSV = HERE / "particle_locations_0.csv"

float	HORIZON = 6.0e-4

float	DENSITY = 1200.0

float	R_IN = 0.02

	R_SMALL

	R_LARGE

float	OMEGA = -20.0 * math.pi

float	GRAVITY = 10.0

dict	R_REF = {0: R_SMALL, 1: R_SMALL, 2: R_SMALL, 3: R_LARGE, 4: R_LARGE, 5: R_LARGE}

list	MESH_FILES

list	WALL_PARAMS

dict	KN

dict	K_MAT = {0: 1.0e4, 1: 1.0e5, 2: 1.0e5}

float	KN_FACTOR = 1.0

list	GRAIN_SHAPES

dict	PRESETS

list	TAGS

float	MESH_SIZE = R_SMALL / 5.0

float	R_OUT = R_IN + 1.5 * MESH_SIZE

float	L_BAR = 0.005

float	W_BAR = 1.5 * MESH_SIZE

tuple	ROT_CENTER = (-0.2 * R_IN, 0.2 * R_IN, 0.0)

float	BETA_N_FACTOR = 100.0

float	EPSILON = 0.95

float	CONTACT_RADIUS_FACTOR = 0.95

float	FRICTION_COEFF = 0.5

int	SEARCH_INTERVAL = 40

float	SEARCH_FACTOR = 10.0

float	R = 0.001

float	K = 2.16e7

float	NU = 0.25

float	GC = 50.0

float	WALL_VY = -0.06

	NCOLS

	NROWS

float	FINAL_TIME = 0.004

int	NUM_STEPS = 20000

float	W = 0.0008

float	H = 0.0012

	A_OUT

	B_OUT

	A_IN

	B_IN

float	ELL_THETA = 0.0

float	TIP_GAP = 0.00008

	RHO_T

	K_T

	NU_T

	GC_T

	RHO_E

	K_E

	NU_E

	GC_E

float	R_CONTACT_FACTOR = 0.95

float	IC_VY = -2.5

str	MESH_TRI = "mesh_triangle.msh"

str	MESH_ELL = "mesh_hollow_ellipse.msh"

float	NOTCH_DEPTH = 0.050

float	NOTCH_HALF = 0.025

float	NOTCH_W = 0.0015

float	THICKNESS = 0.009

	IW

	IH

float	IMPACT_V = 32.0

float	IMPACTOR_MASS = 1.57

float	RC_FACTOR = 0.95

float	RHO = 8000.0

float	E = 191.0e9

float	K_BULK = 159.2e9

float	DT = 2.5e-9

	QUICK

float	R1 = 0.001

float	PARTICLE_DIST = 0.001

float	RADIUS = 0.003

float	PULL_RATE = 0.005

str	MESH_FILE = "mesh_cir_1_0.msh"

float	SIDE = 0.01

Detailed Description

Attrition sim1: a rotating drum with an inward protrusion.

Circles, triangles and drums in two sizes tumble inside a cylinder that is
spun by a rotation displacement BC; the protrusion grinds them. This is the
deck behind ``attrition_test_sim1.gif`` in the top-level README.

The deck is assembled through :class:`peridem.Deck`. The wall is a complex
geometry, an outer circle with an inner circle removed and a protrusion
rectangle added. It is placed at its signed-volume centroid, computed by the
geometry object. The mesh axis is at the origin, so any other site translates
the mesh.

``gen_input.py`` next to this file writes the same decks as JSON for
``bin/PeriDEM``. ``python/tests/test_example_parity.py`` compares the two.

Attrition sim2: thin rotating container with an off-centre spin axis.

Circles, triangles, drums and hexagons in two sizes inside a thin-walled
cylinder with an inward bar. The container spins twice as fast as sim1 and
about a point offset from the origin, which throws the pack against the bar.

The deck is assembled through :class:`peridem.Deck`. Every parameter whose
default changes the result is set here. ``INPUT_DEFAULTS.md`` in this folder
records which those are.

Compression of a small circular-grain pack, set up in Python.

Parameters from Jha et al., J. Mech. Phys. Solids 151 (2021) 104376, section
4.3: material M1, lc = R/5, horizon = 3 lc, Rc = 0.95 h, C-bar = 100 and plate
velocity -0.06 m/s. The pack here is 4 by 3 grains and not the 502 of the
paper. The grain positions, the open channel that contains them and the moving
plate are derived from R and the pack size, and no deck file is read.

``test/test_data/peridem/jha2021_comp_n50/main.cpp`` builds the same deck in
C++. ``python/tests/test_example_parity.py`` compares the two decks key by
key.

Hollow ellipse dropped onto a short tip-up triangle, set up in Python.

``main.cpp`` in this folder builds the same deck in C++, with the same
geometry, materials and contact parameters.
``python/tests/test_example_parity.py`` compares the two decks key by key, so
that the comparison of the two runs is a comparison of the interfaces and not
of two problems.

Contact at the tip opens a crack. The ring separates into two pieces.

Silling 2003 Kalthoff-Winkler notched-plate impact, 2D and 3D, from Python.

A 200 x 100 mm maraging-steel plate with two open 1.5 mm notches is struck
edge-on by a rigid 1.57 kg cylinder at 32 m/s. Silling reports cracks running
from the notch tips at roughly 900 m/s and about 70 degrees to the notch.

This example uses the parts of the interface a pure JSON deck cannot reach:

* ``sim.setup()`` builds the model, then ``sim.break_bonds_in_slots(...)``
  seeds the pre-notch by cutting the peridynamic bonds that span each slot --
  the same operation the C++ driver does between ``init()`` and the time loop;
* ``sim.integrate()`` then runs the loop *without* re-initialising, or
  :func:`run_with_arrival_times` drives the loop step by step from Python to
  record when each node first becomes damaged, which gives the crack speed.

``Test_PeriDEM_notched_impact_inbuilt`` builds the same deck in C++.
``python/tests/test_example_parity.py`` compares the two decks key by key.

Two deformable circles: one fixed, one dropped onto it.

Two grains, one contact pair, gravity and an initial velocity. Both meshes are
built by Gmsh in the calling process and nothing is read from disk.

``test/test_data/peridem/twop_circ_inbuilt/main.cpp``, compiled with
TWOP_CONTACT_EXAMPLE as the target example_twop_circ_contact, builds the same
deck in C++.

Single circular particle: fixed SW patch, linear pull on the NE patch.

The deck is built in Python. ``input.json`` next to this file states the same
problem as a deck file. ``python/tests/test_example_parity.py`` runs the Python
deck through ``bin/PeriDEM`` and through the in-process interface and compares
the node fields.

Single rectangular particle pulled diagonally.

The south-west corner is fixed and the north-east corner is pulled at a
constant rate. The deck is built in Python and the mesh is generated in the
calling process, so no ``.msh`` file is read or written.

Function Documentation

◆ _nu()

float problem._nu	(	float	E_,
		float	K_
	)

protected

Definition at line 80 of file problem.py.

def _nu(E_: float, K_: float) -> float:
    return 0.5 * (1.0 - E_ / (3.0 * K_))

Referenced by build_deck().

Here is the caller graph for this function:

◆ build_deck() [1/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"out/"`,
		*float	final_time = `FINAL_TIME`,
		int	num_steps = `NUM_STEPS`,
		float \| None	mesh_size = `None`,
		float \| None	horizon = `None`,
		bool	damping_on = `False`,
		float	eps = `EPSILON`,
		bool	zero_ic = `False`,
		str \| os.PathLike[str] \| None	mesh_dir = `None`,
		bool	write_mesh = `True`,
		bool	in_process_mesh = `False`
	)

Definition at line 58 of file problem.py.

               in_process_mesh: bool = False) -> Deck:
    mesh_size = min(R1, R2) / 5.0 if mesh_size is None else mesh_size
    horizon = 3.0 * mesh_size if horizon is None else horizon
 
    E = to_E(K, NU)
    G = to_G(E=E, nu=NU)
 
    d = Deck(dim=2, t_final=final_time, n_steps=num_steps)
    d.set_output(output_path, tags=TAGS, interval=num_steps // 10, debug=2,
                 pvd_collection=True)
    d.set_gravity(0.0, -GRAVITY)
    d.set_neighbor(update_criteria="simple_all", s_factor=10.0, update_interval=40,
                   near_bd_nodes_tol=0.5)
 
    if in_process_mesh:
        meshes = [MeshSpec(size=mesh_size), MeshSpec(size=mesh_size)]
    else:
        base = Path(mesh_dir) if mesh_dir is not None else HERE / "inp"
        meshes = [MeshSpec(file=base / f"mesh_cir_{i}.msh", size=mesh_size,
                           write=write_mesh) for i in (1, 2)]
 
    g1 = d.add_particle_type(Geometry("circle", [R1, 0.0, 0.0, 0.0]), meshes[0])
    g2 = d.add_particle_type(Geometry("circle", [R2, 0.0, 0.0, 0.0]), meshes[1])
 
    m1 = d.add_material(horizon=horizon, density=DENSITY, K=K, G=G, Gc=GC,
                        influence_fn_type=1)
    m2 = d.add_material(horizon=horizon, density=DENSITY, K=K, G=G, Gc=GC,
                        influence_fn_type=1)
 
    kn = contact_stiffness(K, K, horizon)
    for i, j in ((0, 0), (0, 1), (1, 1)):
        d.add_contact_pair(i, j, contact_radius_factor=R_CONTACT_FACTOR, Kn=kn,
                           eps=eps, mu=FRICTION_COEFF, damping_on=damping_on,
                           friction_on=False, beta_n_factor=BETA_N_FACTOR,
                           K=K)
    d.set_contact_laws(damping_law="com_and_node",
                       friction_law="coulomb_simple")
 
    # Particle 0 is fixed and particle 1 falls onto it.
    d.add_displacement_bc(particles=[0], direction=[1, 2],
                          zero_displacement=True)
    vy = 0.0
    if not zero_ic:
        fallen = PARTICLE_DIST - horizon
        if fallen > 0.0:
            vy = -math.sqrt(2.0 * GRAVITY * fallen)
    d.add_initial_velocity([0.0, vy, 0.0], particles=[1])
 
    top_gap = PARTICLE_DIST
    d.place(R1, R1, 0.0, geometry=g1, material=m1, contact=0)
    d.place(R1, 2.0 * R1 + R2 + top_gap, 0.0, geometry=g2, material=m2,
            contact=1, theta=math.pi)
    return d
 
 

References build_deck().

Here is the call graph for this function:

◆ build_deck() [2/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"runs/"`,
		*float	final_time = `FINAL_TIME`,
		int	num_steps = `NUM_STEPS`,
		float	mesh_size = `MESH_SIZE`,
		float \| None	horizon = `None`,
		str \| os.PathLike[str] \| None	mesh_dir = `None`,
		tuple[str, str] \| None	mesh_files = `None`,
		bool	write_mesh = `True`,
		bool	in_process_mesh = `False`
	)

Assemble the deck.

``in_process_mesh`` keeps Gmsh output in memory and writes no ``.msh``,
which is what the README demo uses. The parity runs instead point both the
Python and the C++ side at the same ``.msh`` files under ``mesh_dir``.

Definition at line 66 of file problem.py.

               in_process_mesh: bool = False) -> Deck:
    """Assemble the deck.
 
    ``in_process_mesh`` keeps Gmsh output in memory and writes no ``.msh``,
    which is what the README demo uses. The parity runs instead point both the
    Python and the C++ side at the same ``.msh`` files under ``mesh_dir``.
    """
    horizon = 3.0 * mesh_size if horizon is None else horizon
 
    E_t = to_E(K_T, NU_T)
    G_t = to_G(E=E_t, nu=NU_T)
    E_e = to_E(K_E, NU_E)
    G_e = to_G(E=E_e, nu=NU_E)
 
    d = Deck(dim=2, t_final=final_time, n_steps=num_steps,
             bond_break="tension",
             # Without this, self-contact across broken bonds closes the
             # crack.
             self_contact="none")
 
    dt_out_n = max(1, num_steps // 40)
    d.set_output(output_path,
                 tags=["Displacement", "Velocity", "Force", "Damage_Z",
                       "Damage", "Particle_ID", "Fixity"],
                 interval=dt_out_n, debug=1, pvd_collection=True)
 
    d.set_gravity(0.0, -GRAVITY)
    d.set_neighbor(update_criteria="simple_all", s_factor=8.0, update_interval=5,
                   near_bd_nodes_tol=0.5)
 
    # --- particle types
    tri = Geometry("triangle", [-0.5 * W, 0.0, 0.0,
                                0.5 * W, 0.0, 0.0,
                                0.0, H, 0.0])
    ell = Geometry("ellipse_minus_ellipse",
                   [A_OUT, B_OUT, A_IN, B_IN, ELL_THETA, 0.0, 0.0, 0.0])
 
    if in_process_mesh:
        tri_mesh = MeshSpec(size=mesh_size)
        ell_mesh = MeshSpec(size=mesh_size)
    else:
        base = Path(mesh_dir) if mesh_dir is not None else HERE / "inp"
        names = mesh_files or (MESH_TRI, MESH_ELL)
        tri_mesh = MeshSpec(file=base / names[0], size=mesh_size,
                            write=write_mesh)
        ell_mesh = MeshSpec(file=base / names[1], size=mesh_size,
                            write=write_mesh)
 
    g_tri = d.add_particle_type(tri, tri_mesh)
    g_ell = d.add_particle_type(ell, ell_mesh)
 
    # --- materials
    m_tri = d.add_material(horizon=horizon, density=RHO_T, K=K_T, G=G_t,
                           Gc=GC_T, influence_fn_type=1)
    m_ell = d.add_material(horizon=horizon, density=RHO_E, K=K_E, G=G_e,
                           Gc=GC_E, influence_fn_type=1)
 
    # --- contact: group 0 is the triangle tip (wall), group 1 the ellipse
    d.add_contact_pair(0, 0, contact_radius_factor=R_CONTACT_FACTOR,
                       Kn=contact_stiffness(K_T, K_T, horizon), eps=0.95,
                       beta_n_factor=100.0, K=K_T)
    # Beta_n at the tip is large enough to start a crack, and damping is
    # on to limit the node velocities after contact.
    d.add_contact_pair(0, 1, contact_radius_factor=0.90,
                       Kn=contact_stiffness(K_T, K_E, horizon), eps=0.4,
                       beta_n_factor=8.0,
                       K=peridem.harmonic_mean(K_T, K_E))
    d.add_contact_pair(1, 1, contact_radius_factor=R_CONTACT_FACTOR,
                       Kn=contact_stiffness(K_E, K_E, horizon), eps=0.95,
                       beta_n_factor=100.0, K=K_E)
    d.set_contact_laws(damping_law="com_and_node",
                       friction_law="coulomb_simple")
 
    # --- boundary and initial conditions
    # The triangle is fixed in x and y at every node.
    d.add_displacement_bc(particles=[0], direction=[1, 2],
                          zero_displacement=True)
    d.add_initial_velocity([0.0, IC_VY, 0.0], particles=[1])
 
    # The triangle is placed at its centroid and the ellipse above the tip.
    tip_y = H
    ell_cy = tip_y + TIP_GAP + B_OUT
    cx, cy, cz = tri.center
    d.place(cx, cy, 0.0, geometry=g_tri, material=m_tri, contact=0, wall=True)
    d.place(0.0, ell_cy, 0.0, geometry=g_ell, material=m_ell, contact=1)
 
    return d
 
 

References build_deck().

Here is the call graph for this function:

◆ build_deck() [3/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"runs/"`,
		*float	final_time = `FINAL_TIME`,
		int	num_steps = `NUM_STEPS`,
		int	output_interval = `1000`,
		str \| os.PathLike[str] \| None	mesh = `None`,
		float \| None	mesh_size = `None`,
		list[str] \| None	tags = `None`
	)

Assemble the deck.

Pass ``mesh_size`` to have Gmsh build the disc in memory (no ``.msh`` on
disk); otherwise the committed ``mesh_cir_1_0.msh`` is read, which is what
the parity run against ``bin/PeriDEM`` uses.

Definition at line 53 of file problem.py.

               tags: list[str] | None = None) -> Deck:
    """Assemble the deck.
 
    Pass ``mesh_size`` to have Gmsh build the disc in memory (no ``.msh`` on
    disk); otherwise the committed ``mesh_cir_1_0.msh`` is read, which is what
    the parity run against ``bin/PeriDEM`` uses.
    """
    d = Deck(dim=2, t_final=final_time, n_steps=num_steps,
             particle_sim_type="Single_Particle")
    d.set_comment("Single-particle circle: fixed SW patch, linear pull on NE "
                  "(PDState), set up in Python")
    d.set_output(output_path, tags=tags or TAGS, interval=output_interval,
                 debug=1, tag_pp="0", pvd_collection=True)
 
    if mesh_size is not None:
        mesh_spec = MeshSpec(size=mesh_size)
    else:
        mesh_spec = MeshSpec(file=Path(mesh) if mesh is not None
                             else HERE / MESH_FILE)
    d.add_particle_type(Geometry("circle", [RADIUS, 0.0, 0.0, 0.0]), mesh_spec)
    d.add_material(horizon=HORIZON, density=DENSITY, K=K,
                   G=to_G(E=peridem.to_E(K, NU), nu=NU), Gc=GC,
                   influence_fn_type=1)
 
    # A square in the third quadrant is held and an equal square in the
    # first quadrant is pulled.
    hold = Geometry("rectangle", [-0.001, -0.001, 0.0, -0.0005, -0.0005, 0.0])
    pull = Geometry("rectangle", [0.0005, 0.0005, 0.0, 0.001, 0.001, 0.0])
    d.add_displacement_bc(region=hold, direction=[1, 2],
                          time_fn_type="constant", time_fn_params=[0.0],
                          spatial_fn_type="constant", zero_displacement=True)
    d.add_displacement_bc(region=pull, direction=[1, 2],
                          time_fn_type="linear", time_fn_params=[PULL_RATE],
                          spatial_fn_type="constant")
 
    d.set_test("test_peridynamics")
    return d
 
 

References build_deck().

Here is the call graph for this function:

◆ build_deck() [4/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"runs/"`,
		*float	final_time = `FINAL_TIME`,
		int	num_steps = `NUM_STEPS`,
		int	output_interval = `2000`,
		float	mesh_size = `MESH_SIZE`,
		float	horizon = `HORIZON`,
		list[str] \| None	tags = `None`
	)

Definition at line 51 of file problem.py.

               tags: list[str] | None = None) -> Deck:
    d = Deck(dim=2, t_final=final_time, n_steps=num_steps,
             particle_sim_type="Single_Particle")
    d.set_comment("Single-particle rectangle: uniform in-process mesh; fixed "
                  "SW corner, linear pull on NE, set up in Python")
    d.set_output(output_path, tags=tags or TAGS, interval=output_interval,
                 debug=1, tag_pp="1", pvd_collection=True)
 
    # The uniform mesh is generated in memory and no file is read.
    d.add_particle_type(Geometry("rectangle", [0.0, 0.0, 0.0, SIDE, SIDE, 0.0]),
                        MeshSpec(size=mesh_size, info="uniform"))
    d.add_material(horizon=horizon, density=DENSITY, K=K,
                   G=to_G(E=peridem.to_E(K, NU), nu=NU), Gc=GC,
                   influence_fn_type=1)
 
    corner = 0.2 * SIDE
    hold = Geometry("rectangle", [0.0, 0.0, 0.0, corner, corner, 0.0])
    pull = Geometry("rectangle",
                    [SIDE - corner, SIDE - corner, 0.0, SIDE, SIDE, 0.0])
    d.add_displacement_bc(region=hold, direction=[1, 2],
                          time_fn_type="constant", time_fn_params=[0.0],
                          spatial_fn_type="constant", zero_displacement=True)
    d.add_displacement_bc(region=pull, direction=[1, 2],
                          time_fn_type="linear", time_fn_params=[PULL_RATE],
                          spatial_fn_type="constant")
 
    d.set_test("test_peridynamics")
    return d
 
 

References build_deck().

Here is the call graph for this function:

◆ build_deck() [5/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"runs/"`,
		*int	ncols = `NCOLS`,
		int	nrows = `NROWS`,
		float	final_time = `FINAL_TIME`,
		int	num_steps = `NUM_STEPS`,
		float	mesh_size = `MESH_SIZE`,
		int	search_interval = `SEARCH_INTERVAL`,
		str \| os.PathLike[str] \| None	mesh_dir = `None`,
		bool	write_mesh = `True`,
		bool	in_process_mesh = `False`
	)

Definition at line 91 of file problem.py.

               in_process_mesh: bool = False) -> Deck:
    g = pack_geometry(ncols, nrows, R, mesh_size)
    horizon = float(g["horizon"])
    sites = g["sites"]
    n_pack = len(sites)
 
    E = to_E(K, NU)
    G = to_G(E=E, nu=NU)
    Kn = contact_stiffness(K, K, horizon)
 
    d = Deck(dim=2, t_final=final_time, n_steps=num_steps)
    d.set_comment("jha2021_comp_contact")
 
    dt_out = max(1, num_steps // 4)
    d.set_output(output_path, tags=TAGS, interval=dt_out, debug=1,
                 dt_test_out=max(1, dt_out // 10), tag_pp="0",
                 pvd_collection=True)
    d.set_gravity(0.0, -GRAVITY)
    d.set_neighbor(update_criteria="simple_all", s_factor=5.0, update_interval=search_interval,
                   near_bd_nodes_tol=0.5)
 
    grain = Geometry("circle", [R, 0.0, 0.0, 0.0])
    cup = Geometry("open_rect_channel_2d", list(g["cup"]))
    plate = Geometry("rectangle", list(g["plate"]))
 
    if in_process_mesh:
        meshes = [MeshSpec(size=mesh_size) for _ in range(3)]
    else:
        base = Path(mesh_dir) if mesh_dir is not None else HERE
        meshes = [MeshSpec(file=base / name, size=mesh_size, write=write_mesh)
                  for name in MESH_FILES]
 
    g_grain = d.add_particle_type(grain, meshes[0])
    g_cup = d.add_particle_type(cup, meshes[1])
    g_plate = d.add_particle_type(plate, meshes[2])
 
    m_grain = d.add_material(horizon=horizon, density=DENSITY, K=K, G=G, Gc=GC,
                             influence_fn_type=1)
    m_wall = d.add_material(horizon=horizon, density=DENSITY, K=K, G=G, Gc=GC,
                            influence_fn_type=1)
 
    for i, j in ((0, 0), (0, 1), (1, 1)):
        d.add_contact_pair(i, j, contact_radius_factor=0.95, Kn=Kn, eps=0.95,
                           beta_n_factor=100.0, K=K)
 
    # Particle ids: 0..n_pack-1 grains, then the cup, then the plate.
    id_cup, id_plate = n_pack, n_pack + 1
    d.add_displacement_bc(particles=[id_cup], direction=[1, 2],
                          zero_displacement=True)
    d.add_displacement_bc(particles=[id_plate], direction=[2],
                          time_fn_type="linear", time_fn_params=[WALL_VY],
                          spatial_fn_type="constant")
 
    for x, y in sites:
        d.place(x, y, geometry=g_grain, material=m_grain, contact=0)
    cx, cy, _ = cup.center
    d.place(cx, cy, geometry=g_cup, material=m_wall, contact=1, wall=True)
    px, py, _ = plate.center
    d.place(px, py, geometry=g_plate, material=m_wall, contact=1, wall=True)
 
    d.set_test("compressive_test", wall_id=id_plate, wall_force_direction=2)
    return d
 
 

References build_deck(), and pack_geometry().

Here is the call graph for this function:

◆ build_deck() [6/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"runs/out/"`,
		*bool	dim3 = `False`,
		bool	quick = `False`,
		float \| None	final_time = `None`,
		float \| None	dt = `None`,
		int \| None	num_steps = `None`,
		str \| os.PathLike[str] \| None	mesh_dir = `None`,
		float	v_impact = `IMPACT_V`
	)

Definition at line 95 of file problem.py.

               v_impact: float = IMPACT_V) -> Deck:
    cfg = dict(W=W, H=H, notch_half=NOTCH_HALF, mesh_size=MESH_SIZE, Iw=IW,
               Ih=IH, rho=RHO, E=E, K_bulk=K_BULK, Gc=GC, dt=DT,
               final_time=FINAL_TIME)
    if quick:
        cfg.update(QUICK)
        cfg["notch_depth"] = 0.5 * cfg["H"]
        cfg["notch_w"] = max(cfg["mesh_size"], 0.1 * cfg["notch_half"])
    else:
        cfg["notch_depth"] = NOTCH_DEPTH
        cfg["notch_w"] = NOTCH_W
 
    if final_time is not None:
        cfg["final_time"] = final_time
    if dt is not None:
        cfg["dt"] = dt
    n_steps = (int(round(cfg["final_time"] / cfg["dt"])) if num_steps is None
               else num_steps)
 
    w, h = cfg["W"], cfg["H"]
    mesh_size = cfg["mesh_size"]
    horizon = 3.0 * mesh_size
    notch_half, notch_w, notch_depth = (cfg["notch_half"], cfg["notch_w"],
                                        cfg["notch_depth"])
    iw, ih = cfg["Iw"], cfg["Ih"]
    gap = 1.5 * RC_FACTOR * mesh_size
    rho, e_mod, k_bulk, gc = cfg["rho"], cfg["E"], cfg["K_bulk"], cfg["Gc"]
    nu = _nu(e_mod, k_bulk)
    g_mod = e_mod / (2.0 * (1.0 + nu))
 
    d = Deck(dim=3 if dim3 else 2, t_final=cfg["final_time"], n_steps=n_steps,
             # Element-node connectivity is only used for strain output and
             # does not handle the hexahedra the 3D structured grid produces.
             populate_element_node_connectivity=not dim3,
             self_contact="none", bond_break="tension", wall_contact="meshed")
 
    d.set_output(output_path, tags=TAGS, interval=max(1, n_steps // 10),
                 debug=1, perform_fe_out=False, dt_test_out=n_steps,
                 pvd_collection=False)
    d.set_neighbor(update_criteria="simple_all", s_factor=5.0, update_interval=1,
                   near_bd_nodes_tol=0.5)
 
    # --- bodies
    if dim3:
        plate = Geometry("cuboid", [-0.5 * w, -0.5 * h, -0.5 * THICKNESS,
                                    0.5 * w, 0.5 * h, 0.5 * THICKNESS])
        # Cylinder axis along the impact direction, flat face striking the edge.
        impactor = Geometry("cylinder", [0.5 * iw, 0.0, -0.5 * ih, 0.0,
                                         0.0, ih, 0.0])
        z_lo, z_hi = -0.5 * THICKNESS - 1.0e-9, 0.5 * THICKNESS + 1.0e-9
    else:
        plate = Geometry("rectangle",
                         [-0.5 * w, -0.5 * h, 0.0, 0.5 * w, 0.5 * h, 0.0])
        impactor = Geometry("rectangle",
                            [-0.5 * iw, -0.5 * ih, 0.0, 0.5 * iw, 0.5 * ih, 0.0])
        z_lo, z_hi = -1.0e-9, 1.0e-9
 
    voids = notch_void_boxes(h, notch_half, notch_w, notch_depth, z_lo, z_hi)
    base = Path(mesh_dir) if mesh_dir is not None else HERE / "runs/inp"
    # The plate sits on Silling's equally spaced structured grid with the notch
    # slots carved out, so they are real gaps rather than cut material.
    g_plate = d.add_particle_type(
        plate, MeshSpec(file=base / "mesh_plate.msh", size=mesh_size,
                        info="uniform", voids=voids))
    # In 2D the impactor is a rectangle, so it goes on the same grid; in 3D it
    # is a cylinder, which a uniform grid cannot represent.
    g_imp = d.add_particle_type(
        impactor,
        MeshSpec(file=base / "mesh_impactor.msh", size=mesh_size,
                 info="gmsh_builtin_mesh" if dim3 else "uniform"))
 
    # The impactor has the plate material with Gc = 0. Its nodal forces are
    # replaced by the rigid-body acceleration, so its bonds carry no load.
    m_plate = d.add_material(material_type="PMBBond", horizon=horizon,
                             density=rho, K=k_bulk, G=g_mod, Gc=gc, E=e_mod,
                             influence_fn_type=0, influence_fn_params=[1.0])
    m_imp = d.add_material(material_type="PDElasticBond", horizon=horizon,
                           density=rho, K=k_bulk, G=g_mod, Gc=0.0, E=e_mod,
                           influence_fn_type=0, influence_fn_params=[1.0])
 
    # Contact force density is Kn * V_j * overlap, so Kn carries one power of
    # the horizon per spatial dimension of the nodal weight: 5 in 3D, 4 in 2D.
    kn = contact_stiffness(k_bulk, k_bulk, horizon,
                           horizon_power=5 if dim3 else 4)
    for i, j in ((0, 0), (0, 1), (1, 1)):
        d.add_contact_pair(i, j, contact_radius=RC_FACTOR * mesh_size,
                           Kn=kn, eps=1.0, damping_on=False, friction_on=False,
                           beta_n_factor=0.0, K=k_bulk)
    d.set_contact_laws(damping_law="off", friction_law="coulomb_simple")
 
    # The plate carries no load on its boundary. The displacement condition
    # constrains the impactor to move along y.
    d.add_displacement_bc(particles=[1], direction=[1],
                          time_fn_type="constant", time_fn_params=[0.0],
                          spatial_fn_type="constant", zero_displacement=True)
    d.add_initial_velocity([0.0, -v_impact, 0.0], particles=[1])
 
    # In two dimensions the mass is per unit thickness of the 9 mm plate.
    d.add_rigid_particle(1, IMPACTOR_MASS if dim3 else IMPACTOR_MASS / THICKNESS)
 
    d.place(0.0, 0.0, 0.0, geometry=g_plate, material=m_plate, contact=0)
    d.place(0.0, 0.5 * h + 0.5 * ih + gap, 0.0, geometry=g_imp,
            material=m_imp, contact=1)
 
    # Used by seed_prenotch and crack_speed below.
    d.extra_info = {"H": h, "notch_half": notch_half, "notch_w": notch_w,
                    "notch_depth": notch_depth, "mesh_size": mesh_size,
                    "horizon": horizon, "n_steps": n_steps,
                    "dt": cfg["dt"], "final_time": cfg["final_time"]}
    return d
 
 

References _nu(), and notch_void_boxes().

Here is the call graph for this function:

◆ build_deck() [7/7]

Deck problem.build_deck	(	str \| os.PathLike[str]	output_path = `"runs/out/"`,
		*str	preset = `"short"`,
		float \| None	final_time = `None`,
		int \| None	num_steps = `None`,
		int \| None	output_interval = `None`,
		float	omega = `OMEGA`,
		str \| os.PathLike[str]	mesh_dir = `MESH_DIR`,
		str \| os.PathLike[str]	csv = `CSV`
	)

Definition at line 125 of file problem.py.

               csv: str | os.PathLike[str] = CSV) -> Deck:
    t_default, n_default, out_default = PRESETS[preset]
    final_time = t_default if final_time is None else final_time
    num_steps = n_default if num_steps is None else num_steps
    output_interval = out_default if output_interval is None else output_interval
 
    d = Deck(dim=2, t_final=final_time, n_steps=num_steps,
             bond_break="tension", self_contact="none", wall_contact="meshed")
    d.set_output(output_path, tags=TAGS, interval=output_interval, debug=1,
                 perform_fe_out=False,
                 dt_test_out=max(1, output_interval // 10), tag_pp="0",
                 pvd_collection=True)
    d.set_gravity(0.0, -GRAVITY)
    d.set_neighbor(update_criteria="simple_all", s_factor=10.0, update_interval=40,
                   near_bd_nodes_tol=0.5)
 
    base = Path(mesh_dir)
    for i, (name, params) in enumerate(GRAIN_SHAPES):
        d.add_particle_type(Geometry(name, params),
                            MeshSpec(file=(base / MESH_FILES[i]).resolve()))
    wall = Geometry("complex", WALL_PARAMS,
                    vec_type=["circle", "circle", "rectangle"],
                    vec_flag=["plus", "minus", "plus"])
    g_wall = d.add_particle_type(wall,
                                 MeshSpec(file=(base / MESH_FILES[6]).resolve()))
 
    # Material 0 is the small grains, material 1 the large grains and drum.
    m_small = d.add_material(horizon=HORIZON, density=DENSITY, K=1.0e4, G=6.0e3,
                             Gc=50.0, influence_fn_type=1)
    m_large = d.add_material(horizon=HORIZON, density=DENSITY, K=1.0e5, G=6.0e4,
                             Gc=100.0, influence_fn_type=1)
 
    # Contact groups: 0 small grains, 1 large grains, 2 the drum.
    for (i, j), kn in KN.items():
        d.add_contact_pair(
            i, j, contact_radius_factor=0.95, Kn=kn,
            K=peridem.harmonic_mean(K_MAT[i], K_MAT[j]),
            damping_on=False, friction_on=False, eps=0.95, mu=0.5,
            Kn_factor=KN_FACTOR,
            # Damping is off, so Beta_n_Factor does not enter the force.
            raw={"Beta_n_Factor": 100.0})
    d.set_contact_laws(damping_law="off", friction_law="coulomb_simple",
                       correct_volume=False)
 
    sites = read_sites(csv)
    for s in sites:
        zone = int(s["zone"])
        large = zone >= 3
        d.place(s["x"], s["y"], s["z"], geometry=zone,
                material=m_large if large else m_small,
                contact=1 if large else 0,
                theta=s["theta"], scale=s["r"] / R_REF[zone])
 
    wcx, wcy, wcz = wall.center
    wall_id = len(sites)
    d.place(wcx, wcy, wcz, geometry=g_wall, material=m_large, contact=2,
            wall=True)
 
    # The rotation centre is the last three time function parameters.
    d.add_displacement_bc(particles=[wall_id], direction=[1, 2],
                          time_fn_type="rotation",
                          time_fn_params=[omega, 0.0, 0.0, 0.0],
                          spatial_fn_type="rotation")
    return d
 
 

References read_sites().

Referenced by build_deck(), build_deck(), build_deck(), build_deck(), build_deck(), crack_speed(), main(), and validate_sites().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ crack_speed()

tuple[float, int, float] \| None problem.crack_speed	(		sim,
		np.ndarray	arrival,
		dict[str, float]	info
	)

Least-squares crack-tip speed from the arrival-time field.

Fits distance-from-the-right-notch-tip against arrival time over damaged
plate nodes below the tip (where the crack runs) and within half the plate
height of it. Returns (m/s, points used, correlation).

The returned correlation states how well the fit holds. If damage spreads
through the plate instead of advancing as a front, the fit has no meaning
and the correlation is low. The ``--quick`` settings produce that, on a
plate of 40 by 20 mm with E = 1.23 GPa.

Definition at line 264 of file problem.py.

                info: dict[str, float]) -> tuple[float, int, float] | None:
    """Least-squares crack-tip speed from the arrival-time field.
 
    Fits distance-from-the-right-notch-tip against arrival time over damaged
    plate nodes below the tip (where the crack runs) and within half the plate
    height of it. Returns (m/s, points used, correlation).
 
    The returned correlation states how well the fit holds. If damage spreads
    through the plate instead of advancing as a front, the fit has no meaning
    and the correlation is low. The ``--quick`` settings produce that, on a
    plate of 40 by 20 mm with E = 1.23 GPa.
    """
    y_tip = 0.5 * info["H"] - info["notch_depth"]
    tip = np.array([info["notch_half"], y_tip])
    x = np.asarray(sim.reference)
    seen = (arrival >= 0.0) & (np.asarray(sim.particle_id) == 0)
    dist = np.linalg.norm(x[:, :2] - tip, axis=1)
    # The crack runs into the plate, below the notch tip.
    band = seen & (dist < 0.5 * info["H"]) & (x[:, 1] < y_tip)
    n = int(band.sum())
    if n < 20:
        return None
    t = arrival[band]
    r_dist = dist[band]
    if np.ptp(t) <= 0.0:
        return None
    slope, _ = np.polyfit(t, r_dist, 1)
    r = float(np.corrcoef(t, r_dist)[0, 1])
    return float(slope), n, r
 
 

References build_deck(), crack_speed(), run_with_arrival_times(), and seed_prenotch().

Referenced by crack_speed().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ main()

int problem.main ( list[str] | None argv = None )

Definition at line 198 of file problem.py.

def main(argv: list[str] | None = None) -> int:
    import argparse
 
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument("-o", "--output", default=str(HERE / "runs_py/out"))
    p.add_argument("--preset", choices=sorted(PRESETS), default="short")
    p.add_argument("--final-time", type=float, default=None)
    p.add_argument("--num-steps", type=int, default=None)
    p.add_argument("--output-interval", type=int, default=None)
    p.add_argument("--nthreads", type=int,
                   default=int(os.environ.get("NTHREADS", "4")))
    p.add_argument("--snapshot", metavar="PNG", default=None)
    p.add_argument("--write-deck", metavar="PATH", default=None)
    args = p.parse_args(argv)
 
    out = Path(args.output).resolve()
    d = build_deck(str(out) + "/", preset=args.preset,
                   final_time=args.final_time, num_steps=args.num_steps,
                   output_interval=args.output_interval)
    if args.write_deck:
        print(d.write(args.write_deck))
        return 0
 
    out.mkdir(parents=True, exist_ok=True)
    peridem.init(n_threads=args.nthreads)
    sim = d.run(workdir=HERE)
    if peridem.mpi_rank() != 0:
        return 0
    print(f"done: grains={sim.n_particles} walls={sim.n_walls} "
          f"nodes={sim.n_nodes} t={sim.time:g}")
    worst = max(sim.particles, key=lambda q: float(q.damage.max()), default=None)
    if worst is not None:
        print(f"most damaged grain: id={worst.id} "
              f"max Z={float(worst.damage.max()):.4f} "
              f"shape={worst.geometry_name}")
    if args.snapshot:
        print(peridem.snapshot(peridem.last_vtu(out), args.snapshot,
                               color="Damage_Z",
                               title="attrition sim1: damage"))
    print(f"VTU under {out}")
    return 0
 
 

References build_deck().

Here is the call graph for this function:

◆ notch_void_boxes()

list[list[float]] problem.notch_void_boxes	(	float	h,
		float	notch_half,
		float	notch_w,
		float	notch_depth,
		float	z_lo,
		float	z_hi
	)

The two notch slots as boxes, open at the top edge.

Definition at line 84 of file problem.py.

                     z_hi: float) -> list[list[float]]:
    """The two notch slots as boxes, open at the top edge."""
    hw = 0.5 * notch_w
    y_tip = 0.5 * h - notch_depth
    y_hi = 0.5 * h + 1.0e-9
    return [[-notch_half - hw, y_tip, z_lo, -notch_half + hw, y_hi, z_hi],
            [notch_half - hw, y_tip, z_lo, notch_half + hw, y_hi, z_hi]]
 
 

Referenced by build_deck().

Here is the caller graph for this function:

◆ pack_geometry()

dict[str, object] problem.pack_geometry	(	int	ncols = `NCOLS`,
		int	nrows = `NROWS`,
		float	r = `R`,
		float	mesh_size = `MESH_SIZE`
	)

Pack spacing and container dimensions, all derived from R and the mesh.

The gap between grains is slightly larger than the contact radius, so that
no pair is in contact at t = 0. Gravity and the plate bring pairs into
contact.

Definition at line 61 of file problem.py.

                  mesh_size: float = MESH_SIZE) -> dict[str, object]:
    """Pack spacing and container dimensions, all derived from R and the mesh.
 
    The gap between grains is slightly larger than the contact radius, so that
    no pair is in contact at t = 0. Gravity and the plate bring pairs into
    contact.
    """
    horizon = 3.0 * mesh_size
    h_est = 0.7 * mesh_size          # realized hmin on a Gmsh disc is ~0.7 lc
    rc_est = 0.95 * h_est
    padding = 1.15 * rc_est
    rwp = horizon + padding
    wall_t = rwp - padding
 
    step = 2.0 * r + padding
    lin = 2.0 * padding + 2.0 * r + (ncols - 1) * step
    win = 2.0 * padding + 2.0 * r + (nrows - 1) * step
 
    sites = [(padding + r + i * step, padding + r + j * step)
             for j in range(nrows) for i in range(ncols)]
    return {
        "horizon": horizon, "padding": padding, "rwp": rwp, "wall_t": wall_t,
        "Lin": lin, "Win": win, "sites": sites,
        # The fixed open channel and the plate that compresses the pack.
        "cup": [-rwp, -rwp, lin + rwp, win + wall_t, wall_t, 0.0],
        "plate": [-padding, win, 0.0, lin + padding, win + wall_t, 0.0],
    }
 
 

Referenced by build_deck().

Here is the caller graph for this function:

◆ read_sites()

list[dict[str, float]] problem.read_sites ( str | os.PathLike[str] csv = CSV )

Packed grain sites: ``zone, x, y, z, r, theta`` per row.

Definition at line 105 of file problem.py.

def read_sites(csv: str | os.PathLike[str] = CSV) -> list[dict[str, float]]:
    """Packed grain sites: ``zone, x, y, z, r, theta`` per row."""
    sites: list[dict[str, float]] = []
    with Path(csv).open() as f:
        next(f)
        for line in f:
            parts = [p.strip() for p in line.split(",")]
            if len(parts) < 6:
                continue
            zone = int(float(parts[0]))
            x, y, z, r, theta = (float(v) for v in parts[1:6])
            if math.hypot(x, y) + r > R_IN - 1.0e-6:
                raise ValueError(
                    f"initial pack overlaps the drum: zone={zone} at "
                    f"({x}, {y}) r={r} reaches past R_in={R_IN}")
            sites.append({"zone": zone, "x": x, "y": y, "z": z, "r": r,
                          "theta": theta})
    return sites
 
 

Referenced by build_deck(), and validate_sites().

Here is the caller graph for this function:

◆ run_with_arrival_times()

np.ndarray problem.run_with_arrival_times	(		sim,
		dict[str, float]	info,
		*float	threshold = `0.30`
	)

Drive the time loop from Python, recording first-damage time per node.

The arrival-time field is what gives the crack speed Silling reports.
Returns an array of length n_nodes, -1 where the node never damaged.

Definition at line 228 of file problem.py.

                           threshold: float = 0.30) -> np.ndarray:
    """Drive the time loop from Python, recording first-damage time per node.
 
    The arrival-time field is what gives the crack speed Silling reports.
    Returns an array of length n_nodes, -1 where the node never damaged.
    """
    n_nodes = sim.n_nodes
    arrival = np.full(n_nodes, -1.0)
    plate = np.asarray(sim.particle_id) == 0
    sample_every = max(1, int(info["n_steps"]) // 400)
 
    def sample() -> None:
        fresh = (arrival < 0.0) & plate & (sim.damage >= threshold)
        arrival[fresh] = sim.time
 
    sim.apply_initial_condition()
    if sim.perform_output:
        sim.write_output()
    sim.set_current_dt(sim.dt)
    sim.apply_displacement_bc()
    sim.compute_forces()
    sim.apply_rigid_body_constraint()
    while sim.step_index < sim.n_steps:
        sim.step()
        if sim.should_output:
            sim.write_output()
        if sim.step_index % sample_every == 0:
            sample()
        sim.check_stop()
        if sim.stopped:
            break
    sample()
    return arrival
 
 

Referenced by crack_speed().

Here is the caller graph for this function:

◆ seed_prenotch()

int problem.seed_prenotch	(		sim,
		dict[str, float]	info
	)

Cut the peridynamic bonds that span the two notch slots.

The slots are already absent from the mesh, but the horizon is twice the
slot width, so bonds still reach across them. Must run after ``setup()``.

Definition at line 212 of file problem.py.

def seed_prenotch(sim, info: dict[str, float]) -> int:
    """Cut the peridynamic bonds that span the two notch slots.
 
    The slots are already absent from the mesh, but the horizon is twice the
    slot width, so bonds still reach across them. Must run after ``setup()``.
    """
    y_top = 0.5 * info["H"]
    y_tip = y_top - info["notch_depth"]
    n = sim.break_bonds_in_slots(
        [-info["notch_half"], info["notch_half"]], info["notch_w"],
        y_tip, y_top + 0.01 * info["H"], 0)
    if n < 10:
        raise RuntimeError(f"expected pre-notch bonds to cut, got {n}")
    return n
 
 

Referenced by crack_speed().

Here is the caller graph for this function:

◆ validate_sites()

None problem.validate_sites ( list[dict[str, float]] rows )

Reject a pack that starts inside the wall, the bar, or another grain.

Definition at line 129 of file problem.py.

def validate_sites(rows: list[dict[str, float]]) -> None:
    """Reject a pack that starts inside the wall, the bar, or another grain."""
    bar = (R_IN - L_BAR, -0.5 * W_BAR, R_IN, 0.5 * W_BAR)
    errors: list[str] = []
    for i, a in enumerate(rows):
        if a["zone"] not in R_REF:
            errors.append(f"row {i}: unknown shape id {a['zone']}")
            continue
        if math.hypot(a["x"], a["y"]) + a["r"] > R_IN - 1.0e-9:
            errors.append(f"row {i}: reaches past R_in")
        if not (a["x"] + a["r"] < bar[0] or a["x"] - a["r"] > bar[2]
                or a["y"] + a["r"] < bar[1] or a["y"] - a["r"] > bar[3]):
            errors.append(f"row {i}: overlaps the protrusion")
    for i in range(len(rows)):
        for j in range(i + 1, len(rows)):
            a, b = rows[i], rows[j]
            if math.hypot(a["x"] - b["x"], a["y"] - b["y"]) < \
                    a["r"] + b["r"] - 1.0e-9:
                errors.append(f"rows {i},{j} overlap")
    if errors:
        raise ValueError("initial pack invalid:\n  " + "\n  ".join(errors[:40]))
 
 

References build_deck(), and read_sites().

Here is the call graph for this function:

Variable Documentation

◆ A_IN

problem.A_IN

Definition at line 44 of file problem.py.

◆ A_OUT

problem.A_OUT

Definition at line 43 of file problem.py.

◆ B_IN

problem.B_IN

Definition at line 44 of file problem.py.

◆ B_OUT

problem.B_OUT

Definition at line 43 of file problem.py.

◆ BETA_N_FACTOR

float problem.BETA_N_FACTOR = 100.0

Definition at line 85 of file problem.py.

◆ CONTACT_RADIUS_FACTOR

float problem.CONTACT_RADIUS_FACTOR = 0.95

Definition at line 87 of file problem.py.

◆ CSV

str problem.CSV = HERE / "particle_locations_0.csv"

Definition at line 45 of file problem.py.

◆ DENSITY

float problem.DENSITY = 1200.0

Definition at line 48 of file problem.py.

◆ DT

float problem.DT = 2.5e-9

Definition at line 68 of file problem.py.

◆ E

float problem.E = 191.0e9

Definition at line 64 of file problem.py.

◆ ELL_THETA

float problem.ELL_THETA = 0.0

Definition at line 45 of file problem.py.

◆ EPSILON

float problem.EPSILON = 0.95

Definition at line 86 of file problem.py.

◆ FINAL_TIME

float problem.FINAL_TIME = 0.004

Definition at line 50 of file problem.py.

◆ FRICTION_COEFF

float problem.FRICTION_COEFF = 0.5

Definition at line 88 of file problem.py.

◆ GC

float problem.GC = 50.0

Definition at line 45 of file problem.py.

◆ GC_E

problem.GC_E

Definition at line 51 of file problem.py.

◆ GC_T

problem.GC_T

Definition at line 49 of file problem.py.

◆ GRAIN_SHAPES

list problem.GRAIN_SHAPES

Initial value:

=  [
    ("circle", [R_SMALL, 0.0, 0.0, 0.0]),
    ("triangle", [R_SMALL, 0.0, 0.0, 0.0]),
    ("drum2d", [R_SMALL, 0.0004, 0.0, 0.0, 0.0]),
    ("circle", [R_LARGE, 0.0, 0.0, 0.0]),
    ("triangle", [R_LARGE, 0.0, 0.0, 0.0]),
    ("drum2d", [R_LARGE, 0.0012, 0.0, 0.0, 0.0]),
]

Definition at line 86 of file problem.py.

◆ GRAVITY

float problem.GRAVITY = 10.0

Definition at line 52 of file problem.py.

◆ H

problem.H = 0.0012

Definition at line 40 of file problem.py.

◆ HERE

problem.HERE = Path(__file__).resolve().parent

Definition at line 33 of file problem.py.

◆ HORIZON

float problem.HORIZON = 6.0e-4

Definition at line 47 of file problem.py.

◆ IC_VY

float problem.IC_VY = -2.5

Definition at line 55 of file problem.py.

◆ IH

problem.IH

Definition at line 55 of file problem.py.

◆ IMPACT_V

float problem.IMPACT_V = 32.0

Definition at line 56 of file problem.py.

◆ IMPACTOR_MASS

float problem.IMPACTOR_MASS = 1.57

Definition at line 57 of file problem.py.

◆ IW

problem.IW

Definition at line 55 of file problem.py.

◆ K

float problem.K = 2.16e7

Definition at line 43 of file problem.py.

◆ K_BULK

float problem.K_BULK = 159.2e9

Definition at line 65 of file problem.py.

◆ K_E

problem.K_E

Definition at line 51 of file problem.py.

◆ K_MAT

dict problem.K_MAT = {0: 1.0e4, 1: 1.0e5, 2: 1.0e5}

Definition at line 82 of file problem.py.

◆ K_T

problem.K_T

Definition at line 49 of file problem.py.

◆ KN

dict problem.KN

Initial value:

=  {
    (0, 0): 7.368284e20,
    (0, 1): 1.339688e21,
    (0, 2): 1.339688e21,
    (1, 1): 7.368284e21,
    (1, 2): 7.368284e21,
    (2, 2): 7.368284e21,
}

Definition at line 74 of file problem.py.

◆ KN_FACTOR

float problem.KN_FACTOR = 1.0

Definition at line 83 of file problem.py.

◆ L_BAR

float problem.L_BAR = 0.005

Definition at line 47 of file problem.py.

◆ MESH_DIR

str problem.MESH_DIR = HERE / "meshes"

Definition at line 44 of file problem.py.

◆ MESH_ELL

str problem.MESH_ELL = "mesh_hollow_ellipse.msh"

Definition at line 63 of file problem.py.

◆ MESH_FILE

str problem.MESH_FILE = "mesh_cir_1_0.msh"

Definition at line 50 of file problem.py.

◆ MESH_FILES

tuple problem.MESH_FILES

Initial value:

=  [
    "mesh_cir_small_0.msh",
    "mesh_tri_small_0.msh",
    "mesh_drum2d_small_0.msh",
    "mesh_cir_large_0.msh",
    "mesh_tri_large_0.msh",
    "mesh_drum2d_large_0.msh",
    "mesh_wall_0.msh",
]

Definition at line 57 of file problem.py.

◆ MESH_SIZE

float problem.MESH_SIZE = R_SMALL / 5.0

Definition at line 43 of file problem.py.

◆ MESH_TRI

str problem.MESH_TRI = "mesh_triangle.msh"

Definition at line 62 of file problem.py.

◆ NCOLS

problem.NCOLS

Definition at line 49 of file problem.py.

◆ NOTCH_DEPTH

float problem.NOTCH_DEPTH = 0.050

Definition at line 49 of file problem.py.

◆ NOTCH_HALF

float problem.NOTCH_HALF = 0.025

Definition at line 50 of file problem.py.

◆ NOTCH_W

float problem.NOTCH_W = 0.0015

Definition at line 51 of file problem.py.

◆ NROWS

problem.NROWS

Definition at line 49 of file problem.py.

◆ NU

float problem.NU = 0.25

Definition at line 44 of file problem.py.

◆ NU_E

problem.NU_E

Definition at line 51 of file problem.py.

◆ NU_T

problem.NU_T

Definition at line 49 of file problem.py.

◆ NUM_STEPS

int problem.NUM_STEPS = 20000

Definition at line 51 of file problem.py.

◆ OMEGA

float problem.OMEGA = -20.0 * math.pi

Definition at line 51 of file problem.py.

◆ PARTICLE_DIST

float problem.PARTICLE_DIST = 0.001

Definition at line 39 of file problem.py.

◆ PRESETS

dict problem.PRESETS

Initial value:

=  {
    "short": (0.01, 100_000, 2000),
    "medium": (0.03, 300_000, 3000),
    "paper": (0.1, 1_000_000, 2500),
}

Definition at line 95 of file problem.py.

◆ PULL_RATE

float problem.PULL_RATE = 0.005

Definition at line 44 of file problem.py.

◆ QUICK

problem.QUICK

Initial value:

=  dict(W=0.040, H=0.020, notch_half=0.005, mesh_size=0.020 / 16.0,
             Iw=0.008, Ih=0.004, rho=1200.0, E=1.23e9, K_bulk=2.0e9, Gc=424.0,
             dt=1.0e-8, final_time=1.0e-4)

Definition at line 75 of file problem.py.

◆ R

float problem.R = 0.001

Definition at line 40 of file problem.py.

◆ R1

float problem.R1 = 0.001

Definition at line 38 of file problem.py.

◆ R_CONTACT_FACTOR

float problem.R_CONTACT_FACTOR = 0.95

Definition at line 53 of file problem.py.

◆ R_IN

float problem.R_IN = 0.02

Definition at line 49 of file problem.py.

◆ R_LARGE

problem.R_LARGE

Definition at line 50 of file problem.py.

◆ R_OUT

float problem.R_OUT = R_IN + 1.5 * MESH_SIZE

Definition at line 46 of file problem.py.

◆ R_REF

dict problem.R_REF = {0: R_SMALL, 1: R_SMALL, 2: R_SMALL, 3: R_LARGE, 4: R_LARGE, 5: R_LARGE}

Definition at line 55 of file problem.py.

◆ R_SMALL

problem.R_SMALL

Definition at line 50 of file problem.py.

◆ RADIUS

float problem.RADIUS = 0.003

Definition at line 35 of file problem.py.

◆ RC_FACTOR

float problem.RC_FACTOR = 0.95

Definition at line 60 of file problem.py.

◆ RHO

float problem.RHO = 8000.0

Definition at line 63 of file problem.py.

◆ RHO_E

problem.RHO_E

Definition at line 51 of file problem.py.

◆ RHO_T

problem.RHO_T

Definition at line 49 of file problem.py.

◆ ROT_CENTER

tuple problem.ROT_CENTER = (-0.2 * R_IN, 0.2 * R_IN, 0.0)

Definition at line 53 of file problem.py.

◆ SEARCH_FACTOR

float problem.SEARCH_FACTOR = 10.0

Definition at line 90 of file problem.py.

◆ SEARCH_INTERVAL

int problem.SEARCH_INTERVAL = 40

Definition at line 89 of file problem.py.

◆ SIDE

float problem.SIDE = 0.01

Definition at line 34 of file problem.py.

◆ TAGS

list problem.TAGS

Initial value:

1= ["Displacement", "Velocity", "Force", "Damage_Z", "Damage",

2 "Particle_ID", "Fixity", "Contact_Nodes"]

Definition at line 101 of file problem.py.

◆ THICKNESS

float problem.THICKNESS = 0.009

Definition at line 52 of file problem.py.

◆ TIP_GAP

float problem.TIP_GAP = 0.00008

Definition at line 46 of file problem.py.

◆ W

problem.W = 0.0008

Definition at line 39 of file problem.py.

◆ W_BAR

float problem.W_BAR = 1.5 * MESH_SIZE

Definition at line 48 of file problem.py.

◆ WALL_PARAMS

list problem.WALL_PARAMS

Initial value:

=  [0.021, 0.0, 0.0, 0.0,
               0.02, 0.0, 0.0, 0.0,
               0.014, -0.0015, 0.0, 0.02, 0.0015, 0.0]

Definition at line 68 of file problem.py.

◆ WALL_VY

float problem.WALL_VY = -0.06

Definition at line 46 of file problem.py.

Functions

Variables

Detailed Description

Function Documentation

◆ _nu()

◆ build_deck() [1/7]

◆ build_deck() [2/7]

◆ build_deck() [3/7]

◆ build_deck() [4/7]

◆ build_deck() [5/7]

◆ build_deck() [6/7]

◆ build_deck() [7/7]

◆ crack_speed()

◆ main()

◆ notch_void_boxes()

◆ pack_geometry()

◆ read_sites()

◆ run_with_arrival_times()

◆ seed_prenotch()

◆ validate_sites()

Variable Documentation

◆ A_IN

◆ A_OUT

◆ B_IN

◆ B_OUT

◆ BETA_N_FACTOR

◆ CONTACT_RADIUS_FACTOR

◆ CSV

◆ DENSITY

◆ DT

◆ E

◆ ELL_THETA

◆ EPSILON

◆ FINAL_TIME

◆ FRICTION_COEFF

◆ GC

◆ GC_E

◆ GC_T

◆ GRAIN_SHAPES

◆ GRAVITY

◆ H

◆ HERE

◆ HORIZON

◆ IC_VY

◆ IH

◆ IMPACT_V

◆ IMPACTOR_MASS

◆ IW

◆ K

◆ K_BULK

◆ K_E

◆ K_MAT

◆ K_T

◆ KN

◆ KN_FACTOR

◆ L_BAR

◆ MESH_DIR

◆ MESH_ELL

◆ MESH_FILE

◆ MESH_FILES

◆ MESH_SIZE

◆ MESH_TRI

◆ NCOLS

◆ NOTCH_DEPTH

◆ NOTCH_HALF

◆ NOTCH_W

◆ NROWS

◆ NU

◆ NU_E

◆ NU_T

◆ NUM_STEPS

◆ OMEGA

◆ PARTICLE_DIST

◆ PRESETS

◆ PULL_RATE

◆ QUICK

◆ R

◆ R1

◆ R_CONTACT_FACTOR

◆ R_IN