DEM Examples
Loading a particle from mesh
#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.