Download :download:`this testcase <3d_simple_pills.py>`.

Simple Pills
============
This example demonstrates the simulation of capsule-shaped particles (pills)
within a cylindrical container. The pills are randomly oriented and packed
inside the drum, which is later removed to observe the settling behavior of
the pills under gravity. The simulation highlights the complex interactions
and dynamics of non-spherical granular materials.

Keywords
--------
DEM, 3D, Pills, Non-spherical Particles

.. code-block:: python
  
  import numpy as np
  import os, shutil, sys
  from migflow import scontact, time_integration
  from migflow.quaternion import (
      Quaternion,
      random_rotation,
      to_quaternion,
      quaternion_between_vectors,
  )
  from cylinder_mod import generate_geometry
  

Output Directory
----------------
Prepare the directory that will store the transient simulation results.

.. code-block:: python
  
  outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
  shutil.rmtree(outputdir, True)
  
  

Pill Geometry Utilities
-----------------------
Helper routines to assemble capsule-shaped bodies and distribute them in the drum.

.. code-block:: python
  
  def add_pill_to_problem(
      _p: scontact.ParticleProblem,
      R: float,
      AR: float,
      density: float,
      position: np.ndarray,
      direction: np.ndarray,
  ):
      """Create a capsule made of a cylinder and two hemispherical caps, then attach it to the problem."""
  
      # Pills default generation axis
      default_axis = np.array([0, 0, 1])
  
      # cylinder
      cyl_H = R * (AR - 2)
      cyl_mass = cyl_H * np.pi * R**2 * density
      cyl_inertia = (
          0.5
          * cyl_mass
          * R**2
          * np.array([0.5 + (cyl_H / R) ** 2 / 6, 0.5 + (cyl_H / R) ** 2 / 6, 1])
          * np.eye(3)
      )
  
      # semi spherical cap
      rel_pos = np.array([[0, 0, cyl_H / 2], [0, 0, -cyl_H / 2]])
      cap_mass = 2 / 3 * np.pi * R**3 * density
      cap_shift = cyl_H / 2 + 3 / 8 * R
      cap_inertia = cap_mass * R**2 * np.array([83 / 320, 83 / 320, 2 / 5]) * np.eye(3)
      cap_inertia_total = cap_inertia + cap_mass * (
          cap_shift**2 * np.eye(3)
          - np.outer(np.array([0, 0, cap_shift]), np.array([0, 0, cap_shift]))
      )
      cap_inertia_total += cap_inertia + cap_mass * (
          cap_shift**2 * np.eye(3)
          - np.outer(np.array([0, 0, -cap_shift]), np.array([0, 0, -cap_shift]))
      )
  
      inertia = cyl_inertia + cap_inertia_total
      mass = cyl_mass + 2 * cap_mass
  
      # Pill orientation
      qrot = quaternion_between_vectors(default_axis, direction)
      rot_pos = np.vstack([qrot.rotate_3Dvec(rel_pos[0]), qrot.rotate_3Dvec(rel_pos[1])])
      rotmat = qrot.to_rot_mat()
      rotmat[np.abs(rotmat) < 1e-14] = 0
      inertia_rot = rotmat @ inertia @ rotmat.T
      eigenvalues, eigenvectors = np.linalg.eigh(inertia_rot)
      if np.linalg.det(eigenvectors) < 0:
          eigenvectors[:, 0] = -eigenvectors[:, 0]
      qtheta = Quaternion.from_basis(eigenvectors)
      norm_theta = np.linalg.norm(qtheta.to_array())
      if not np.allclose(norm_theta, 1.0):
          print("Quaternion norm is too small", norm_theta)
  
      # Body initiation
      body = _p.add_body(position, 1 / mass, 0)
      p.body_invert_inertia()[body] = 1 / eigenvalues
      # Adding spheres and cylinder to body
      _p.add_segment_to_body(
          rot_pos[0], rot_pos[1], body, tag="bnd", material="default", radius=R
      )
      _p.add_particle_to_body(rot_pos[0], R, body, material="default")
      _p.add_particle_to_body(rot_pos[1], R, body, material="default")
      p.body_theta()[body] = to_quaternion(direction)
      return body
  
  
  def gen_mesh_geometry(cyl_radius, cyl_height):
      """Build or reuse a Gmsh mesh describing the confining cylinder."""
      print(">> Trying generating mesh...")
      meshName = "mesh.msh"
      # meshPath = os.path.join(os.getcwd(), meshName)
      meshPath = meshName
      if os.path.exists(meshPath):
          print(f" -> Mesh already exists at {meshPath} \(ÔÔ)/")
          return meshPath
      else:
          meshPath = generate_geometry(cyl_radius, cyl_height, meshPath)
          print(f" -> Mesh created at {meshPath} \(ÔÔ)/")
      return meshPath
  
  
  def gen_random_layer(pillR, cylR):
      """Generate a non-overlapping set of 2D positions inside the drum cross section."""
      d = pillR
      polypod_r = d / 2
      itmax = 2000
      it = 0
      positions = []
      while it < itmax:
          it += 1
          pos = np.array([np.random.uniform(-cylR, cylR), np.random.uniform(-cylR, cylR)])
          if np.linalg.norm(pos) < cylR - 1.05 * polypod_r:
              valid = True
              for p in positions:
                  if np.linalg.norm(pos - p) < d:
                      valid = False
                      break
              if valid:
                  positions.append(pos)
      positions = np.array(positions)
      return positions
  
  
  def add_pill_layer(p, insert_height, pill_r, pill_ar, rho, cylr):
      """Insert a horizontal layer of pills with random orientations."""
      positions = gen_random_layer(pill_r * pill_ar, cylr)
      hh = np.ones_like(positions[:, 0]) * insert_height
      positions = np.hstack((positions, hh[:, None]))
      for pos in positions:
          bid = add_pill_to_problem(p, pill_r, pill_ar, rho, pos, random_rotation())
          # p.body_omega()[bid] = np.random.uniform(-1,1,3) * 2*np.pi
          p.body_velocity()[bid] = (
              np.array(
                  [
                      np.random.uniform(-1, 1),
                      np.random.uniform(-1, 1),
                      np.random.uniform(-0.5, -1),
                  ]
              )
              * pill_r
              * pill_ar
          )
      return positions.shape[0]
  
  

Particle Problem
----------------
We begin by defining the 3D 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.

.. code-block:: python
  
  p = scontact.ParticleProblem(3)
  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.

.. code-block:: python
  
  rho = 1000
  pill_r = 0.05
  pill_ar = 4.0
  pill_R = pill_r * pill_ar
  g = np.array([0.0, 0.0, -9.81])
  
  H = 1.0
  R = 1.0
  
  material_wall = "wall"
  material_bottom = "bottom"
  mu_pill_pill = 0.4
  mu_pill_wall = 0.3
  mu_pill_bottom = 0.5
  p.set_friction_coefficient(mu_pill_pill, "default", "default")
  p.set_friction_coefficient(mu_pill_wall, "default", material_wall)
  p.set_friction_coefficient(mu_pill_bottom, "default", material_bottom)
  

Pill Packing
------------
Fill the container with randomly oriented pills and report the number of created bodies.

.. code-block:: python
  
  add_pill_layer(p, pill_r * pill_ar, pill_r, pill_ar, rho, 0.5 * R)
  print(">> Done %d pills generation" % p.n_bodies())
  

Mesh And Boundaries
-------------------
Generate the confining mesh, attach it as a rigid body, and mark the boundary patches.

.. code-block:: python
  
  mesh = generate_geometry(R, H)
  p.load_msh_boundaries(mesh, tags=["bottom_inner"], material=material_bottom)
  
  # Write initial configuration
  p.write_mig(outputdir, 0, write_contacts=True)
  

Body Forces
-----------
Pre-compute gravitational force contributions for each rigid body.

.. code-block:: python
  
  fp = np.zeros((p.n_particles(), 3))
  inv_mass = p.body_invert_mass().ravel()
  mass = np.where(inv_mass > 0, 1.0 / inv_mass, 0.0)
  bf = mass[:, None] * g
  pills_weight = bf.copy()
  

Time Integration
Select constant time step, solver tolerance, and output frequency.
Advance the system in time, periodically exporting states and removing the drum once settled.

.. code-block:: python
  
  tEnd = 5.0
  dt = 1e-3
  ctol = pill_r / 100
  outf = int(0.01 / dt)
  t = 0
  i = 0
  
  while t < tEnd:
      sub = time_integration._advance_particles(p, fp, dt, 1, ctol, fbody=bf)
      q, _ = p.get_current_quality(ctol)
  
      fbottom = p.get_boundary_forces("bottom_inner")[2]
      qratio = abs(q) / ctol / sub
      if qratio < 1:
          print(
              f"Quality ratio : {abs(q)/ctol/sub:.6f}[-], NSUB {sub:d} Bottom reaction force: {fbottom:.6g} N, Pills weight: {np.abs(pills_weight[:,2].sum()):.6g} N"
          )
      else:
          print(f"Error, solver tolerance is not respected")
  
      t += dt
      i += 1
  
      if i % outf == 0:
          print(f"{i =:4d}, {t:.6g}/{tEnd:.6g}")
          p.write_mig(outputdir, t, write_contacts=True)
