# 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;

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

x =  10.0, 0 , 0;
w =   0.0, 0 , M_PI/10.;
v =  -1.0, 0 , 0;
dom.Particles->v = v;
dom.Particles->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 4th 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->v = v;
dom.Particles->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` 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.Particles->v = v;
dom.Particles->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 `Kn` and `Kt`; two dissipative constants `Gn` and `Gt`; 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:

```                                    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 `Kn` and `Kt` (as `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 3rd 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: