# DEM Examples

## Collision of 2 particles

This tutorial aims to build the simulation for 2 bodies colliding in an elastic way. The user declares a sphero-cube and a sphero-tetrahedron in a collision course. The energy, liner and angular momenta can be measured to check these conservation laws. The whole simulation is described in the following code which will be explained in detail:// MechSys #include <mechsys/dem/domain.h> int main(int argc, char **argv) try { DEM::Domain dom; // add cube Vec3_t x(-10,0,0); Vec3_t v(1.,0,0); Vec3_t w(0,M_PI/5,0); dom.AddCube (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.3,/*Length*/3.,/*Density*/1.); dom.Particles[0]->v = v; dom.Particles[0]->w = w; // add tetrahedron x = 10.0, 0 , 0; w = 0.0, 0 , M_PI/10.; v = -1.0, 0 , 0; dom.AddTetra (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.5,/*Length*/5.,/*Density*/1.); dom.Particles[1]->v = v; dom.Particles[1]->w = w; // particle parameters Dict D; D.Set(-1,"Gn Gt Mu",0.0,0.0,0.0); dom.SetProps(D); // initial constants Vec3_t l0; Vec3_t p0; double Ek0,Ep0,E0; dom.LinearMomentum (p0); dom.AngularMomentum (l0); E0 = d.CalcEnergy (Ek0,Ep0); // solve dom.CamPos = 0.0,30.0,0.0; dom.Solve(/*tf*/30.0, /*dt*/1.0e-4, /*dtOut*/0.3, NULL, NULL, /*filekey*/"test_dynamics"); // final constants Vec3_t l1; // initial linear momentum Vec3_t p1; // initial angular momentum double Ek1,Ep1,E1; // initial energy dom.LinearMomentum (p1); dom.AngularMomentum (l1); E1 = dom.CalcEnergy (Ek1,Ep1); // check double tol = 1.0e-3; double err_l = norm(l1-l0); double err_p = norm(p1-p0); double err_E = fabs(E1-E0); double error = err_l + err_p + err_E; cout << "Error in energy = " << err_E << endl; cout << "Error in angular momentum = " << err_l << endl; cout << "Error in linear momentum = " << err_p << endl; // results if (error > tol) return 1; else return 0; } MECHSYS_CATCH

### Initialization

MechSys/DEM simulations are prepared with a C++ source code. The first step is to include the required MechSys libraries. In this case the file domain.h contains the DEM domain definitions:

#include <mechsys/dem/domain.h>

### Simulation script structure

The simulation script structure is:

#include ... int main(int argc, char **argv) try { ... } MECHSYS_CATCH

After including the required MechSys' libraries (and C++ standard libraries, if necessary), the `main`

function is declared.

For MechSys, it is important to add the `try`

keyword just after
the `int main()`

and before the first opening brace `{`

. Then, after the last closing
brace the `MECHSYS_CATCH`

is added.

### Declaring the DEM Domain and Particles

The line,

DEM::Domain d;

declares the container `dom`

of the class
`DEM::Domain`

. This class defines the DEM universe and also has many functions
to include particles, some of which will be used in this tutorial.

To declare the cube, some 3D vectors are needed. The `Vec3_t`

data type is used
for fast vectorial operations:

Vec3_t x(-10,0,0); Vec3_t v(1.,0,0); Vec3_t w(0,M_PI/5,0);

The vectors `x`

`v`

`w`

define respectively the cube center position, its linear velocity and its angular velocity. The declared angular velocity ensures that the cube will
spin around the y-axis counter clockwise. Once these three vectors are declared, the cube is introduced into the domain:

dom.AddCube (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.3,/*Length*/3.,/*Density*/1.);

This function creates a cube into the memory. The first argument is a Tag which will be useful later on. The second is the position vector declared
above. Next is the Sphero-Radius which is the radius of the rounded corners, and is required for the collision law. Ideally this spheroradius should
be large enough so the overlapping between the colliding particles is small compared to the particles dimensions. The 4^{th} argument is the length of
the cube's edges and finally the density of the cube is also declared.

The method, however, initializes the particle at rest and without rotation. Therefore, to add the vector velocities, internal variables must be accessed:

dom.Particles[0]->v = v; dom.Particles[0]->w = w;

The object `dom.Particles`

is an array that contains the memory pointers to all the particles created in the domain `dom`

. To access
the first particle, which is the cube, the index `0`

must be introduced. Therefore `dom.Particles[0]`

is a pointer to the cube. The internal vector variables `v`

and `w`

are the cube's linear and angular velocity vectors.

In a similar way the tetrahedron is included at the right side of the cube moving towards it:

x = 10.0, 0 , 0; w = 0.0, 0 , M_PI/10.; v = -1.0, 0 , 0; dom.AddTetra (/*Tag*/-1, /*Position*/x,/*SpheroRadius*/0.5,/*Length*/5.,/*Density*/1.); dom.Particles[1]->v = v; dom.Particles[1]->w = w;

With both particles properly defined, some collision parameters must be set. MechSys uses the spheropolyhedra method to model the collision of
solid objects. The collision law is defined by two spring constants `K`

and _{n}`K`

; two dissipative constants _{t}`G`

and _{n}`G`

; and a friction coefficient μ. Since this simulation deals with an elastic collision, the
two dissipative constants and the friction coefficient must be set to zero. In MechSys this is achieved by the use of dictionaries as follows:_{t}

Dict D; D.Set(-1,"Gn Gt Mu",0.0,0.0,0.0); dom.SetProps(D);

Where the dictionary `D`

contains the user information on the new parameters. The first argument of the
function `D.Set`

is the particle's tag. Hence different groups of particles can be associated with different
collision parameters. The elastic parameters `K`

and _{n}`K`

(as _{t}`Kn`

and `Kt`

in the dictionary string) can also be included into the dictionary, however for this example we will work
with their default values.

### Starting the simulation

To check the conservation laws some values must be calculated before the simulation starts:

Vec3_t l0; Vec3_t p0; double Ek0,Ep0,E0; dom.LinearMomentum (p0); dom.AngularMomentum (l0); E0 = d.CalcEnergy (Ek0,Ep0);

The vectors `l0`

and `p0`

will store the linear and angular momenta given by the functions `dom.LinearMomentum`

and `dom.AngularMomentum`

. The three numbers `Ek0`

, `Ep0`

and `E0`

will be respectively the total kinetic energy, the
total potential energy and the total energy as given by the function `dom.CalcEnergy`

. This variable will be kept in memory to be compared with the situation after the simulation.

With the DEM universe defined, the simulation will begin. First, for future visualization, is useful to define the position of a camera:

dom.CamPos = 0.0,30.0,0.0;

Then the system evolves for as long as defined by the user by means of the `dom.Solve`

:

dom.Solve(/*tf*/30.0, /*dt*/1.0e-4, /*dtOut*/0.3, NULL, NULL, /*filekey*/"test_dynamics");

The system will evolve in this case up to `tf = 30.0`

seconds of simulation time. The next
argument is the integration time step `dt`

, the smaller the more accurate the simulation is but also the more
time consuming. The 3^{rd} argument is the report time step, in other words, the time span between movie frames. The next two arguments are
pointers to functions which are not explained in this simple example. This functions may be called from inside the Solve function, and hence offer
some control over the evolution of the system from outside. Finally, the last argument is a file key which is used to identify the report and
animation files form this simulation.

At the end of the simulation, the same variables that were measured before it are measured again:

Vec3_t l1; Vec3_t p1; double Ek1,Ep1,E1; dom.LinearMomentum (p1); dom.AngularMomentum (l1); E1 = d.CalcEnergy (Ek1,Ep1);

And compared with their initial values:

double tol = 1.0e-3; double err_l = norm(l1-l0); double err_p = norm(p1-p0); double err_E = fabs(E1-E0); double error = err_l + err_p + err_E; cout << "Error in energy = " << err_E << endl; cout << "Error in angular momentum = " << err_l << endl; cout << "Error in linear momentum = " << err_p << endl;

After compiling and running, the simulation should give this results

Error in energy = 0.128839 Error in angular momentum = 0.000122599 Error in linear momentum = 8.30147e-10

Unfortunately the conservation laws are not exactly met in the DEM and therefore there is a difference between the initial and final values. As mentioned before this error can be tuned down with the integration time step

### Simulation visualization

MechSys uses VisIt to visualize simulation results. A full tutorial can be found here. After following the instructions the dynamics of the system should resemble the following video: