# %%
# Arc-boutement
# =============
# This example demonstrates the arc-boutement phenomenon in a simple
# two-dimensional granular configuration. A small chain of circular particles
# is compressed horizontally against a fixed particle, forming a self-locking
# structure due to frictional contacts. This effect illustrates how granular
# materials can support compressive loads through contact force chain.
# %%
# Keywords
# --------
# DEM, Rigid Body, Friction
import sys, os, shutil
import numpy as np
from migflow import scontact, time_integration

# %%
# Output Directory
# ----------------
# The output directory is prepared for storing the simulation results.
outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)


# %%
# Particle Problem
# ----------------
# We begin by defining the 2D particle problem and configuring the contact
# solver. The contact geometry is updated dynamically, and the prediction
# basis is activated to improve convergence of the contact iterations.
p = scontact.ParticleProblem(2)
p.set_fixed_contact_geometry(0)
p.set_predict_basis(1)

# %%
# Material Properties and Parameters
# ----------------------------------
# The particle density, radius, gravity vector and the friction coefficient are defined.
# A large friction coefficient is assigned between particles to ensure interlocking and prevent sliding.

rho = 1000
r = 0.01
g = np.array([0.0, -9.81])
friction_coefficient = 5
p.set_friction_coefficient(friction_coefficient, "particle", "particle")
p.set_friction_coefficient(0.0, "particle", "wall")
# %%
# Body Definition
# ---------------
# A small chain of five particles is created and initialized as a rigid body.
# The body is placed slightly above the ground, with each particle in contact
# with its neighbors. The high inverse inertia value ensures the body behaves
# as a nearly rigid object.
n = 5  # number of particles in the chain
xp = np.array([0, 3 * r])  # body center position
x_rel = np.arange(n) * 2 * r  # relative particle spacing
x_rel = np.vstack((x_rel, np.zeros_like(x_rel))).T  # particle coordinates
radii = np.full(x_rel.shape[0], r)  # uniform radii
masses = np.full(x_rel.shape[0], np.pi * rho * r**2)  # uniform masses
body = p.init_body(xp, x_rel, masses, radii, material="particle")

# %%
# Fixed Particle
# --------------
# A fixed particle is added to act as a support against which the body
# will push. This setup allows the arc-boutement effect to develop between
# the free body and the fixed particle through frictional contact.
fixed_body = p.add_body(np.array([0.0, r]), 0.0, 0.0)
p.add_particle_to_body(np.array([0.0, 0.0]), r, fixed_body, material="particle")

# %%
# Boundary Conditions
# -------------------
# A horizontal wall is added below the particles to prevent vertical motion.
p.add_boundary_segment(
    np.array([-10 * r, 0.0]), np.array([50 * r, 0.0]), tag="wall", material="wall"
)

# %%
# Time Integration
# ----------------
# The simulation runs for a short time interval with a fixed time step.
# At each iteration, the gravitational forces are applied to the particles,
# and the system evolves according to the contact dynamics solver.
t = 0.0
i = 0
dt = 1e-3
tEnd = 250 * dt
outf = 10

while t < tEnd:
    print(f"{i =:4d}, {t:.6g}/{tEnd:.6g}")
    if (i % outf) == 0:
        p.write_mig(outputdir, t)
    forces = np.pi * p.r() ** 2 * rho * g
    time_integration._advance_particles(p, forces, dt, 1, 1e-6 * r)
    t += dt
    i += 1

# %%
# Plot
# ----
# .. code-block:: shell
#
#  python3 -m migflow.plot.migplot output --actors particles
# %%
# Artifacts
# ---------
# - output/animation.mp4