tutorial/two_particles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e70892dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import sys\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "07b8fa92",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-a4d10d561234>:17: MatplotlibDeprecationWarning: Support for setting the 'text.latex.preamble' or 'pgf.preamble' rcParam to a list of strings is deprecated since 3.3 and will be removed two minor releases later; set it to a single string instead.\n",
      "  plt.rcParams['text.latex.preamble']=[r\"\\usepackage{amsmath}\"]\n"
     ]
    }
   ],
   "source": [
    "# for plotting\n",
    "import matplotlib\n",
    "# matplotlib.use('TkAgg') \n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set()\n",
    "# sns.set_style(\"whitegrid\")\n",
    "sns.set(font_scale=1.5, rc={'text.usetex' : True})\n",
    "sns.set_context(\"paper\", font_scale=4, rc={\"lines.linewidth\": 4})\n",
    "mpl.rcParams['lines.linewidth'] = 3\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "plt.rcParams['text.latex.preamble']=[r\"\\usepackage{amsmath}\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d0d9dc7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import util functions\n",
    "from util import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d27b4868",
   "metadata": {},
   "source": [
    "# 1 Two-particle test setup\n",
    "Particle $p_2$ is falling freely due to gravity $g = [0, -10, 0]$ in units of m/s$^2$ (i.e., gravity acts along downward direction). The initial velocity of $p_2$ is zero and the initial distance (boundary-to-boundary) is $d$ m, see Figure 1(a). In Figure 1(b), we show the equivalent setup where we start $p_2$ much closer to the $p_1$, with $\\delta$ satisfying $0 \\leq \\delta < d$, and assign $p_2$ the initial velocity $v_0 = \\sqrt{2 |g| \\delta}$ that it would have acquired if it were falling freely $\\delta$ distance due to the influence of gravity $g$.   \n",
    "\n",
    "| <img src=\"assets/two-particles.png\" width=\"400\"> | \n",
    "| :---: | \n",
    "| **Figure 1.** Setup of a two-particle test |\n",
    "\n",
    "To simulate this in PeriDEM, we need to do following:\n",
    "- describe model setting such as\n",
    "  * final time and time step\n",
    "  * geometry of particles along with $d$ and $\\delta$\n",
    "  * discretization of two particles \n",
    "  * horizon (nonlocal length scale) in peridynamics\n",
    "- material properties of two particles\n",
    "- contact properties of two particles\n",
    "- boundary conditions\n",
    "\n",
    "Before we list the set of parameters, there is one key idea that we need to discuss. In the PeriDEM library, we allow arbitrary grouping of particles and walls. We refer to such groups by `zone`. \n",
    "\n",
    "## 1.1 Grouping particles - Zones\n",
    "Suppose we are interested in multi-particle test. Consider following scenarios:\n",
    "- **Case 1.** Suppose all particles have same shape. Suppose first $m$ particles, denoted by set $S_1 = \\{1, 2, ..., m\\}$, out of total $N$ particles have material properties of steel $M_1$ and remaining particles, in set $S_2 = \\{m+1, m+2, ..., N\\}$, of glass $M_2$. Suppose for all pair of particles, contact parameters are same.\n",
    "- **Case 2.** Consider **case 1** but now also assume that\n",
    "  * pair of particles in $S_1$ have contact properties $C_{11}$\n",
    "  * pair of particles in $S_2$ have contact properties $C_{22}$\n",
    "  * remaining pair of particles have contact properties $C_{12}$\n",
    "- **Case 3.** Start from scratch and suppose now first $r$ particles, in set $\\bar{S}_1 = \\{1, 2, ..., r\\}$, are of circular shape, and remaining particles, in set $\\bar{S}_2 = \\{r+1, r+2, ..., N\\}$, are hexagon shaped. Also suppose all particles have same material properties and same contact propertis.\n",
    "- **Case 4.** Let **case 3** holds. We now ask that some particles of set $\\bar{S}_1$ have material properties $\\bar{M}_{11}$ and some have material properties $\\bar{M}_{12}$. Similarly, we ask that some subset of $\\bar{S}_2$ have material properties $\\bar{M}_{21}$ and remaining have $\\bar{M}_{22}$\n",
    "- **Case 5.** Let **case 1** holds. We now introduce a fixed wall to the particulate media. We assume that particles and walls have same material properties and same contact properties.\n",
    "- **Case 6.** Start from scract, and assume that all particles are of same shape, same material properties, contact is same for any pair of particles. But now we assume that first $s$ particles, in set $\\tilde{S}_1 = \\{1, 2, ..., s\\}$, are assigned a fine mesh, and remaining, in set $\\tilde{S}_2 = \\{s+1, s+2, ..., N\\}$, are assigned a coarse mesh. \n",
    "\n",
    "The reason for outlining some of the possibility above is to show that in a multi-particle system, we often have to deal of subgroup specified properties. When properties change from one subset to another, it automatically generates a subgrouping. \n",
    "\n",
    "### 1.1.1 Case 1\n",
    "Consider **case 1**. Here it is natural to create two zones (with zone we mean the subgroup); we will assign all particles in set $S_1$ the zone id $1$ and remaining particles zone id $2$. \n",
    "\n",
    "Once we have decided that simulation consists of two zones, we need to prescribe parameters such as\n",
    "- Reference particle shape (see Remark 1 below on reference particle)\n",
    "- Material properties\n",
    "- Mesh of a reference particle (since all other particles are affine-transform related to reference particle, we only need mesh of a reference particle!)\n",
    "- Contact defined for pair of zones\n",
    "- Initial condition and boundary conditions (force or displacement) can either be applied by listing **global id** (i.e. 0, 1, 2, ..., N-1) or **zone id**. In the latter case, the conditions are applied to all particles in the specified **zone**.\n",
    "\n",
    "> **Remark 1.** Reference particle is such that all particles are related by affine transformation of reference particle. You can think of finite element method where shape functions and their derivatives over any arbitrary element is computed using the map taking reference element to the arbitrary element. This reference element of course has a same shape as the arbitrary element.\n",
    "\n",
    "For **case 1**, input file will have two zones and the properties as said above will be zone specific. Contact properties will be specified for the pair of zones, i.e., *<zone 1, zone 1>*, *<zone 2, zone 2>*, *<zone 1, zone 2>*. Note that contact properties are same for *<zone 1, zone 2>* and *<zone 2, zone 1>*. \n",
    "\n",
    "> **Remark 2.** When there is only one group (in situation when particles are of same shape, same material, same contact properties, and are obtained from same mesh of a reference particle), we will have just one **zone** in the simulation and all particles will belong to this zone.\n",
    "\n",
    "### 1.1.2 Case 2\n",
    "In this case, having just two zones (repeating that **zone** simply means subgroup) is enough as the contact in our simulation is always assigned for a pair of zones. So we can have the same input file as in **case 1** and simply modify the contact block of the input file where we can have properties $C_{11}, C_{22}, C_{12}$ for *<zone 1, zone 1>*, *<zone 2, zone 2>*, and *<zone 1, zone 2>* pairs.\n",
    "\n",
    "### 1.1.3 Case 3\n",
    "In this case, since the shape of particles differ and one group of particles defined by set $\\bar{S}_1$ are of circular shape and remainig particles are of hexagon shape, it makes sense to again create two **zones**:\n",
    "- add particles in set $\\bar{S}_1$ to **zone 1**\n",
    "- add remaining particles to **zone 2**.\n",
    "\n",
    "We then specify the properties, as before, in a `zone-specific` fashion. We provide the details of reference particle (i.e., circle) and associated mesh in **zone 1** area of input file and do same of hexagon reference particle and associated mesh in **zone 2** area of input file. \n",
    "\n",
    "As before, we will define material properties for all particles in zone 1 and zone 2 and define the contact properties using pair of zones. We can have different material properties for particles in zone 1 and zone 2 and similary have different contact properties for different pair of two zones.\n",
    "\n",
    "### 1.1.4 Case 4\n",
    "In this case, having just two zones is not enough. If we look at **case 3**, we see that **case 4** further subdivides the **case 3** in the basis of material properties. We need **four zones** to handle this case:\n",
    "- **zone 1.** circle shaped particles with material properties $\\bar{M}_{11}$\n",
    "- **zone 2.** circle shaped particles with material properties $\\bar{M}_{12}$\n",
    "- **zone 3.** hexagon shaped particles with material properties $\\bar{M}_{21}$\n",
    "- **zone 4.** hexagon shaped particles with material properties $\\bar{M}_{22}$.\n",
    "\n",
    "We then define the material properties, reference particle geometry and mesh in a `zone-specific` fashion and define the contact properties for pair of zones. In this case we have to specify 10 contact properties as there are 10 unique pairs of zones (we take *<i,j>* pair same as *<j,i>* for the purpose of contact).\n",
    "\n",
    "### 1.1.5 Case 5\n",
    "In our implementation, both particle and wall inherit from a base particle implementation, see [BaseParticle](https://github.com/prashjha/PeriDEM/blob/main/src/particle/baseParticle.h), however, conceptually both are treated differently. One of the major difference is as follows:\n",
    "\n",
    "**Particle have associated reference particle so that they can be obtained by affine transformation of reference particle. However, for walls we do not assume or require to have reference particle. In the simulation, one has to describe all geometric details and mesh file for each wall seperately.**\n",
    "\n",
    "Wall always belong to a different and unique zone. If we have n walls, we will atleast have n zones. If we have n walls and N particles, in the simplest case, we will have n+1 zone where n zones will correspond to n walls, and remaining zone will correspond to all the N particles. If particles are required to be grouped in smaller groups, say 4 groups, we will have n+4 zones.\n",
    "\n",
    "Now that walls are assigned unique zones such that only one wall can be in one zone, and since the properties are specified in a `zone-specific` fashion, it is easy to specify properties of walls and group of particles in different zones. \n",
    "\n",
    "For the specific example in **case 5**, we will have 3 zones:\n",
    "- **zone 1.** For the wall\n",
    "- **zone 2.** For particles in set $S_1$ (see description of **case 1**)\n",
    "- **zone 3.** For remaining particles.\n",
    "We then have properties such as mesh, reference particle (only for zone 2 and 3), material properties for each zone, and have contact properties for pair of zones (in this case we need to assign contact properties for 6 unique pair of zones). \n",
    "\n",
    "### 1.1.6 Case 6\n",
    "This example shows why concept of zone is handy. In this case, all particles are identical except their meshes. In this case it is natural to have two zones:\n",
    "- **zone 1.** For particles with fine mesh\n",
    "- **zone 2.** For particles with coarse mesh."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c8703e5",
   "metadata": {},
   "source": [
    "#  2 Preparation of input files for the two-particle test \n",
    "We now discuss the preparation of two-particle setup one by one.\n",
    "\n",
    "## 2.1 Zones\n",
    "Since we assume $p_1$ and $p_2$ have different material properties, we create two zones:\n",
    "- **zone 1.** Particle 1 belongs to this zone\n",
    "- **zone 2.** Particle 2 belongs to this zone.\n",
    "\n",
    "\n",
    "## 2.2 Model setting\n",
    "These are defined as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4c630fb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## main units are meter and seconds\n",
    "\n",
    "## set path for files created\n",
    "fpath = './two_particles/'\n",
    "Path(fpath).mkdir(parents=True, exist_ok=True)\n",
    "\n",
    "## origin (this center is used to discretize the particle and\n",
    "## later discretization is translated and rotated appropriately)\n",
    "center = [0., 0., 0.]\n",
    "\n",
    "## gravity\n",
    "g = [0., -10., 0.]\n",
    "\n",
    "## radius of p1 and p2\n",
    "R1 = 0.0015\n",
    "R2 = 0.001\n",
    "\n",
    "## mesh size (use smaller radius to decide)\n",
    "h = R1 / 8.\n",
    "if R2 < R1:\n",
    "    h = R2 / 8.\n",
    "    \n",
    "## peridynamics horizon (typically taken as 2 to 4 times the mesh size)\n",
    "horizon = 3. * h\n",
    "\n",
    "\n",
    "## initial distance between two particles\n",
    "d = 0.001\n",
    "# we take smallest delta while keeping p1 and p2 seperated by distance 'horizon'\n",
    "delta = d - horizon \n",
    "v0 = [0, -np.sqrt(2. * np.abs(g[1]) * delta), 0.]\n",
    "\n",
    "## final time and time step\n",
    "T = 0.04\n",
    "N_steps = 200000 \n",
    "# if dt is large, the simulation could produce garbage results so care is needed\n",
    "dt = T / N_steps \n",
    "# steps per output, i.e., simulation will generate output every 50 steps\n",
    "N_out = N_steps / 50 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ec59d81",
   "metadata": {},
   "source": [
    "## 2.3 Material properties\n",
    "The deformation of particles are governed by peridynamics model. To fully describe the properties of material, we need density $\\rho$, bulk modulus $K$, Poisson ratio $\\nu$, and critical energy release rate $G_c$. We define material properties of $p_1$ and $p_2$ below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "971b0e60",
   "metadata": {},
   "outputs": [],
   "source": [
    "## material\n",
    "# for zone 1 (in this case zone 1 has just one particle p1)\n",
    "nu1 = 0.25\n",
    "rho1 = 1200.\n",
    "K1 = 2.16e+7\n",
    "E1 = get_E(K1, nu1)  # see util.py\n",
    "G1 = get_G(E1, nu1)\n",
    "Gc1 = 50.\n",
    "\n",
    "# for zone 2 (in this case zone 2 has just one particle p2)\n",
    "nu2 = 0.25\n",
    "rho2 = 1200.\n",
    "K2 = 2.16e+7\n",
    "E2 = get_E(K2, nu2)\n",
    "G2 = get_G(E2, nu2)\n",
    "Gc2 = 50."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0340fe49",
   "metadata": {},
   "source": [
    "## 2.4 Contact parameters\n",
    "We now have to prescribe the contact parameters. As discussed earlier, contact is defined for pair of zones. In this case we have to assign contact properties for 3 pair of zones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1e646213",
   "metadata": {},
   "outputs": [],
   "source": [
    "## contact radius is taken as R_contact_factor * h_{min} where\n",
    "## h_{min} is the minimum mesh size over all meshes (it is computed at the begining \n",
    "## of simulation)\n",
    "## we only specify the R_contact_factor\n",
    "R_contact_factor = 0.95\n",
    "\n",
    "## define Kn parameter for normal contact force (see PeriDEM paper for more details)\n",
    "Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))\n",
    "Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))\n",
    "Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))\n",
    "\n",
    "## define beta_n parameter for normal damping\n",
    "beta_n_eps = 0.9\n",
    "# factor for damping (Constant C in the PeriDEM paper)\n",
    "beta_n_factor = 100.\n",
    "damping_active = True \n",
    "\n",
    "## friction coefficient\n",
    "friction_coeff = 0.5\n",
    "friction_active = False # friction is not active"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0faeaa50",
   "metadata": {},
   "source": [
    "## 2.5 Particle locations\n",
    "We create a `.csv` file that lists out the `zone id` (0-numbered), `radius`, `center location` (x, y, z), `orientation` (the initial confguration is obtained by rotating the particle mesh by amount specified in `orientation`). \n",
    "\n",
    "Function below does this job:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8115de54",
   "metadata": {},
   "outputs": [],
   "source": [
    "def two_particle_locations(filename, radius, centers, pp_tag = None):\n",
    "    \"\"\"\n",
    "    Creates .csv file containing the center, radius, orientation, and zone id of each particle in two-particle setup\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    filename : str\n",
    "      Filename of .csv file to be created\n",
    "    radius: list\n",
    "      List containing radius of two particles\n",
    "    centers: list of list\n",
    "      List of list containing the centers of two particles\n",
    "    pp_tag: str, optional\n",
    "      Postfix .geo file with this tag\n",
    "    \"\"\"\n",
    "    pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else ''\n",
    "    inpf = open(filename + pp_tag_str + '.csv','w')\n",
    "    inpf.write('i, x, y, z, r, o\\n')\n",
    "    for i in range(2):\n",
    "        inpf.write('%d, ' % i)\n",
    "        inpf.write(print_list(centers[i], '%Lf'))\n",
    "        inpf.write(', %Lf' % radius[i])\n",
    "        o = 0. if i == 0 else np.pi * 0.5\n",
    "        inpf.write(', %Lf\\n' % o)\n",
    "    inpf.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "363c368a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create particle_locations.csv file \n",
    "centers = []\n",
    "centers.append([R1, R1, 0.])\n",
    "centers.append([R1, 2.*R1 + R2 + d - delta, 0.])\n",
    "two_particle_locations(fpath + 'particle_locations', [R1, R2], centers)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d93de9c3",
   "metadata": {},
   "source": [
    "## 2.6 Discretization of reference particle\n",
    "We have two zones with each containing single circular particle. We need to define reference particles for both zones and provide mesh file. \n",
    "\n",
    "We let both reference particle (circle in this case) be at center $[0,0,0]$. Radius of reference particle for zone 1 is same as radius of particle $p_1$, i.e., $R_1$. Similar is true for other reference particle.\n",
    "\n",
    "We use `Gmsh` to create the mesh of particles. First step to mesh is to create a `.geo` file that `Gmsh` understands. This is done next:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "23e9756d",
   "metadata": {},
   "outputs": [],
   "source": [
    "## file below creates .geo file (see util.py)\n",
    "# we need to provide the name of file (function will add .geo suffix)\n",
    "# and also center, radius, and mesh size\n",
    "# center is the center of reference particle. In the simulation, reference particle will affinely transformed\n",
    "# to get the actual particle.\n",
    "# Note center is defined as [0., 0., 0.] somewhere above\n",
    "generate_circle_gmsh_input(fpath + 'p1', center, R1, h) # for p1\n",
    "generate_circle_gmsh_input(fpath + 'p2', center, R2, h) # for p2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e633d1d0",
   "metadata": {},
   "source": [
    "We can now discretize the circles using `.geo` file as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c8404652",
   "metadata": {},
   "outputs": [],
   "source": [
    "fp = fpath + 'p1'\n",
    "# create .msh file\n",
    "!gmsh $fp'.geo' -2 -o $fp'.msh' &> /dev/null\n",
    "# .vtk is not required. But it may allow you to check if mesh looks okay\n",
    "!gmsh $fp'.geo' -2 -o $fp'.vtk' &> /dev/null  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4e4bd6bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fp = fpath + 'p2'\n",
    "# create .msh file\n",
    "!gmsh $fp'.geo' -2 -o $fp'.msh' &> /dev/null\n",
    "# .vtk is not required. But it may allow you to check if mesh looks okay\n",
    "!gmsh $fp'.geo' -2 -o $fp'.vtk' &> /dev/null  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5d31293",
   "metadata": {},
   "source": [
    "## 2.7 Putting all together\n",
    "We now have defined and obtained all components for the simulation and we only need to create a instruction file `<your filename>.yaml`. \n",
    "\n",
    "We do this in one go below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b2cccba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## open a input file (we choose pretty obvious name input.yaml)\n",
    "inpf = open(fpath + 'input.yaml','w')\n",
    "\n",
    "## provide model specific details\n",
    "inpf.write(\"Model:\\n\")\n",
    "inpf.write(\"  Dimension: 2\\n\")\n",
    "inpf.write(\"  Discretization_Type:\\n\")\n",
    "inpf.write(\"    Spatial: finite_difference\\n\") # we refer to meshfree discretization by finite_difference\n",
    "inpf.write(\"    Time: central_difference\\n\")   # you can use either central_difference or velocity_verlet\n",
    "inpf.write(\"  Final_Time: %4.6e\\n\" % (T))\n",
    "inpf.write(\"  Time_Steps: %d\\n\" % (N_steps))\n",
    "\n",
    "## not used (I should get rid of this soon)\n",
    "inpf.write(\"Policy:\\n\")\n",
    "inpf.write(\"  Enable_PostProcessing: true\\n\")\n",
    "\n",
    "## container info\n",
    "## Container is the region in 2d (or in 3d) where you expect all activites \n",
    "## to take place. If not sure, you can take a wild guess of bounding box of\n",
    "## this simulation.\n",
    "inpf.write(\"Container:\\n\")\n",
    "inpf.write(\"  Geometry:\\n\")\n",
    "inpf.write(\"    Type: rectangle\\n\")\n",
    "# we define container to be rectangle. \n",
    "# in the code, we need to provide 6 element vector to define a rectangle\n",
    "# first 3 in the vector correspond to coordinate of the left-bottom corner\n",
    "# and remaining correspond to coordinates of the right-top corner\n",
    "contain_params = [0., 0., 0., 2.*R1, 2.*R1 + 2.*R2 + d, 0.]\n",
    "if R2 > R1:\n",
    "    contain_params[3] = 2.*R2\n",
    "inpf.write(\"    Parameters: \" + print_dbl_list(contain_params))\n",
    "\n",
    "## zone info\n",
    "## we have two zones\n",
    "inpf.write(\"Zone:\\n\")\n",
    "inpf.write(\"  Zones: 2\\n\")\n",
    "\n",
    "## zone 1 (bottom particle)\n",
    "inpf.write(\"  Zone_1:\\n\")\n",
    "# specify that this is not a zone containing wall\n",
    "inpf.write(\"    Is_Wall: false\\n\")\n",
    "\n",
    "## zone 2 (top particle)\n",
    "inpf.write(\"  Zone_2:\\n\")\n",
    "inpf.write(\"    Is_Wall: false\\n\")\n",
    "\n",
    "## particle info\n",
    "## as we said earlier, the properites are specific in a `zone-specific` manner\n",
    "## what I really mean by `zone-specific` will be clear now\n",
    "inpf.write(\"Particle:\\n\")\n",
    "# for two-particle test, we not only output the standard .vtu simulation file\n",
    "# but also some other quantities such as maximum distance between two particles\n",
    "# this is to debug the code and to account for the damping effects\n",
    "inpf.write(\"  Test_Name: two_particle\\n\")\n",
    "# specify reference particle details for zone 1 particles (in this case zone 1 has just one particle)\n",
    "inpf.write(\"  Zone_1:\\n\")\n",
    "# reference particle is of circular type\n",
    "inpf.write(\"    Type: circle\\n\")\n",
    "# in the code, circle is defined using 4 vector array\n",
    "# first element in the vector is radius, and rest are the coordinates of the center\n",
    "# recall that we fixed referene particles at center = [0,0,0]\n",
    "p1_geom = [R1, center[0], center[1], center[2]]\n",
    "inpf.write(\"    Parameters: \" + print_dbl_list(p1_geom)) \n",
    "inpf.write(\"  Zone_2:\\n\")\n",
    "# zone 2 reference particle is also circle. The radius is different and center is same.\n",
    "inpf.write(\"    Type: circle\\n\")\n",
    "p2_geom = [R2, center[0], center[1], center[2]]\n",
    "inpf.write(\"    Parameters: \" + print_dbl_list(p2_geom)) \n",
    "\n",
    "## particle generation\n",
    "## We generate particles by reading the particle_locations.csv file which provides list of particles\n",
    "## in the simulation (each row correspond to unique particle)\n",
    "## for each particle, particle_locations.csv provides information such as \n",
    "## <zone id>, <center>, <radius>, <orientation>\n",
    "inpf.write(\"Particle_Generation:\\n\")\n",
    "inpf.write(\"  From_File: particle_locations.csv\\n\")\n",
    "# loc_rad_orient basically says we provide center location, radiua, and orientation\n",
    "inpf.write(\"  File_Data_Type: loc_rad_orient\\n\") \n",
    "\n",
    "## Mesh info\n",
    "inpf.write(\"Mesh:\\n\")\n",
    "\n",
    "# zone 1\n",
    "inpf.write(\"  Zone_1:\\n\")\n",
    "inpf.write(\"    File: p1.msh \\n\")\n",
    "\n",
    "# zone 2\n",
    "inpf.write(\"  Zone_2:\\n\")\n",
    "inpf.write(\"    File: p2.msh \\n\")\n",
    "\n",
    "## Contact info\n",
    "## we have to provide contact parameters for pair of zones\n",
    "inpf.write(\"Contact:\\n\")\n",
    "\n",
    "# 11\n",
    "# this fixes the contact between two particles in zone 1\n",
    "inpf.write(\"  Zone_11:\\n\")\n",
    "inpf.write(\"    Contact_Radius_Factor: %4.6e\\n\" % (R_contact_factor))\n",
    "\n",
    "if damping_active == False:\n",
    "    inpf.write(\"    Damping_On: false\\n\")\n",
    "if friction_active == False:\n",
    "    inpf.write(\"    Friction_On: false\\n\")\n",
    "\n",
    "inpf.write(\"    Kn: %4.6e\\n\" % (Kn_11))\n",
    "inpf.write(\"    Epsilon: %4.6e\\n\" % (beta_n_eps))\n",
    "inpf.write(\"    Friction_Coeff: %4.6e\\n\" % (friction_coeff))\n",
    "inpf.write(\"    Kn_Factor: 1.0\\n\")\n",
    "inpf.write(\"    Beta_n_Factor: %4.6e\\n\" % (beta_n_factor))\n",
    "\n",
    "# 12\n",
    "# this fixes the contact between two particles in zone 1 and zone 2\n",
    "inpf.write(\"  Zone_12:\\n\")\n",
    "inpf.write(\"    Contact_Radius_Factor: %4.6e\\n\" % (R_contact_factor))\n",
    "\n",
    "if damping_active == False:\n",
    "    inpf.write(\"    Damping_On: false\\n\")\n",
    "if friction_active == False:\n",
    "    inpf.write(\"    Friction_On: false\\n\")\n",
    "\n",
    "inpf.write(\"    Kn: %4.6e\\n\" % (Kn_12))\n",
    "inpf.write(\"    Epsilon: %4.6e\\n\" % (beta_n_eps))\n",
    "inpf.write(\"    Friction_Coeff: %4.6e\\n\" % (friction_coeff))\n",
    "inpf.write(\"    Kn_Factor: 1.0\\n\")\n",
    "inpf.write(\"    Beta_n_Factor: %4.6e\\n\" % (beta_n_factor))\n",
    "\n",
    "# 22\n",
    "inpf.write(\"  Zone_22:\\n\")\n",
    "inpf.write(\"    Contact_Radius_Factor: %4.6e\\n\" % (R_contact_factor))\n",
    "\n",
    "if damping_active == False:\n",
    "    inpf.write(\"    Damping_On: false\\n\")\n",
    "if friction_active == False:\n",
    "    inpf.write(\"    Friction_On: false\\n\")\n",
    "\n",
    "inpf.write(\"    Kn: %4.6e\\n\" % (Kn_22))\n",
    "inpf.write(\"    Epsilon: %4.6e\\n\" % (beta_n_eps))\n",
    "inpf.write(\"    Friction_Coeff: %4.6e\\n\" % (friction_coeff))\n",
    "inpf.write(\"    Kn_Factor: 1.0\\n\")\n",
    "inpf.write(\"    Beta_n_Factor: %4.6e\\n\" % (beta_n_factor))\n",
    "\n",
    "## Neighbor info\n",
    "## at present, this block has no impact in the simulation. We have it for later extensions of the library.\n",
    "inpf.write(\"Neighbor:\\n\")\n",
    "inpf.write(\"  Update_Criteria: simple_all\\n\")\n",
    "inpf.write(\"  Search_Factor: 5.0\\n\")\n",
    "\n",
    "## Material info\n",
    "## material properties have to be specified in a `zone-specific` manner\n",
    "inpf.write(\"Material:\\n\")\n",
    "\n",
    "# zone 1\n",
    "inpf.write(\"  Zone_1:\\n\")\n",
    "inpf.write(\"    Type: PDState\\n\")\n",
    "inpf.write(\"    Horizon: %4.6e\\n\" % (horizon))\n",
    "inpf.write(\"    Density: %4.6e\\n\" % (rho1))\n",
    "inpf.write(\"    Compute_From_Classical: true\\n\")\n",
    "inpf.write(\"    K: %4.6e\\n\" % (K1))\n",
    "inpf.write(\"    G: %4.6e\\n\" % (G1))\n",
    "inpf.write(\"    Gc: %4.6e\\n\" % (Gc1))\n",
    "inpf.write(\"    Influence_Function:\\n\")\n",
    "inpf.write(\"      Type: 1\\n\")\n",
    "\n",
    "# zone 2\n",
    "inpf.write(\"  Zone_2:\\n\")\n",
    "inpf.write(\"    Type: PDState\\n\")\n",
    "inpf.write(\"    Horizon: %4.6e\\n\" % (horizon))\n",
    "inpf.write(\"    Density: %4.6e\\n\" % (rho2))\n",
    "inpf.write(\"    Compute_From_Classical: true\\n\")\n",
    "inpf.write(\"    K: %4.6e\\n\" % (K2))\n",
    "inpf.write(\"    G: %4.6e\\n\" % (G2))\n",
    "inpf.write(\"    Gc: %4.6e\\n\" % (Gc2))\n",
    "inpf.write(\"    Influence_Function:\\n\")\n",
    "inpf.write(\"      Type: 1\\n\")\n",
    "\n",
    "## Force\n",
    "# `Gravity` force is of body force type and all particles and walls experience this force\n",
    "inpf.write(\"Force_BC:\\n\")\n",
    "inpf.write(\"  Gravity: \" + print_dbl_list(g))\n",
    "\n",
    "## IC\n",
    "## Initial condition and boundary conditions can be applied in many ways\n",
    "## method 1: apply conditions on all particles in particular zone in which case one has to provide the \n",
    "##           Zone_List\n",
    "## method 2: apply conditions on specific particles by providing their global id. This is the current case.\n",
    "##           In Particle_List we prove `1` and say that apply constant velocity initial condition v0\n",
    "##           to particle 1. Particles are numbered starting from 0, so 1 correspons to p2 and 0 corresponds to \n",
    "##           p1\n",
    "## method 3: apply conditions on any particle withing specified region. So one can specify a rectangular region\n",
    "##           and apply conditions on any particle that falls in that region.\n",
    "inpf.write(\"IC:\\n\")\n",
    "inpf.write(\"  Constant_Velocity:\\n\")\n",
    "inpf.write(\"    Velocity_Vector: \" + print_dbl_list(v0))\n",
    "inpf.write(\"    Particle_List: [1]\\n\")\n",
    "\n",
    "## Displacement\n",
    "inpf.write(\"Displacement_BC:\\n\")\n",
    "# it is typical to have many different type of boundary conditions so it is conventient to \n",
    "# consider a sets of boundary condition (force boundary condition is similar)\n",
    "# in this case we just have 1 set\n",
    "inpf.write(\"  Sets: 1\\n\")\n",
    "inpf.write(\"  Set_1:\\n\")\n",
    "# this set is applicable to `0` particle (bottom)\n",
    "inpf.write(\"    Particle_List: [0]\\n\")\n",
    "# displacement components along x and y are impacted\n",
    "inpf.write(\"    Direction: [1,2]\\n\")\n",
    "# displacement is fixed by a function which is constant in time and space. Further, the value of this \n",
    "# constant is 0 (see Parameters block which expects single parameter for constant type function)\n",
    "inpf.write(\"    Time_Function:\\n\")\n",
    "inpf.write(\"      Type: constant\\n\")\n",
    "inpf.write(\"      Parameters:\\n\")\n",
    "inpf.write(\"        - 0.0\\n\")\n",
    "inpf.write(\"    Spatial_Function:\\n\")\n",
    "inpf.write(\"      Type: constant\\n\")\n",
    "# for `0` particle we have fixed all degree of freedoms and the displacement is fixed to 0\n",
    "# so we explicitly say that for this particle all defs are zero\n",
    "# this extra information allows us to not compute this boundary condition every time step as there is nothing\n",
    "# really to compute (displacemnts are always fixed to zero)\n",
    "inpf.write(\"    Zero_Displacement: true\\n\")\n",
    "\n",
    "## Output info\n",
    "## in this block, we take care of simulation output\n",
    "inpf.write(\"Output:\\n\")\n",
    "# where to write files\n",
    "inpf.write(\"  Path: %s\\n\" %(fpath))\n",
    "inpf.write(\"  Tags:\\n\")\n",
    "# list the variable that you want to output in the .vtu file\n",
    "inpf.write(\"    - Displacement\\n\")\n",
    "inpf.write(\"    - Velocity\\n\")\n",
    "inpf.write(\"    - Force\\n\")\n",
    "inpf.write(\"    - Force_Density\\n\")\n",
    "inpf.write(\"    - Damage_Z\\n\")\n",
    "inpf.write(\"    - Damage\\n\")\n",
    "inpf.write(\"    - Nodal_Volume\\n\")\n",
    "inpf.write(\"    - Zone_ID\\n\")\n",
    "inpf.write(\"    - Particle_ID\\n\")\n",
    "inpf.write(\"    - Fixity\\n\")\n",
    "inpf.write(\"    - Force_Fixity\\n\")\n",
    "inpf.write(\"    - Contact_Nodes\\n\")\n",
    "inpf.write(\"    - No_Fail_Node\\n\")\n",
    "inpf.write(\"    - Boundary_Node_Flag\\n\")\n",
    "inpf.write(\"    - Theta\\n\")\n",
    "# how frequent the output should take place\n",
    "inpf.write(\"  Output_Interval: %d\\n\" % (N_out))\n",
    "inpf.write(\"  Compress_Type: zlib\\n\")\n",
    "inpf.write(\"  Perform_FE_Out: false\\n\")\n",
    "inpf.write(\"  Perform_Out: true\\n\")\n",
    "# how frequent we dump the test specific outputs\n",
    "inpf.write(\"  Test_Output_Interval: %d\\n\" % (N_out))\n",
    "# what is verbosity level (last time I checked, the maximum is 3 so if you have 3 or above it will\n",
    "# produce maximum debug information)\n",
    "inpf.write(\"  Debug: 3\\n\")\n",
    "# assign this simulation a tag\n",
    "inpf.write(\"  Tag_PP: 0\\n\")\n",
    "\n",
    "## close file\n",
    "## damn I feel tired already!!\n",
    "inpf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76a56432",
   "metadata": {},
   "source": [
    "## 3 Running the simulation\n",
    "By now, we should have following files\n",
    "- `particle_locations.csv`\n",
    "- `p1.msh` and `p2.msh`\n",
    "- `input.yaml` .\n",
    "\n",
    "We run the simulation as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "504d9260",
   "metadata": {},
   "outputs": [],
   "source": [
    "PeriDEM = <path to peridem executable>/PeriDEM\n",
    "Nthreads = 4\n",
    "# comment out below once path to peridem is set correctly\n",
    "#!$PeriDEM -i input.yaml -nThreads $Nthreads"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}