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

2D Hexagonal Polygon Deposition
===============================
This example demonstrates the deposition of hexagonal polygon-shaped
particles in a rectangular box. It uses the oriented face contact algorithm
with a fine contact detection distance to handle hexagonal grain interactions.

Keywords
--------
DEM, Polygons, Hexagonal, Deposition

.. code-block:: python
  
  import os
  import shutil
  
  import numpy as np
  from migflow import scontact, time_integration
  

Initialization
--------------

.. code-block:: python
  
  
  np.random.seed(0)
  
  outputdir = "output"
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir, exist_ok=True)
  

Physical parameters
-------------------

.. code-block:: python
  
  
  w = 1.0
  r = 1e-2 * w
  tol = 1e-4 * r
  
  mass = 1.0
  inertia = (mass * w**2) / 2.0
  
  mu = 0.4
  g = np.array([0.0, -9.81])
  

DEM problem
-----------

.. code-block:: python
  
  
  dem = scontact.ParticleProblem(2)
  dem.set_oriented_faces(True)
  dem.set_friction_coefficient(mu)
  dem.set_separation_threshold(r / 2)
  dem.contact_detection_d = 3 * r
  

Geometry utilities
------------------

.. code-block:: python
  
  
  
  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
  

Bodies
------

.. code-block:: python
  
  
  hexa = dem.add_body((0, 0), 1 / mass, 1 / inertia)
  
  x_poly, edges_poly = add_polygon((0, 0), w, w, 12, 0)
  x_smoothed, circle_centers = crop_polygon(x_poly, edges_poly, r)
  
  for center in circle_centers:
      dem.add_particle_to_body(center, r, hexa)
  
  for i in range(len(edges_poly)):
      dem.add_segment_to_body(
          x_smoothed[(2 * i + 1) % len(x_smoothed)],
          x_smoothed[(2 * i + 2) % len(x_smoothed)],
          hexa,
          "line",
      )
  
  y_min = np.min(x_smoothed[:, 1])
  
  octo = dem.add_body((0, -2 * y_min + w), 1 / mass, 1 / inertia)
  
  x_poly, edges_poly = add_polygon((0, 0), w, w, 64, 0)
  x_smoothed, circle_centers = crop_polygon(x_poly, edges_poly, r)
  
  for center in circle_centers:
      dem.add_particle_to_body(center, r, octo)
  
  for i in range(len(edges_poly)):
      dem.add_segment_to_body(
          x_smoothed[(2 * i + 1) % len(x_smoothed)],
          x_smoothed[(2 * i + 2) % len(x_smoothed)],
          octo,
          "line",
      )
  

Boundaries
----------

.. code-block:: python
  
  
  l = 10 * w
  bnd = dem.add_body((0, 0), 0, 0)
  
  dem.add_segment_to_body(np.array([+l, +y_min]), np.array([-l, +y_min]), bnd, "line")
  dem.add_segment_to_body(
      np.array([+l, +4 * y_min]), np.array([-l, +4 * y_min]), bnd, "line"
  )
  dem.add_segment_to_body(
      np.array([+l, -2 * y_min]), np.array([-l, -2 * y_min]), bnd, "line"
  )
  

Disk
----

.. code-block:: python
  
  
  disk = dem.add_body((0, 4 * y_min + w), 1 / mass, 1 / inertia)
  dem.add_particle_to_body((0, 0), w, disk)
  

Initial conditions
------------------

.. code-block:: python
  
  
  nonzero = dem.body_invert_mass().reshape(-1) != 0
  dem.body_velocity()[nonzero] = np.array([1.0, 0.0])
  
  dem.set_friction_coefficient(mu)
  

Time integration
----------------

.. code-block:: python
  
  
  i = 0
  t = 0.0
  dt = 1e-1
  t_end = 5.0
  outf = 1
  
  particles_forces = np.zeros_like(dem.velocity())
  body_forces = np.zeros_like(dem.body_velocity())
  body_forces[nonzero] = mass * g
  
  while t < t_end:
      print(f"Writing iteration {i} at time {t:.4f}")
  
      if i % outf == 0:
          dem.write_mig(outputdir, t)
  
      time_integration._advance_particles(
          dem,
          particles_forces,
          dt,
          1,
          tol,
          fbody=body_forces,
      )
  
      i += 1
      t += dt
