Download :download:`this example <../../examples/shear.py>`.

.. code-block:: python
  
  # MigFlow - Copyright (C) <2010-2026>
  # <Universite catholique de Louvain (UCL), Belgium
  #  Universite de Montpellier, France>
  #
  # List of the contributors to the development of MigFlow: see AUTHORS file.
  # Description and complete License: see LICENSE file.
  #
  # This program (MigFlow) is free software:
  # you can redistribute it and/or modify it under the terms of the GNU Lesser General
  # Public License as published by the Free Software Foundation, either version
  # 3 of the License, or (at your option) any later version.
  #
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  # GNU Lesser General Public License for more details.
  #
  # You should have received a copy of the GNU Lesser General Public License
  # along with this program (see COPYING and COPYING.LESSER files).  If not,
  # see <http://www.gnu.org/licenses/>.
  

Bidimensional shear of a dry granular material
==============================================

.. code-block:: python
  
  import os
  import time
  import gmsh
  import shutil
  import sys
  import numpy as np
  from migflow import scontact
  from migflow import time_integration
  
  np.random.seed(0)

First, let's create the output directory.
Define output directory

.. code-block:: python
  
  outputdir = "output-shear"
  shutil.rmtree(outputdir,ignore_errors=True)
  if not os.path.isdir(outputdir) :
      os.makedirs(outputdir)

Particle Problem
----------------
The particle problem is created and its dimension is provided.

.. code-block:: python
  
  p = scontact.ParticleProblem(2)
  p.set_fixed_contact_geometry(0)

The domain geometry is defined.

.. code-block:: python
  
  h = 0.1           # domain height
  w = 0.2           # domain width

The physical parameters are defined.

.. code-block:: python
  
  g =  np.array([0.,0.])                              # gravity
  rhop = 1000                                         # particle density
  r = 1e-3                                            # particle radius
  gamma = 100                                         # shear rate
  Pp = 0.05                                           # confining pressure
  eps = 1e-4
  lx = w/2-r                                          # particle area width
  ly = h/2-r-eps                                      # particle area height
  cmax=0.2                                            # compacity of the initial particle area

Top and bottom boundaries are added to the particle problem.

.. code-block:: python
  
  x0 = np.array([-w/2-w/10, 0])
  x1 = np.array([ w/2+w/10, 0])
  x_segment = np.row_stack([x0, x1])
  top_body = p.add_entities(xbody=np.array([0, h/2]), invertM=1/Pp, segments_coordinates=x_segment[None,:,:], segments_tags="Top", segments_materials="wall")
  bottom_body = p.add_entities(xbody=np.array([0, -h/2]), segments_coordinates=x_segment[None,:,:], segments_tags="Bottom", segments_materials="wall")

Periodic left and right boundaries are added.

.. code-block:: python
  
  trans = np.array([w,0])
  entity0 = p.add_boundary_periodic_entity(1, trans)
  entity1 = p.add_boundary_periodic_entity(1,-trans)
  p.add_boundary_periodic_segment((-w/2, h/2), (-w/2,-h/2), entity0)
  p.add_boundary_periodic_segment(( w/2, h/2), ( w/2,-h/2), entity1)

The confining pressure, the shear rate and the friction coefficients are prescribed

.. code-block:: python
  
  p.body_velocity()[top_body, 0] = gamma*h 
  p.set_friction_coefficient(0.5, "wall", "sand")
  p.set_friction_coefficient(0.1 , "sand", "sand")

The particles are added at random positions in the given area with the given compacity.

.. code-block:: python
  
  c=0
  pbodies = np.asarray([], np.int32)
  while c < cmax:
    xp = np.random.uniform(np.array([-lx, -ly]), np.array([lx, ly]))
    if p.n_particles() == 0:
      d = 4*r
    else:
      d = np.min(np.linalg.norm(p.position()-xp[None,:],axis=1))
    if d > 2*(r+eps):
      ri = np.random.uniform(0.8,1.0)*r
      body = p.add_particle(xp, ri, ri**2 * np.pi * rhop,"sand")
      pbodies = np.append(pbodies, body)
    c = np.sum(p.volume())/(h*w)
    print(f"Initialisation %1.1f"%(100*c/cmax) + r'%' + " ready", end="\r")
  print('')

Time integration
----------------
The numerical parameters are defined given and the initial conditions are written in the output directory.

.. code-block:: python
  
  outf = 25                                           # number of iterations between output files
  dt   = 2.5e-5                                       # time step
  tEnd =  0.5                                         # final time
  t    =  0                                           # initial time
  ii   =  0                                           # initial iteration
  p.write_mig(outputdir, 0)

Useful functions. 
The first one moves particles that cross the periodic boundaries accordingly, and the second one allows to measure the force on boundaries.

.. code-block:: python
  
  def periodic_map(w, bodies):
      p.body_position()[bodies,0] = np.remainder(p.body_position()[bodies,0]+3*w/2,w)-w/2
  
  def accumulate(f_contacts, n_divide):
      f_contacts[0] += p.get_boundary_forces("Top")/n_divide
      f_contacts[1] += p.get_boundary_forces("Bottom")/n_divide
  
  tic = time.time()
  forces_particles = np.zeros_like(p.velocity())
  forces_body = np.zeros_like(p.body_velocity())
  forces_body[top_body, 1] = -500 
  xold = p.body_position()[0,:].copy()

The computational loop can start.
The external forces given to the particles have to be given at each time step.
A lambda function is given as after_sub_iter argument so that the forces on the boundary can be computed even if the time interval is split.

.. code-block:: python
  
  while t < tEnd:
      # map the positions of the particles which crossed the periodic boundaries
      periodic_map(w, pbodies)
      # reset the velocity and position of the Top boundary
      p.body_position()[top_body,0] = xold[0]
      p.body_velocity()[top_body, 0] = gamma*h 
      # initialise the vector of the forces on the boundaries
      f_contacts = np.zeros((2,2))
      # solve the particle problem
      time_integration._advance_particles(p,forces_particles,dt,min_nsub=1,contact_tol=1e-6,after_sub_iter=lambda n_divide: accumulate(f_contacts, n_divide), max_nsub=5, fbody=forces_body)
      t += dt
      # Output files writing
      if ii%outf == 0:
          print(f"iteration : {ii} ---- t : {t} / {tEnd}")
          p.write_mig(outputdir, t)
      ii += 
