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

3D Mill
=======
This example demonstrates the workings of a 3-dimensional rolling drum.
Spherical particles are added to a rotating drum, demonstratin.

Keywords
--------
DEM, 3D, Rotating Drum, Granular Mixing

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

Output Directory
----------------
The output directory is prepared for storing the simulation results.

.. code-block:: python
  
  outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
  shutil.rmtree(outputdir, ignore_errors=True)
  os.makedirs(outputdir)
  meshfile = "3d_mesh_mill"
  subprocess.call(
      ["gmsh", "-3", f"{meshfile}.geo", "-o", meshfile + ".msh", "-clscale", "1"]
  )
  

Geometrical Parameters
----------------------
We define the area in which the particles will be created inside the mill.

.. code-block:: python
  
  rout = 0.5  # Outer radius of the mill
  h = 0.6  # Height
  w = 1  # Width
  d = 0.05  # Depth
  

Material Properties and Parameters
----------------------------------

.. code-block:: python
  
  rho = 2500  # Density
  r = 9e-3  # Particle radius
  g = np.array([0.0, -9.81, 0.0])  # Gravity vector
  friction_coefficient_particle_particle = 0.3  # Friction coefficient between particles
  friction_coefficient_particle_wall = (
      0.6  # Friction coefficient between particles and wall
  )
  

Particle Problem
----------------
We begin by defining the 3D particle problem.
A friction coefficient is set for both particle-particle and particle-wall interactions based on material types.

.. code-block:: python
  
  p = scontact.ParticleProblem(3)
  p.set_friction_coefficient(friction_coefficient_particle_particle, "Glass", "Glass")
  p.set_friction_coefficient(friction_coefficient_particle_wall, "Glass", "PVC")
  
  

Particles generation
--------------------
We generate the particles on a rectangular grid inside the mill

.. code-block:: python
  
  def gen_rect(origin, w, h, d, step):
      """Generate coodinates on a rectangular grid
  
      Args:
          origin (list[float, float, float]): origin of top left corner
          w (float): width
          h (float): height
          d (float): depth
          step (float): maximal radius
  
      Returns:
          Positions where to place the
      """
  
      eps = 1e-8
      x = np.arange(-w / 2 + step - eps, w / 2 - step + eps, 2 * step) + origin[0]
      y = np.arange(step, h - step, 2 * step) + origin[1]
      z = np.arange(step, d - step, 2 * step) + origin[2]
      x, y, z = np.meshgrid(x, y, z)
      return x.reshape(-1), y.reshape(-1), z.reshape(-1)
  
  
  origin = [0, -rout, 0]
  x, y, z = gen_rect(origin, w, h, d, 1.3 * r)
  for xi, yi, zi in zip(x, y, z):
      if xi**2 + yi**2 < (rout - r) ** 2:
          p.add_particle((xi, yi, zi), r, 4 * r**3 * np.pi * rho / 3, "Glass")
  

Boundary Conditions
-------------------
Load the cylindrical mill mesh as an additional body

.. code-block:: python
  
  bnd_body = p.add_body((0, 0, 0), invertM=0, invertI=0)
  p.load_mesh_to_body(
      bnd_body, f"{meshfile}.msh", ["Outer", "Front", "Rear"], material="PVC"
  )
  

Mill Rotation Parameters

.. code-block:: python
  
  rpm = 30  # revolutions per minute
  omega = rpm * 2 * np.pi / 60  # rad/s
  

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.

.. code-block:: python
  
  p.body_omega()[bnd_body, 2] = omega
  
  t = 0.0
  i = 0
  dt = 2e-4
  tEnd = 1
  outf = 50
  mass = 4 * np.pi * p.r() ** 3 * rho / 3
  forces = mass * g
  while t < tEnd:
      print(f"{i = :5d}, {t:.6g}/{tEnd:.6g}")
      if (i % outf) == 0:
          p.write_mig(outputdir, t)
      time_integration._advance_particles(p, forces, dt, 1, 1e-4 * r)
      t += dt
      i += 1
