Download :download:`this testcase <2d_depot.py>`.

2D Polygon Particle Deposition
==============================
This example demonstrates the deposition of hexagonal polygon-shaped
particles in a rectangular box. The particles interact through oriented
contact faces, simulating the packing behavior of non-spherical granular
materials.

Keywords
--------
DEM, Polygons, Deposition, Non-spherical Particles

.. code-block:: python
  
  from migflow import scontact, time_integration
  import numpy as np
  import shutil, os, sys
  
  
  import time
  
  tic = time.time()
  
  n_particles = 100
  use_polygon = True
  # use_polygon = False
  # Seed for the rng
  np.random.seed(0)
  

output directory
----------------

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

Physical parameters of the particles
------------------------------------

.. code-block:: python
  
  rho = 2650  # [kg.m**-3]
  friction_coefficient_particle_particle = 0.4
  friction_coefficient_particle_wall = 1.4
  
  mm = 1e-3  # [m] reference length
  d_min = 0.6 * mm  # [m]
  r_min = d_min / 2
  d_max = 2 * d_min
  n_vertices = 6
  r_roundness = 0.05 * d_min
  r_roundness = 0.1 * d_min
  tol = 1e-3 * 0.05 * d_min  # [m] tolerance for the contact solver
  
  # Domain size in 2d of the cut of the cylinder
  domain_length = 50 * mm  # [m]
  domain_length = 25 * mm  # [m]
  domain_height = 700 * mm  # [m]
  x_shift = domain_length
  x_bounds = np.array([-domain_length + x_shift, domain_length + x_shift])
  

Problem definition
------------------

.. code-block:: python
  
  dem = scontact.ParticleProblem(2)
  dem.set_oriented_faces(True)
  dem.set_predict_basis(False)
  dem.set_separation_threshold(d_min / 3)
  # dem.contact_detection_d = d_min
  dem.contact_detection_d = d_min / 100
  
  # Geometry = an immovable body
  floor = dem.add_body(
      (0, 0), 0, 0
  )  # [0, 0] = origin of the reference coordinates of the body
  wall_left = dem.add_body((0, 0), 0, 0)
  wall_right = dem.add_body((0, 0), 0, 0)
  
  dem.add_segment_to_body([x_bounds[1], 0], [x_bounds[0], 0], floor, "Wall")
  dem.add_segment_to_body(
      [-domain_length / 2 + x_shift, 0],
      [-domain_length / 2 + x_shift, domain_height],
      wall_left,
      "Wall",
  )
  dem.add_segment_to_body(
      [+domain_length / 2 + x_shift, domain_height],
      [+domain_length / 2 + x_shift, 0],
      wall_right,
      "Wall",
  )
  

Particle initialisation
-----------------------

.. code-block:: python
  
  dem_particles = np.asarray([], int)
  x0 = np.array([-domain_length / 2 + x_shift, 0.8 * domain_height])
  x1 = np.array([-domain_length / 2 + x_shift, 0])
  x2 = np.array([domain_length / 2 + x_shift, 0])
  x3 = np.array([domain_length / 2 + x_shift, 0.8 * domain_height])
  bnd = np.array([x0, x1, x2, x3])
  radii = np.random.uniform(d_min / 2, d_max / 2, n_particles)
  x, r = scontact.fastdeposit_2d(bnd, radii)
  x[:, 1] += 2 * d_max
  r = 0.9 * r
  
  
  def add_polygon(x_center, w, h, n_vertices, theta_0=0.0):
      if n_vertices < 3:
          raise ValueError("n_vertices must be at least 3")
      theta = np.linspace(0, 2 * np.pi, n_vertices, endpoint=False) + theta_0
      x = w * np.cos(theta) + x_center[0]
      y = h * np.sin(theta) + x_center[1]
      vertices = np.array([x, y]).T
      nodes = np.arange(n_vertices)
      edges_0 = nodes.copy()
      edges_1 = np.append(nodes[1:], nodes[0])
      edges = np.array((edges_0, edges_1)).T
      return vertices, edges
  
  
  def crop_polygon(vertices, edges, r):
      n_vertices = len(edges) * 2
      n_centers = len(vertices)
      new_vertices = np.zeros((n_vertices, 2))
      circle_centers = np.zeros((n_centers, 2))
      for i in range(len(edges)):
          v0 = vertices[edges[i][0]]
          v1 = vertices[edges[i][1]]
          v2 = vertices[edges[(i + 1) % len(edges)][1]]
          du = (v0 - v1) / np.linalg.norm(v0 - v1)
          dv = (v2 - v1) / np.linalg.norm(v2 - v1)
          cos_theta = np.dot(du, dv)
          alpha = np.arccos(cos_theta)
          tan_theta = np.tan(alpha / 2)
          a = r / tan_theta
          x0 = v1 + a * du
          x1 = v1 + a * dv
          t_dv = np.array([-dv[1], dv[0]])
          x_center = x1 + r * t_dv
          new_vertices[2 * i + 0] = x0
          new_vertices[2 * i + 1] = x1
          circle_centers[i] = x_center
      return new_vertices, circle_centers
  
  
  x_ref, edges_ref = add_polygon((0, 0), r_min, r_min, n_vertices, 0)
  x_smoothed, circle_centers = crop_polygon(x_ref, edges_ref, r_roundness)
  for xi, ri in zip(x, r):
      polygon_mass = np.pi * ri**2 * rho  # [kg] mass of the particle
      polygon_inertia = 0.5 * polygon_mass * ri**2
      body = dem.add_body((xi[0], xi[1]), 1 / polygon_mass, 1 / polygon_inertia)
      if use_polygon:
          for center in circle_centers:
              dem.add_particle_to_body(
                  center / r_min * ri, r_roundness / r_min * ri, body
              )
          for i in range(len(edges_ref)):
              x0 = x_smoothed[(2 * i + 1) % len(x_smoothed)] / r_min * ri
              x1 = x_smoothed[(2 * i + 2) % len(x_smoothed)] / r_min * ri
              dem.add_segment_to_body(x0, x1, body, "line")
      else:
          dem.add_particle_to_body((0, 0), ri, body)
  

Simulation parameters
---------------------

.. code-block:: python
  
  dt = 1e-4
  tEnd = 1000 * dt
  outf = 10
  iit = 0
  t = 0
  
  particles_forces = np.zeros_like(dem.velocity())
  body_forces = np.divide(
      -9.81, dem.body_invert_mass(), where=dem.body_invert_mass() > 0
  ).reshape(-1, 1) * np.array([0, 1])
  nonzero = dem.body_invert_mass().reshape(-1) != 0
  dem.body_theta()[nonzero, :] = np.random.uniform(
      0, 2 * np.pi, (dem.body_theta()[nonzero, :].shape)
  )
  
  while t < tEnd:
      print(
          f"{iit=: 04d} -- t/tEnd: {t:.3f}/{tEnd: .3f} -- {t/tEnd*100:02.2f}% completed"
      )
      if iit % outf == 0:
          dem.write_mig(outputdir, t)
      dem.iterate(dt, particles_forces, body_forces, tol, 1)
      contacts = dem.contacts()
      reaction = contacts["reaction"]
      active_contacts = np.sum(np.linalg.norm(reaction, axis=1) > 0)
      print(f"  Active contacts: {active_contacts} / {len(reaction)}")
      # time_integration._advance_particles(dem, particles_forces, dt, 1, tol, fbody = body_forces)
      t += dt
      iit += 1
      dem.contacts()[:] = 0.0  # reset contact reactions for next iteration
  
  toc = time.time()
  print(f"Simulation time: {toc - tic:.2f} seconds")
  

Plot
----
.. code-block:: shell

 python3 -m migflow.plot.migplot output --actors particles

