# DEM Examples

```                                    #include <mechsys/dem/domain.h>
#include <mechsys/util/fatal.h>
#include <mechsys/util/util.h>

using std::cout;
using std::endl;
using Util::PI;
using Util::FmtErr;
using DEM::Domain;

int main(int argc, char **argv) try
{
// test the mesh particle generation
double R = 0.1;
Domain d;
d.AddFromJson(-1, "test.msh", R, 1.0, 1.0);
d.Particles[0]->Position(OrthoSys::O);
double epsilon = R/pow(d.Particles[0]->Props.V,1.0/3.0);

std::cout << "[1;33m\n--- Results for mass properties from surface integrals --------------------------------------[0m\n";
cout << "  Volume            = " << d.Particles[0]->Props.V << "\n";
cout << "  Center of mass    = " << d.Particles[0]->x(0)    << ", " << d.Particles[0]->x(1) << ", " << d.Particles[0]->x(2) << "\n";
cout << "  Moment of inertia = " << d.Particles[0]->I(0)    << ", " << d.Particles[0]->I(1) << ", " << d.Particles[0]->I(2) << "\n";
cout << "  Quaternion        = " << d.Particles[0]->Q(0)    << ", " << d.Particles[0]->Q(1) << ", " << d.Particles[0]->Q(2) << ", " << d.Particles[0]->Q(3) << "\n";
cout << endl;

d.Particles[0]->Props.V *= 1.0 + 3.0*epsilon;
d.Particles[0]->I       *= 1.0 + 5.0*epsilon;

std::cout << "[1;33m\n--- Results for mass properties from surface integrals corrected by spheroradius ------------[0m\n";
cout << "  Volume            = " << d.Particles[0]->Props.V << "\n";
cout << "  Center of mass    = " << d.Particles[0]->x(0)    << ", " << d.Particles[0]->x(1) << ", " << d.Particles[0]->x(2) << "\n";
cout << "  Moment of inertia = " << d.Particles[0]->I(0)    << ", " << d.Particles[0]->I(1) << ", " << d.Particles[0]->I(2) << "\n";
cout << "  Quaternion        = " << d.Particles[0]->Q(0)    << ", " << d.Particles[0]->Q(1) << ", " << d.Particles[0]->Q(2) << ", " << d.Particles[0]->Q(3) << "\n";
cout << endl;

d.Particles[0]->CalcProps(10000);
// output
std::cout << "[1;33m\n--- Results for mass properties from Monte Carlo integration -------------------------------[0m\n";
cout << "  Volume            = " << d.Particles[0]->Props.V << "\n";
cout << "  Center of mass    = " << d.Particles[0]->x(0)    << ", " << d.Particles[0]->x(1) << ", " << d.Particles[0]->x(2) << "\n";
cout << "  Moment of inertia = " << d.Particles[0]->I(0)    << ", " << d.Particles[0]->I(1) << ", " << d.Particles[0]->I(2) << "\n";
cout << "  Quaternion        = " << d.Particles[0]->Q(0)    << ", " << d.Particles[0]->Q(1) << ", " << d.Particles[0]->Q(2) << ", " << d.Particles[0]->Q(3) << "\n";
cout << endl;

// draw
d.WriteXDMF ("test_mesh");
return 0;
}
MECHSYS_CATCH
```

The example is simpler than the previous one on the collision of two particles. The key component of this example is the function `Domain::AddFromJson`,

```                                        d.AddFromJson(-1, "test.msh", R, 1.0, 1.0);
```

The first argument is the tag. As explained in the previous example, the tags are used mainly to assign different properties to different sets of particles. The second argument is the name of the file that contains the mesh data to generate the particle. Please open 'test.msh' to see its structure. If you wish to produce your own mesh file using the open source Blender software, please follow Pei Zhang's excellent tutorial that you can find here. Additionally you will need the python plugin to properly export the blender mesh into the Mechsys format .The third argument is the spheroradius which is the key parameter for the collision law (please see the references on the method). The fourth argument is the density which is irrelevant in this example since the particle is not going to move. The fifth argument is the scale factor which is the multiplier of the polyhedron dimensions. There are some more arguments taking default values, please see `mechsys/lib/domain.h` for details.

The next two lines of code determine the position and correction factors for the particles mass properties,

```                                        d.Particles[0]->Position(OrthoSys::O);
double epsilon = R/pow(d.Particles[0]->Props.V,1.0/3.0);
```

The first one puts the center of mass of the particle in the position `Orthosys::O` which is the zero vector in Mechsys (Alternatively `Vec3_t(0,0,0)` could be used). The second line declares a correction parameter epsilon which accounts for the spheroradius dilation. The variable `Props.V` Is the volume of the Particle.

Next some of the mass properties are reported. Mechsys calculates the mass properties of DEM particles in two ways. The first one and fastest is by using Green's theorem to turn the volume integrals into surface integrals over each face of the polyhedron. Mass properties are the volume, the center of mass `x`, the principal values of the moment of inertia tensor `I` (represented as a Vec3_t) and the quaternion representing the current orientation `Q`

```                                        std::cout << "[1;33m\n--- Results for mass properties from surface integrals --------------------------------------[0m\n";
cout << "  Volume            = " << d.Particles[0]->Props.V << "\n";
cout << "  Center of mass    = " << d.Particles[0]->x(0)    << ", " << d.Particles[0]->x(1) << ", " << d.Particles[0]->x(2) << "\n";
cout << "  Moment of inertia = " << d.Particles[0]->I(0)    << ", " << d.Particles[0]->I(1) << ", " << d.Particles[0]->I(2) << "\n";
cout << "  Quaternion        = " << d.Particles[0]->Q(0)    << ", " << d.Particles[0]->Q(1) << ", " << d.Particles[0]->Q(2) << ", " << d.Particles[0]->Q(3) << "\n";
cout << endl;
```

Such quantities can be corrected by the correction factor `epsilon` to account for the dilation of the spheroradious. These corrections should be accurate enough provided the spheroradious is small compared with the dimensions of the polyhedral skeleton.

```                                        d.Particles[0]->Props.V *= 1.0 + 3.0*epsilon;
d.Particles[0]->I       *= 1.0 + 5.0*epsilon;
```

Finally, after reporting them, the last method of mass properties calculation is used, the Monte Carlo integration which is slow but already integrates over the volume of the polyhedral skeleton dilated by the spheroradious. Such calculation is carried out by the `Particle::CalcProps` method of the Particle class with an argument giving the number of iterations.

```                                        d.Particles[0]->CalcProps(10000);
```

After reporting the mass properties from all these methods, the polyhedra is draw so it can be visualized by the function `Domain::WriteXDMF`:

```                                        d.WriteXDMF ("test_mesh");
```

This function is called internally when the dynamic loop (function `Solve` of the previous example) is called as well.