Beginner tutorials

This page contains the following tutorials:

# T1: Particle motion in outer space

Particle moving with a constant velocity in outer space.

#### Problem description:

The first tutorial is setup to simulate a particle moving with a constant velocity in the absence of gravity. You can find the simulation setup in Tutorial1_ParticleInOuterSpace.cpp The detailed description of this tutorial is presented below.

### Headers:

The headers are necessary to set up the simulations. They come from the kernel and standard libraries (see Developing a kernel feature). Headers files usually have a .h or .hpp extensions. For example:

### Before main function:

A class is the guide for objects to connect with handlers.
You can find in the code:

class Tutorial1 : public Mercury3D
This adds on the hierarchical grid code for 3D problems.
Definition: Mercury3D.h:37
[T1:headers]
Definition: Tutorial1_ParticleInOuterSpace.cpp:44

Here is defined the class Tutorial1, which is inherited from the main class Mercury3D. Then, you can call the function to fit the initial conditions of the problem:

void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: Tutorial1_ParticleInOuterSpace.cpp:48

It creates the problems’ initial conditions to define the radius, position, velocity of particles. In order to create a particle, the following structure must be used:

p0.setSpecies(speciesHandler.getObject(0));
p0.setRadius(0.05); // sets particle radius
p0.setPosition(Vec3D(0.1 * getXMax(), 0.1 * getYMax(), 0.1 * getZMax())); // sets particle position
p0.setVelocity(Vec3D(0.5, 0.1, 0.1));// sets particle velocity
particleHandler.copyAndAddObject(p0);
void setVelocity(const Vec3D &velocity)
set the velocity of the BaseInteractable.
Definition: BaseInteractable.cc:350
virtual void setPosition(const Vec3D &position)
Sets the position of this BaseInteractable.
Definition: BaseInteractable.h:239
virtual void setRadius(Mdouble radius)
Sets the particle's radius_ (and adjusts the mass_ accordingly, based on the particle's species)
Definition: BaseParticle.cc:553
void setSpecies(const ParticleSpecies *species)
Definition: BaseParticle.cc:818
A spherical particle is the most simple particle used in MercuryDPM.
Definition: SphericalParticle.h:37
Definition: Vector.h:51

The particle data members as positions, radius and velocities are pointed by using an object (P0) to their class member. For the case of the mentioned data members, the class member is called BaseParticle. Then, the defined information is sent to its corresponding handler: particleHandler by using the function copyAndAddObject

### Main function:

In the main function, the global parameters of the problem are defined. It includes gravity, spatial dimensions (x,y,z), total run time, the type of contact law, etc.
The object 'problem' is an instance of the defined class Tutorial1. Global parameters and species will point to the object 'problem'.

Tutorial1 problem;

Then, the data members are defined as:

problem.setName("Tutorial1");
problem.setSystemDimensions(3);
problem.setGravity(Vec3D(0.0, 0.0, 0.0));
problem.setXMax(1.0);
problem.setYMax(1.0);
problem.setZMax(1.0);
problem.setTimeMax(2.0);

Subsequently, the Species of the problem are defined. It means, the properties (density,stiffness) and the corresponding contact law (for this tutorial, see LinearViscoelasticSpecies). Initially, when a particle is created, it attains the properties of a default species type with ‘0’ as its index.

//Set the species of particles and walls
//The normal spring stiffness and normal dissipation is computed and set as
//For collision time tc=0.005 and restitution coefficeint rc=1.0
species.setDensity(2500.0); //sets the species type_0 density
species.setStiffness(258.5);//sets the spring stiffness.
species.setDissipation(0.0); //sets the dissipation.
problem.speciesHandler.copyAndAddObject(species);
void setDissipation(Mdouble dissipation)
Allows the normal dissipation to be changed.
Definition: LinearViscoelasticNormalSpecies.cc:117
void setStiffness(Mdouble new_k)
Allows the spring constant to be changed.
Definition: LinearViscoelasticNormalSpecies.cc:93
void setDensity(Mdouble density)
Definition: ParticleSpecies.cc:108

Here, 'species' is an object of the class: LinearViscoelasticSpecies. Data members are fitted to point to this object, and then, coppied and added to their corresponding handler speciesHandler.

Outputs:
Data output is vital to analyse simulations, which leads to defining ways to save the simulation data for post processing.
The simulations generate several types of data files.

problem.setSaveCount(10);
problem.dataFile.setFileType(FileType::ONE_FILE);
problem.restartFile.setFileType(FileType::ONE_FILE);
problem.fStatFile.setFileType(FileType::NO_FILE);
problem.eneFile.setFileType(FileType::NO_FILE);
logger(INFO, "run number: %", problem.dataFile.getCounter());
@ NO_FILE
file will not be created/read
@ ONE_FILE
all data will be written into/ read from a single file called name_
LL< Log::INFO > INFO
Info log level.
Definition: Logger.cc:55
Logger< MERCURYDPM_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here.

Thereby, data members and member functions are pointed to the object 'problem'.
For XBalls users, additional display options can be set as below

problem.setXBallsAdditionalArguments("-solidf -v0");

After all the simulation parameters are set, we reach a point where we put all the above bits of code in action by the following statements

problem.setTimeStep(0.005 / 50.0); // (collision time)/50.0
problem.solve(argc, argv);

# T2: Particle motion on earth

Particle falling due to gravity.

#### Problem description:

In Tutorial2_ParticleUnderGravity.cpp, we simulate a particle when dropped under the influence of gravity, $$g = 9.81 m/s^2$$. Basically, this tutorial is an extension of Tutorial1_ParticleInOuterSpace.cpp with few minor changes.
All we need to do is to change the initial particle position and velocity in the inherited class Tutorial2.

p0.setSpecies(speciesHandler.getObject(0));
p0.setRadius(0.05);
p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), getZMax()));
p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
particleHandler.copyAndAddObject(p0);

Also, the gravity vector in the main function:

problem.setName("Tutorial2");
problem.setSystemDimensions(3);
problem.setGravity(Vec3D(0.0, 0.0, -9.81));
problem.setXMax(1.0);
problem.setYMax(1.0);
problem.setZMax(5.0);
problem.setTimeMax(1.5);

# T3: Bouncing ball (elastic)

Particle bouncing off the blue wall.

#### Problem description:

The Tutorial3_BouncingBallElastic.cpp simulates a particle bouncing off a wall. It is assumed the collistion between the particle and the wall is elastic by implying that the restitution coefficient is unity. It means that the particle velocity before and after collision remains the same. Additionally, we will learn how to add a wall over which the ball bounces.

### Headers:

The header InfiniteWall.h is included.

### Before the main class:

The class InfiniteWall is included in the inherited class Tutorial3:

class Tutorial3 : public Mercury3D
{
public:
void setupInitialConditions() override {
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), getZMax()));
p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0, 0, getZMin()));
}
};
std::enable_if<!std::is_pointer< U >::value, U * >::type copyAndAddObject(const U &object)
Creates a copy of a Object and adds it to the BaseHandler.
Definition: BaseHandler.h:379
T * getObject(const unsigned int id)
Gets a pointer to the Object at the specified index in the BaseHandler.
Definition: BaseHandler.h:613
void setSpecies(const ParticleSpecies *species)
Defines the species of the current wall.
Definition: BaseWall.cc:169
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax.
Definition: DPMBase.h:626
SpeciesHandler speciesHandler
A handler to that stores the species type i.e. LinearViscoelasticSpecies, etc.
Definition: DPMBase.h:1427
WallHandler wallHandler
An object of the class WallHandler. Contains pointers to all the walls created.
Definition: DPMBase.h:1447
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created.
Definition: DPMBase.h:1437
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax.
Definition: DPMBase.h:638
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax.
Definition: DPMBase.h:650
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin.
Definition: DPMBase.h:644
A infinite wall fills the half-space {point: (position_-point)*normal_<=0}.
Definition: InfiniteWall.h:48
void set(Vec3D normal, Vec3D point)
Defines a standard wall, given an outward normal vector s.t. normal*x=normal*point for all x of the w...
Definition: InfiniteWall.cc:118
[T3:headers]
Definition: Tutorial3_BouncingBallElastic.cpp:44
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: Tutorial3_BouncingBallElastic.cpp:47

Where the object "w0" is an instance of the class InfiniteWall. Then, the function "copyAndAddObject" coppy and add the object to its corresponding handler. The above set of statements, create and place the wall at $$Z_{min}$$.
Note: Don’t forget to include the InfiniteWall.h header, as shown in the header section. In some sense, creation and addition of a wall is similar to creation and addition of a particle.

# T4: Bouncing ball with dissipation (inelastic)

Particle bouncing off the blue wall with restitution coefficient.

#### Problem description:

In Tutorial4_BouncingBallInelastic.cpp, the difference between an elastic and inelastic collision between a particle and a wall is illustrated. The only difference between Tutorial3_BouncingBallElastic.cpp and Tutorial4_BouncingBallInelastic.cpp is the value of the restitution coefficient:

double rc = 0.88; // restitution coefficient

See Bouncing ball - inelastic (code) for more details.

# T5: Elastic collision (2 particles)

Particles moving towards each other.

Problem description:
So far, in the above tutorials, we have seen how a particle and a wall interact during a collision. In this tutorial, we illustrate how two particles interact using Tutorial5_ElasticCollision.cpp. For this purpose, we need two particles. The particles may or may not be of the same species type. But, here we shall assume they are of same species and same size.

### Before the main function

class Tutorial5 : public Mercury3D
{
public:
void setupInitialConditions() override {
p0.setRadius(0.005); //particle-1 radius
p0.setPosition(Vec3D(0.25 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
p0.setVelocity(Vec3D(0.25, 0.0, 0.0));
particleHandler.copyAndAddObject(p0); // 1st particle created
p0.setRadius(0.005); // particle-2 radius
p0.setPosition(Vec3D(0.75 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
p0.setVelocity(Vec3D(-0.25, 0.0, 0.0));
particleHandler.copyAndAddObject(p0); // 2nd particle created
}
};
[T5:headers]
Definition: Tutorial5_ElasticCollision.cpp:43
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: Tutorial5_ElasticCollision.cpp:46

On comparison between the above class and class Tutorial1, we see how an additional particle is added. Two particles are created, and positioned oppositely apart at a certain distance between them. Both the particles, have a certain velocity directing them towards each other.

# T6: Elastic collisions with periodic boundaries

(a) Illustrates the idea behind periodic boundaries, particle exiting boundary b2 re-enters through boundary b1 (b) Illustrates the problem setup.

#### Problem description:

In order to have multiple collisions, Tutorial5_ElasticCollision.cpp is fitted with periodic boundaries in X.

### Headers:

The header PeriodicBoundary.h is included.

### Before the main function:

class Tutorial6 : public Mercury3D
{
public:
void setupInitialConditions() override {
p0.setRadius(0.005);//particle-1
p0.setPosition(Vec3D(0.25 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
p0.setVelocity(Vec3D(0.25, 0.0, 0.0));
p0.setRadius(0.005);// particle-2
p0.setPosition(Vec3D(0.75 * getXMax(), 0.5 * getYMax(), 0.5 * getZMax()));
p0.setVelocity(Vec3D(-0.25, 0.0, 0.0));
b0.set(Vec3D(1.0, 0.0, 0.0), getXMin(), getXMax());
}
};
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin.
Definition: DPMBase.h:619
BoundaryHandler boundaryHandler
An object of the class BoundaryHandler which concerns insertion and deletion of particles into or fro...
Definition: DPMBase.h:1452
Defines a pair of periodic walls. Inherits from BaseBoundary.
Definition: PeriodicBoundary.h:41
void set(Vec3D normal, Mdouble distanceLeft, Mdouble distanceRight)
Defines a PeriodicBoundary by its normal and positions.
Definition: PeriodicBoundary.cc:84
[T6:headers]
Definition: Tutorial6_ElasticCollisionPeriodic.cpp:44
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: Tutorial6_ElasticCollisionPeriodic.cpp:47

In the Class Tutorial6
(i) we create two particles of same type and different sizes. (ii) we setup periodic boundaries in X-direction as

b0.set(Vec3D(1.0, 0.0, 0.0), getXMin(), getXMax());
boundaryHandler.copyAndAddObject(b0);

Where "b0" is an instance of the class PeriodicBoundary. Then, the object is copied and added to its corresponding handler: boundaryHandler.

# T7: Motion of a particle in a two dimensional (2D) box

Particle motion in a box (blue and black denote the walls).

#### Problem description:

In previous tutorials, we have seen how a particle interacts with a wall and itself. In this tutorial, we will design boxes of different shapes by using more than one wall. As an example, in absence of gravity, we will simulate a particle moving in a two dimensional square shaped box. We consider two dimensions only for the sake of simplicity. The same idea was introduced in Tutorial3_BouncingBallElastic.cpp.

### Before the main function:

class Tutorial7 : public Mercury3D
{
public:
void setupInitialConditions() override {
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.5 * getXMax(), 0.5 * getYMax(), 0.0));
p0.setVelocity(Vec3D(1.2, 1.3, 0.0));
w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0.0, 0.0));
w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0.0, 0.0));
w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, getYMax(), 0.0));
w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, getYMin(), 0.0));
}
};
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin.
Definition: DPMBase.h:632
[T7:headers]
Definition: Tutorial7_ParticleIn2DBox.cpp:44
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: Tutorial7_ParticleIn2DBox.cpp:47

In this class, we setup a 2D square shaped box or a polygon by adding more walls as shown above. In total, we have 4 walls forming our box within which the particle will traverse. Note: As we simulate in 2D, no walls are set in z-direction.

### Main function:

As our simulation is two dimensional, we set the system dimensions as 2

void setSystemDimensions(unsigned int newDim)
Sets the system dimensionality.
Definition: DPMBase.cc:1417

The complete code for the above problem description can be found in Tutorial7_ParticleIn2DBox.cpp.

# T8: Motion of a particle in a box with an obstacle

#### Problem description:

We extend the problem setup of Tutorial 7, by adding a rectangular block as shown in the above figure. To create this block of wall or obstacle, we will use the Class IntersectionOfWalls. Before we go ahead it is advised to know the difference between an infinite wall and finite wall, see Different types of walls. As an example, we create an obstacle using a set of finite walls and place it within the box created using a set of infinite walls. See the above figure.

### Headers:

The header IntersectionOfWalls.h is included.

### Before the main function:

To add intersection walls, we use:

w1.setSpecies(speciesHandler.getObject(0));
w1.addObject(Vec3D(-1.0, 0.0, 0.0), Vec3D(0.75 * getXMax(), 0.0, 0.0));
w1.addObject(Vec3D(1.0, 0.0, 0.0), Vec3D(0.25 * getXMax(), 0.0, 0.0));
w1.addObject(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, 0.75 * getYMax(), 0.0));
w1.addObject(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, 0.25 * getYMax(), 0.0));
wallHandler.copyAndAddObject(w1);
A IntersectionOfWalls is convex polygon defined as an intersection of InfiniteWall's.
Definition: IntersectionOfWalls.h:59
void addObject(Vec3D normal, Vec3D point)
Adds a wall to the set of infinite walls, given a normal vector pointing into the wall (i....
Definition: IntersectionOfWalls.cc:138
void setSpecies(const ParticleSpecies *species)
sets species of subwalls as well
Definition: IntersectionOfWalls.cc:72

where "w1" is an instance of the class IntersectionOfWalls. This surface (wall) will have the properties of the index 0: "getObject(0)". Then, the object is coppied and added to its corresponding handler: wallHandler.

# T9: Motion of a ball over an inclined plane (Sliding + Rolling)

#### Problem description:

the motion of three particles over an inclined plane is presented in Tutorial9_InclinedPlane.cpp. Here, the particles are fitted with different features.

### Headers:

For Tutorial9_InclinedPlane.cpp, the header which implements the elastic properties and contact force modelc has been changed. In order to include the rolling and sliding effects, LinearViscoelasticFrictionSpecies.h is used.

### Before the main function:

class Tutorial9 : public Mercury3D
{
public:
void setupInitialConditions() override {
//sets the particle to species type-1
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.05 * getXMax(), 0.01 * getYMax(), getZMin() + p0.getRadius()));
p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
//sets the particle to species type-2
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.05 * getXMax(), 0.21 * getYMax(), getZMin() + p0.getRadius()));
p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
//sets the particle to species type-3
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.05 * getXMax(), 0.41 * getYMax(), getZMin() + p0.getRadius()));
p0.setVelocity(Vec3D(0.0, 0.0, 0.0));
w0.set(Vec3D(0.0, 0.0, -1.0), Vec3D(0.0, 0.0, getZMin()));
}
};
Mdouble getRadius() const
Returns the particle's radius.
Definition: BaseParticle.h:348
[T9:headers]
Definition: Tutorial9_InclinedPlane.cpp:44
void setupInitialConditions() override
This function allows to set the initial conditions for our problem to be solved, by default particle ...
Definition: Tutorial9_InclinedPlane.cpp:47

Three particles are fitted with radius, initial positions and initial velocities.Then, all particles are coppied and added to their corresponding handler.
It is defined the object 'w0' as an instance of the class InfiniteWall. For this tutorial, the surface will have the same properties of "object0". It means the surface and the first particle will have the same elastic properties.

### Main function:

int main(int argc, char* argv[])
{
// Problem setup
Tutorial9 problem;//instantiate an object of class Tutorial 9
double angle = constants::pi / 180.0 * 20.0;
problem.setName("Tutorial9");
problem.setSystemDimensions(3);
problem.setGravity(Vec3D(sin(angle), 0.0, -cos(angle)) * 9.81);
problem.setXMax(0.3);
problem.setYMax(0.3);
problem.setZMax(0.05);
problem.setTimeMax(0.5);
// The normal spring stiffness and normal dissipation is computed and set as
// For collision time tc=0.005 and restitution coefficeint rc=0.88,
//Properties index 0
species0.setDensity(2500.0);//sets the species type_0 density
species0.setStiffness(259.018);//sets the spring stiffness.
species0.setDissipation(0.0334);//sets the dissipation.
species0.setSlidingStiffness(2.0 / 7.0 * species0.getStiffness());
species0.setRollingStiffness(2.0 / 5.0 * species0.getStiffness());
species0.setSlidingFrictionCoefficient(0.0);
species0.setRollingFrictionCoefficient(0.0);
auto ptrToSp0=problem.speciesHandler.copyAndAddObject(species0);
//Properties index 1
species1.setDensity(2500.0);//sets the species type-1 density
species1.setStiffness(259.018);//sets the spring stiffness
species1.setDissipation(0.0334);//sets the dissipation
species1.setSlidingStiffness(2.0 / 7.0 * species1.getStiffness());
species1.setRollingStiffness(2.0 / 5.0 * species1.getStiffness());
species1.setSlidingFrictionCoefficient(0.5);
species1.setRollingFrictionCoefficient(0.0);
auto ptrToSp1=problem.speciesHandler.copyAndAddObject(species1);
//Combination of properties index 0 and index 1
auto species01 = problem.speciesHandler.getMixedObject(ptrToSp0,ptrToSp1);
species01->setStiffness(259.018);//sets the spring stiffness
species01->setDissipation(0.0334);//sets the dissipation
species01->setSlidingStiffness(2.0 / 7.0 * species01->getStiffness());
species01->setRollingStiffness(2.0 / 5.0 * species01->getStiffness());
species01->setSlidingFrictionCoefficient(0.5);
species01->setRollingFrictionCoefficient(0.0);
//Properties index 2
species2.setDensity(2500.0);//sets the species type-2 density
species2.setStiffness(258.5);//sets the spring stiffness
species2.setDissipation(0.0);//sets the dissipation
species2.setSlidingStiffness(2.0 / 7.0 * species2.getStiffness());
species2.setRollingStiffness(2.0 / 5.0 * species2.getStiffness());
species2.setSlidingFrictionCoefficient(0.5);
species2.setRollingFrictionCoefficient(0.5);
auto ptrToSp2 = problem.speciesHandler.copyAndAddObject(species2);
//Combination of properties index 0 and index 2
auto species02 = problem.speciesHandler.getMixedObject(ptrToSp0, ptrToSp2);
species02->setStiffness(259.018);//sets the stiffness
species02->setDissipation(0.0334);//sets the dissipation
species02->setSlidingStiffness(2.0 / 7.0 * species02->getStiffness());
species02->setRollingStiffness(2.0 / 5.0 * species02->getStiffness());
species02->setSlidingFrictionCoefficient(0.5);
species02->setRollingFrictionCoefficient(0.5);
//Output
problem.setSaveCount(10);
problem.setXBallsAdditionalArguments("-solidf -v0 -s 8");
problem.setTimeStep(0.005 / 50.0);
problem.solve(argc, argv);
return 0;
}
int main(int argc, char *argv[])
Definition: T_protectiveWall.cpp:215
void setSaveCount(unsigned int saveCount)
Sets File::saveCount_ for all files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:408
File eneFile
An instance of class File to handle in- and output into a .ene file.
Definition: DPMBase.h:1488
File fStatFile
An instance of class File to handle in- and output into a .fstat file.
Definition: DPMBase.h:1483
void setName(const std::string &name)
Allows to set the name of all the files (ene, data, fstat, restart, stat)
Definition: DPMBase.cc:422
File dataFile
An instance of class File to handle in- and output into a .data file.
Definition: DPMBase.h:1478
void setYMax(Mdouble newYMax)
Sets the value of YMax, the upper bound of the problem domain in the y-direction.
Definition: DPMBase.cc:1191
File restartFile
An instance of class File to handle in- and output into a .restart file.
Definition: DPMBase.h:1493
void setXBallsAdditionalArguments(std::string newXBArgs)
Set the additional arguments for xballs.
Definition: DPMBase.cc:1347
void setXMax(Mdouble newXMax)
Sets the value of XMax, the upper bound of the problem domain in the x-direction.
Definition: DPMBase.cc:1165
void setZMax(Mdouble newZMax)
Sets the value of ZMax, the upper bound of the problem domain in the z-direction.
Definition: DPMBase.cc:1217
void setTimeStep(Mdouble newDt)
Sets a new value for the simulation time step.
Definition: DPMBase.cc:1234
void setTimeMax(Mdouble newTMax)
Sets a new value for the maximum simulation duration.
Definition: DPMBase.cc:873
void solve()
The work horse of the code.
Definition: DPMBase.cc:4270
void setGravity(Vec3D newGravity)
Sets a new value for the gravitational acceleration.
Definition: DPMBase.cc:1383
void setFileType(FileType fileType)
Sets the type of file needed to write into or read from. File::fileType_.
Definition: File.cc:215
std::enable_if<!std::is_pointer< typename U::MixedSpeciesType >::value, typename U::MixedSpeciesType * >::type getMixedObject(const U *S, const U *T)
Definition: SpeciesHandler.h:74
const Mdouble pi
Definition: ExtendedMath.h:45
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
Mdouble sin(Mdouble x)
Definition: ExtendedMath.cc:44

The main function of Tutorial9_InclinedPlane.cpp presents how the elastic features of particles are defined. Also, a new function is added to combine features of different indexes: getMixedObject.

# T11: Axisymmetric walls inside a cylindrical domain

#### Problem description:

This tutorial explains how to include axisymmetric walls inside a cylindrical domain. Also, the user will find an useful implementation of how to create polydisperse particle size distribution and perform internal validation tests. The full code can be found in Tutorial11_AxisymmetricWalls.cpp

### Headers:

To work with axisymmetric walls, the user needs to add the following headers:

### Before the main function:

The class Tutorial 11 contains the constructor Tutorial 11

Tutorial11(const Mdouble width, const Mdouble height){
logger(INFO, "Tutorial 11");
setName("Tutorial11");
setFileType(FileType::ONE_FILE);
restartFile.setFileType(FileType::ONE_FILE);
fStatFile.setFileType(FileType::NO_FILE);
eneFile.setFileType(FileType::NO_FILE);
setParticlesWriteVTK(true);
wallHandler.setWriteVTK(FileType::ONE_FILE);
setXBallsAdditionalArguments("-v0 -solidf");
//specify body forces
setGravity(Vec3D(0.0, 0.0, -9.8));
//set domain accordingly (domain boundaries are not walls!)
setXMin(0.0);
setXMax(width);
setYMin(0.0);
setYMax(width);
setZMin(0.0);
setZMax(height);
}
double Mdouble
Definition: GeneralDefine.h:34
[T11:headers]
Definition: Tutorial11_AxisymmetricWalls.cpp:38

This constructor is an instance of the main class. Here it's defined the global output and non-changeable parameters.

In this tutorial, it's overrided the member function setupInitialConditions(), which doesn't return any value (defined then as void).

void setupInitialConditions() override {
Vec3D mid = {
(getXMin() + getXMax()) / 2.0,
(getYMin() + getYMax()) / 2.0,
(getZMin() + getZMax()) / 2.0};
const Mdouble halfWidth = (getXMax() - getXMin()) / 2.0;
//Cylinder
w1.setSpecies(speciesHandler.getObject(0));
w1.setPosition(Vec3D(mid.X, mid.Y, 0));
w1.setAxis(Vec3D(0, 0, 1));
w1.addObject(Vec3D(1, 0, 0), Vec3D((getXMax() - getXMin()) / 2.0, 0, 0));
wallHandler.copyAndAddObject(w1);
//Cone
w2.setSpecies(speciesHandler.getObject(0));
w2.setPosition(Vec3D(mid.X, mid.Y, 0));
w2.setAxis(Vec3D(0, 0, 1));
std::vector<Vec3D> points(3);
//define the neck as a prism through corners of your prismatic wall in clockwise direction
points[0] = Vec3D(halfWidth, 0.0, mid.Z + contractionHeight);
points[1] = Vec3D(halfWidth - contractionWidth, 0.0, mid.Z);
points[2] = Vec3D(halfWidth, 0.0, mid.Z - contractionHeight);
w2.createOpenPrism(points);
wallHandler.copyAndAddObject(w2);
//Flat surface
w0.setSpecies(speciesHandler.getObject(0));
w0.set(Vec3D(0, 0, -1), Vec3D(0, 0, mid.Z));
wallHandler.copyAndAddObject(w0);
const int N = 600; //maximum number of particles
p0.setSpecies(speciesHandler.getObject(0));
for (Mdouble z = mid.Z + contractionHeight;
particleHandler.getNumberOfObjects() <= N;
z += 2.0 * maxParticleRadius)
{
for (Mdouble r = halfWidth - maxParticleRadius; r > 0; r -= 1.999 * maxParticleRadius)
{
for (Mdouble c = 2.0 * maxParticleRadius; c <= 2 * constants::pi * r; c += 2.0 * maxParticleRadius)
{
p0.setRadius(random.getRandomNumber(minParticleRadius, maxParticleRadius));
p0.setPosition(Vec3D(mid.X + r * mathsFunc::sin(c / r),
mid.Y + r * mathsFunc::cos(c / r),
z + p0.getRadius()));
const Mdouble vz = random.getRandomNumber(-1, 0);
const Mdouble vx = random.getRandomNumber(-1, 1);
p0.setVelocity(Vec3D(vx,0.0,vz));
particleHandler.copyAndAddObject(p0);
}
}
}
}
Use AxisymmetricIntersectionOfWalls to Screw Screw::read Screw::read Screw::read define axisymmetric ...
Definition: AxisymmetricIntersectionOfWalls.h:126
void setAxis(Vec3D a)
Definition: AxisymmetricIntersectionOfWalls.cc:152
void createOpenPrism(std::vector< Vec3D > points, Vec3D prismAxis)
Creates an open prism which is a polygon between the points, except the first and last point,...
Definition: IntersectionOfWalls.cc:467
Mdouble Y
Definition: Vector.h:66
Mdouble Z
Definition: Vector.h:66
Mdouble X
the vector components
Definition: Vector.h:66

First, it is added the cylinder, then the prism shape or cone that represents the outer part. A flat surface is included in the middle to stop the downforward flow until a specific time. After adding the walls, the species of the particles are included. All particles have the same properties and they are randomly distributed.

An additional member function is included.

//Initially, a wall is inserted in the neck of the hourglass to prevent particles flowing through.
//This wall is moved to form the base of the hourglass at time 0.9
void actionsAfterTimeStep() override {
if (getTime() < 0.9 && getTime() + getTimeStep() > 0.9)
{
logger(INFO, "Shifting bottom wall downward");
dynamic_cast<InfiniteWall*>(wallHandler.getLastObject())->set(Vec3D(0, 0, -1), Vec3D(0, 0, getZMin()));
}
}

This overide function moves the flat surface to the base domain. It happens at a defined time, thus the particles can continue to flow.

### Main function:

The main function of this tutorial includes the setup parameters and the tests for a valid collision time and restitution coefficient.

int main(int argc, char *argv[])
{
Mdouble width = 10e-2; // 10cm
Mdouble height = 60e-2; // 60cm
//Point the object HG to class
Tutorial11 HG(width,height);
//Specify particle radius:
//these parameters are needed in setupInitialConditions()
HG.minParticleRadius = 6e-3; // 6mm
HG.maxParticleRadius = 10e-3; //10mm
//Specify the number of particles
//specify how big the wedge of the contraction should be
const Mdouble contractionWidth = 2.5e-2; //2.5cm
const Mdouble contractionHeight = 5e-2; //5cm
HG.contractionWidth = contractionWidth;
HG.contractionHeight = contractionHeight;
//make the species of the particle and wall
species.setDensity(2000);
species.setStiffness(1e5);
species.setDissipation(0.63);
HG.speciesHandler.copyAndAddObject(species);
//test normal forces
const Mdouble minParticleMass = species.getDensity() * 4.0 / 3.0 * constants::pi * mathsFunc::cubic(HG.minParticleRadius);
//Calculates collision time for two copies of a particle of given dissipation_, k, effective mass
logger(INFO, "minParticleMass = %", minParticleMass);
//Calculates collision time for two copies of a particle of given dissipation_, k, effective mass
const Mdouble tc = species.getCollisionTime(minParticleMass);
logger(INFO, "tc = %", tc);
//Calculates restitution coefficient for two copies of given dissipation_, k, effective mass
const Mdouble r = species.getRestitutionCoefficient(minParticleMass);
logger(INFO, "restitution coefficient = %", r);
//time integration parameters
HG.setTimeStep(tc / 10.0);
HG.setTimeMax(5.0);
HG.setSaveCount(500);
HG.solve(argc, argv);
return 0;
}
Mdouble getCollisionTime(Mdouble mass) const
Calculates collision time for two copies of a particle of given disp, k, mass.
Definition: LinearViscoelasticNormalSpecies.cc:137
Mdouble getRestitutionCoefficient(Mdouble mass) const
Calculates restitution coefficient for two copies of given disp, k, mass.
Definition: LinearViscoelasticNormalSpecies.cc:168
Mdouble getDensity() const
Allows density_ to be accessed.
Definition: ParticleSpecies.cc:118
T cubic(const T val)
calculates the cube of a number
Definition: ExtendedMath.h:115

# T12: Creating objects by using Prisms

#### Problem description:

This tutorial explains how to define objects in space with the Prism function. Also, the user will find that you can presecribe the location and motion of the object. The full code can be found in Tutorial12_PrismWalls.cpp

### Headers:

To work with prismatic walls in this tutorial, the user needs to add the following headers:

### Before the main function:

The walls and Prism object are all predefined in the initial conditions.

void setupInitialConditions() override {
p0.setSpecies(speciesHandler.getObject(0));
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.15 * getXMax(), 0.3335 * getYMax(), 0.0));
p0.setVelocity(Vec3D(1.2, 0.2, 0.0));
particleHandler.copyAndAddObject(p0);
// Creating outer walls
w0.setSpecies(speciesHandler.getObject(0));
w0.set(Vec3D(1.0, 0.0, 0.0), Vec3D(getXMax(), 0.0, 0.0));
wallHandler.copyAndAddObject(w0);
w0.set(Vec3D(-1.0, 0.0, 0.0), Vec3D(getXMin(), 0.0, 0.0));
wallHandler.copyAndAddObject(w0);
w0.set(Vec3D(0.0, 1.0, 0.0), Vec3D(0.0, getYMax(), 0.0));
wallHandler.copyAndAddObject(w0);
w0.set(Vec3D(0.0, -1.0, 0.0), Vec3D(0.0, getYMin(), 0.0));
wallHandler.copyAndAddObject(w0);
// Defining the object in het center of the box
// Create an intersection of walls object w1
// Set the species of the wall
w1.setSpecies(speciesHandler.getObject(0));
std::vector<Vec3D> Points(4);
// Define the vertices of the insert and ensure that points are in clockwise order
Points[0] = Vec3D(0.25 * getXMax(), -0.25 * getYMax(), 0.0);
Points[1] = Vec3D(-0.25 * getXMax(), -0.25 * getYMax(), 0.0);
Points[2] = Vec3D(-0.25 * getXMax(), 0.25 * getYMax(), 0.0);
Points[3] = Vec3D(0.25 * getXMax(), 0.25 * getYMax(), 0.0);
// Creating the object from the vector containing points
w1.createPrism(Points);
/* Define the angular velocity of the object (optional).
* Setting the angular velocity of the object results in rotation of the object around
* the origin (0.0,0.0,0.0). If the object is build such that the center of rotation of the object
* coincides with the origin the object will rotate around its own axis. */
w1.setAngularVelocity(Vec3D(0.0,0.0,1.0));
/* Define the translational velocity of the object (optional) */
w1.setVelocity(Vec3D(0.0,0.0,0.0));
/* Set the desired position of center of the object (optional).
* With this command you can position your object (with the set angular velocity)
* at a specific location in the domain.*/
w1.setPosition(Vec3D(0.5 * getXMax(),0.5 * getYMax(),0.0));
// Add the object to wall w1
wallHandler.copyAndAddObject(w1);
}
void setAngularVelocity(const Vec3D &angularVelocity)
set the angular velocity of the BaseInteractble.
Definition: BaseInteractable.cc:360
void createPrism(std::vector< Vec3D > points, Vec3D prismAxis)
Creates an open prism which is a polygon between the points and extends infinitely in the PrismAxis d...
Definition: IntersectionOfWalls.cc:482

The first part in the the initial conditions is defining the particle.

p0.setSpecies(speciesHandler.getObject(0));
p0.setRadius(0.005);
p0.setPosition(Vec3D(0.15 * getXMax(), 0.3335 * getYMax(), 0.0));
p0.setVelocity(Vec3D(1.2, 0.2, 0.0));
particleHandler.copyAndAddObject(p0);

After defining the particle that outer walls are generated by using four finite walls as seen in previous tutorials

To creat the prism we have to do the following: First one creates an instance of IntersectionOfWalls "w1". Next the species are assigned. This line is followed by defining a 3D vector "Points" with space for four vertices. After that the location of the four vertices are defined. Then, by using w1.createPrism(Points) the prism is formed centered at the origin.

// Defining the object in het center of the box
// Create an intersection of walls object w1
// Set the species of the wall
w1.setSpecies(speciesHandler.getObject(0));
std::vector<Vec3D> Points(4);
// Define the vertices of the insert and ensure that points are in clockwise order
Points[0] = Vec3D(0.25 * getXMax(), -0.25 * getYMax(), 0.0);
Points[1] = Vec3D(-0.25 * getXMax(), -0.25 * getYMax(), 0.0);
Points[2] = Vec3D(-0.25 * getXMax(), 0.25 * getYMax(), 0.0);
Points[3] = Vec3D(0.25 * getXMax(), 0.25 * getYMax(), 0.0);
// Creating the object from the vector containing points
w1.createPrism(Points);

Now the object has been formed you can prescribe several properties of the object such as the angular velocity or translational velocity (using setVelocity). Notice that when using setAngularVelocity the object rotates around the origin, so if an object is located at the origin it will rotate around its own axis while located elsewhere it will rotate around the origin. It is essential that you set velocities and positions in the correct order such that you get the desired motion.

w1.setAngularVelocity(Vec3D(0.0,0.0,1.0));
w1.setVelocity(Vec3D(0.0,0.0,0.0));
w1.setPosition(Vec3D(0.5 * getXMax(),0.5 * getYMax(),0.0));

### Main:

In the main function the problem is defined in the problem instance, an instance of class Tutorial 12.

// Problem setup
Tutorial12 problem; // instantiate an object of class Tutorial 12
problem.setName("Tutorial12");
problem.setSystemDimensions(2);
problem.setGravity(Vec3D(0.0, 0.0, 0.0));
problem.setXMax(0.5);
problem.setYMax(0.5);
problem.setTimeMax(5.0);
// The normal spring stiffness and normal dissipation is computed and set as
// For collision time tc=0.005 and restitution coefficeint rc=1.0,
species.setDensity(2500.0); //sets the species type_0 density
species.setStiffness(258.5);//sets the spring stiffness.
species.setDissipation(0.0); //sets the dissipation.
problem.speciesHandler.copyAndAddObject(species);
problem.setSaveCount(10);
problem.setXBallsAdditionalArguments("-solidf -v0 -s .85");
problem.setTimeStep(.005 / 50.0); // (collision time)/50.0
[T12:headers]
Definition: Tutorial12_PrismWalls.cpp:45

At the end of the main function we call the solve function

problem.solve(argc, argv);

After solving the problem you can visualize the results and see the effect of the motion of the object on the trajectory of the particle.