Tutorial on API usage for bond modelling

Description:

This tutorial shows how to implement a bond model using the Aspherix(R) API.

Introduction:

The previous tutorials on how to use the API only used the API to execute Aspherix(R) input files and single commands. However, the API is actually made for giving access to Aspherix(R) itself in order to extract data directly from it or to extend its features. In the present tutorial the Bond model that is also implemented in Aspherix(R) as cohesion model bond is recreated using the API. This will demonstrate how to get access to material parameters from an input script, define a contact model and many more features.

If you have not read the parallel API tutorial then head over to that page and read through it carefully. Everything that follows here assumes that you are familiar with that tutorial case and know how to compile an API executable.

The bond API example explained

The input script bond_external.asx uses the contact model that we are implementing in the following. It defines a total of 4 particle pairs that are interacting with each other. As well as several material parameters that are required from the contact model defined in the API. Compared to a standard Aspherix(R) input script the most significant difference is the line:

particle_contact_model normal external settings external_model bond

The normal model has to be specified as external and in the settings part of the particle_contact_model command it is necessary to specify the name of the contact model, in this case bond. We will soon see the definition of this model in the main C++ file.

The CMakeLists.txt file is similar to the ones seen in the previous tutorials. Compared to them we now have three C++ files which are

bond_model.cpp: Contains main which contains the commands to Aspherix(R)
cohesion_bond.h: Contains the declaration of the contact model
cohesion_bond.cpp: Contains the actual implementation of the contact model

The bond_model.cpp file:

This file is rather short and mainly instantiates the respective class object and calls Aspherix(R) to execute the above mentioned input script bond_external.asx. As such there are not that many new concepts to be discussed.

// contact model include from this tutorial
#include "cohesion_bond.h"

The only new include required here is for cohesion_bond.h so that we can use the newly defined contact model (details on this file below).

// create bond model
CohesionBond *model = new CohesionBond(asx, "bond");

Directly after the definition of the Aspherix object we define an instance of our new contact model with name CohesionBond and its constructor takes two arguments, the first being a pointer to our Aspherix object and the second being the name of the contact model. This name is the same that has to be used in the settings part of the particle_contact_model command inside the input script.

The only other thing to note here is that we need to delete the contact model before the Aspherix object, i.e.

// end simulation
delete model;
delete asx;

otherwise we will run into a segfault.

The cohesion_bond.h file:

This file declares the class CohesionBond which is our contact model. This contact model shall use the bond model if two particles are in the vicinity of each other and in case of an overlap of these two particles, a Hertz plus Tangential history model will be employed additionally.

The file starts out with an include guard:

#ifndef COHESION_BOND_H
#define COHESION_BOND_H

which is followed by several includes.

// includes from the standard library
#include <string>

// aspherix api includes
#include "aspherix.h"
#include "aspherix_contact_model_external.h"
#include "aspherix_global_properties.h"
#include "aspherix_particle.h"
#include "aspherix_particle_interaction.h"
#include "aspherix_vector.h"

Note, compared to the previous tutorials, we now have to include several more aspherix_*.h files as we not only need access to the Aspherix class, but also others such as Particle, ParticleInteraction, GlobalProperties etc. which are all in the Aspherix_API namespace. A list of all header files and their contents can be found in the API section of the documentation.

The class is then declared as

class CohesionBond : public Aspherix_API::ContactModelExternal

where it is important to note that any contact model must be derived as public from the ContactModelExternal class.

After that we have the list of functions that are public to this class. Besides the constructor and destructor (which is empty) we have several functions that are declared as overrides, i.e. they are virtual functions in our parent class ContactModelExternal and are re-implemented in our class.

public:
    CohesionBond(Aspherix_API::Aspherix *asx, std::string id);

    ~CohesionBond()
    {}

    // this member function is responsible for adding history values for the model
    void addHistoryValues() override;

    // this member function is responsible for connecting the contact model to material or global properties
    void registerGlobalProperties() override;

    // this member function determines the maximum distance between two particles that still requires a contact model interaction
    double maximumParticleDistance(const int itype, const int jtype) override;

    // member function to deal with two particles that are overlapping
    void surfacesIntersect(Aspherix_API::Particle &pi, Aspherix_API::Particle &pj, Aspherix_API::ParticleInteraction &pij, const Aspherix_API::GlobalProperties &properties) override;

    // member functions to deal with two particles that are in each others vicinity (at most maximumParticleDistance apart)
    void surfacesClose(Aspherix_API::Particle &pi, Aspherix_API::Particle &pj, Aspherix_API::ParticleInteraction &pij, const Aspherix_API::GlobalProperties &properties) override;

    // member function that is executed directly before the loop over all particle interactions
    void beginPass(const Aspherix_API::GlobalProperties &properties) override;

The detailed implementation of these functions and their functionality is going to be discussed below.

Finally, we have the private section of our class that defines some internal functions that are not available outside of the current class and will be called from within the public functions shown above. This section also defines some private member variables of the class and their use will be described below in the appropriate functions.

private:
    // this function is specific to this case and determines which contact models we call
    void fullModel(Aspherix_API::Particle &pi, Aspherix_API::Particle &pj, Aspherix_API::ParticleInteraction &pij, const Aspherix_API::GlobalProperties &properties, const bool do_hertz);

    void modelHertzHistory(Aspherix_API::Particle &pi, Aspherix_API::Particle &pj, Aspherix_API::ParticleInteraction &pij, const Aspherix_API::GlobalProperties &properties, Aspherix_API::Vector &force, Aspherix_API::Vector &torque_i, Aspherix_API::Vector &torque_j);

    void modelBond(Aspherix_API::Particle &pi, Aspherix_API::Particle &pj, Aspherix_API::ParticleInteraction &pij, const Aspherix_API::GlobalProperties &properties, Aspherix_API::Vector &force, Aspherix_API::Vector &torque_i, Aspherix_API::Vector &torque_j);

    void breakBond(Aspherix_API::ParticleInteraction &pij);

    void createBond(Aspherix_API::ParticleInteraction &pij, const double distance, Aspherix_API::Vector position_i, Aspherix_API::Vector normal);

    // Hertz/History
    size_t Yeff_index_, Geff_index_;
    size_t betaeff_index_, cf_index_;
    // Bond
    long long time_step_create_bond_;
    size_t create_distance_index_, break_distance_index_;
    size_t normal_bond_stiffness_index_, tangential_bond_stiffness_index_;
    size_t limit_sigma_index_, limit_tau_index_;
    size_t normal_force_damping_index_, tangential_force_damping_index_;
    size_t normal_torque_damping_index_, tangential_torque_damping_index_;
    bool stress_break_, damping_;

The cohesion_bond.cpp file:

This file implements all functions that were previously declared in the header file without a definition. In our case this concerns all functions except for the destructor (which is empty). We will now go through each function of this class.

CohesionBond (Constructor)

The constructor takes two arguments as input, the first is a pointer to an Aspherix object and the second is an string id, i.e. the name with which we want to refer to our contact model in the input script.

CohesionBond::CohesionBond(Aspherix *asx, std::string id) :

this is followed by the constructor of the parent class

ContactModelExternal(asx, id),

which takes exactly the same two parameters. The remainder of the initializer list here concerns the member variables. As they will be discussed in greater detail later on, we will skip them for now. Inside of the body we have

registerOnOff("stress_break", stress_break_);
registerOnOff("damping_bond", damping_);

which registers two settings for the contact model. registerOnOff takes two arguments, the first being the name of the setting and the variable (a boolean) which will contain this setting. Since we use these variables throughout the class they are member variables (note, for clarity all member variables have an appended underscore). The stress_break_ variable will be used in the bond model and if it is true the bond will break if a certain stress is exceeded and if it is false then it will break if a certain distance is exceeded (see below). The damping_ variable is also for the bond model and determines whether the bond will be damped or not. In the input script in the particle_contact_model or wall_contact_model command you can use these two settings as shown in this example:

particle_contact_model normal external settings external_model bond stress_break off damping_bond on

Finally, the constructor is finished with

connectToAspherix();

which is mandatory and needs to be at the end of each contact model constructor. It ensures that the contact model is properly connected to Aspherix(R).

addHistoryValues

void CohesionBond::addHistoryValues()

This function allows to add history values to the contact model. History values are per-interaction values that are stored for use in the contact model calculation persistent across time steps. As an example, you can store the current overlap of the particles and retrieve it in the next time step to calculate the difference in overlaps. There are two functions which can be used to assign history values:

addHistoryValue("name_value", 0);
addHistoryVector("name_vector", 0);

The first one is to store a single value, i.e. scalar, with a unique identifying name ‘name_value’ and the second one is for storing a (three-dimensional) vector with a unique name ‘name_vector’. There is no restriction on the naming of your history values except that they have to be unique. The second argument in these functions is responsible for the behavior of the value in case particle order switches. Let us assume that we have a particle p_i and a particle p_j with $i<j$ then each contact model calculates their interaction only once, so a contact model call looks something like contact_model(p_i, p_j). Now a neighbor list rebuild can and will reorder particles. So suddenly we could have $j<i$ and the call would look like contact_model(p_j, p_i). Now each history value $h$ must know whether $h_{ij} = h_{ji}$ or $h_{ij} = -h_{ji}$ so that the contact model is still calculated properly. An example for a history value that keeps its sign is the overlap. It’s always a positive number and it is the same regardless of the order of the particles. Such a history value requires a 0 in the second argument of the addHistoryValue or addHistoryVector function. The relative velocity is an example for a vector that changes its sign when you switch the order of the particles. In this case a 1 is required in the second argument of these functions.

registerGlobalProperties

The function

void CohesionBond::registerGlobalProperties()

is used to request material properties from the input script. An example would be the Youngs Modulus in the contact model Hertz. There are three different functions that can be called from this function:

registerPairProperty register a material interaction property (e.g. the coefficient of restitution)
registerTypeProperty register a material property (e.g. the Youngs Modulus)
registerScalarProperty register a scalar property (e.g. the time step at which the bond should be created)

Each of these functions need to be called with an appropriate identifier string that is also going to be used in the input script. There are two ways of accessing the values of these strings in the actual contact model, the first is to use exactly this string, e.g. getScalarProperty(“my_scalar_property”) or by using the unique property index getScalarProperty(my_scalar_property_index). This index is unique and is returned by the register*Property functions. Note that the second version is much faster and as a contact model is generally performance critical it is advisable to use this second variant. As can be seen in this model all indices are stored in member variables of our CohesionBond class for later access.

Note also the following code block:

if (getRegisterOnOff("stress_break"))
{
    limit_sigma_index_ = registerPairProperty("limit_sigma");
    limit_tau_index_ = registerPairProperty("limit_tau");
}

where we identify whether the user has specified the stress_break setting in the contact model and whether it is set to on. If this is true the additional limit_sigma and limit_tau pair properties are requested.

maximumParticleDistance

The next function that is overridden is

double CohesionBond::maximumParticleDistance(const int itype, const int jtype)

It determines the maximum distance between two particles of type (material) itype and jtype. This distance is measured relative to the particle surface. As we will see later the contact model is split in two parts, the surfacesIntersect function is called if two particles overlap and the surfacesClose function is called if two particles are close by. “Close by” here refers to the fact that their surfaces are less than maximumParticleDistance apart. It is important to set this distance before the start of the simulation as it determines the maximum interaction distance between particles which is used in the building of the neighbor list. Because of this it is important that this value is as small as allowed by the contact model as very large values (e.g. larger than a particle radius) can lead to significant performance reduction.

surfacesIntersect

This is one of the two main contact model functions. It is called for every particle interaction in case the two particles overlap.

void CohesionBond::surfacesIntersect(Particle &pi, Particle &pj, ParticleInteraction &pij, const GlobalProperties &properties)

The arguments of this function are the particle data for particle i (pi) and j (pj), as well as the interaction data pij and the global properties of the simulation (properties). In our case the body of this function is simple as we call the private member function fullModel. Note, the last argument of the called function is true, with which we will indicate to that function that the Hertz model should be calculated as well, since particles are overlapping.

surfacesIntersect

This is the other main contact model function. It is called for every particle interaction in case the two particles do not overlap.

void CohesionBond::surfacesClose(Particle &pi, Particle &pj, ParticleInteraction &pij, const GlobalProperties &properties)

The arguments of this function are the same as for surfacesIntersect and we similarly call the fullModel function. This time the last argument is false and only the bond model will be calculated as we will see below.

beginPass

The function beginPass

void CohesionBond::beginPass(const GlobalProperties &properties)

is called directly before the loop over all particle interactions. It can be used to update variables outside the interior loops for performance improvements. In our current model we get the new settings from the registry and the scalar property time_step_create_bond using the slow method of access via string.

All functions that will be described below are now private member function that are used to implement the Bond and Hertz model.

fullModel

void CohesionBond::fullModel(Particle &pi, Particle &pj, ParticleInteraction &pij, const GlobalProperties &properties, const bool do_hertz)

This function coordinates the use of the two models. If the do_hertz variable is set it computes the Hertz model on top of the bond model, which is always calculated. The force acting on the two particles as well as the torques are passed to each function via reference and in the end the resulting force and torques are assigned to the interaction.

// set the force for both particles (with opposite signs)
pij.setForce(total_force);
// set torque acting on particle i
pij.setTorqueI(total_torque_i);
if (!pij.isWall())
    // set torque acting on particle j, only required if particle j is not a wall (particle i can never be a wall)
    pij.setTorqueJ(total_torque_j);

The setForce function needs to be called only once as the force, due to Newton’s third law, is defined for both particles by one vector. Note, the Torque does not need to be set if “particle j” is a wall element. Particle i can never be a wall element.

modelHertzHistory

The modelHertzHistory is the equivalent model of the normal model Hertz and tangential model history. We are not going to discuss the exact computation as the formulae can be found in the documentation of the models and instead focus on the API functions.

void CohesionBond::modelHertzHistory(Particle &pi, Particle &pj, ParticleInteraction &pij, const GlobalProperties &properties, Vector &force, Vector &torque_i, Vector &torque_j)

At first we are getting some properties of the particles and the interaction

const bool is_wall = pij.isWall();
const double ri = pi.getRadius();
const double rj = is_wall ? 0. : pj.getRadius();

For a detailed list of all available functions check the ParticleInteraction and Particle documentation.

We also need the actual values from the particle properties for the current particle materials:

const double Yeff = getPairProperty(Yeff_index_, itype, jtype);
const double Geff = getPairProperty(Geff_index_, itype, jtype);
const double betaeff = getPairProperty(betaeff_index_, itype, jtype);
const double cf = getPairProperty(cf_index_, itype, jtype);

where we use the indices from the registerGlobalProperties function discussed above.

After calculating the normal Hertz contact we compute the tangential force by getting the old shear strain from the history

Vector shear = pij.getHistoryVector("shear");

we update the value of the shear with the current tangential velocity and finally set the new history value with

pij.setHistoryVector("shear", shear);

For details on how to use the Vector class refer to its documentation.

This model also requires global simulation data, which can be obtained from the GlobalProperties class, e.g.

const double dt = properties.getDt();

In the end the model writes the calculated force and torques into the force, torque_i and torque_j variables that were passed with reference to this function.

modelBond

The modelBond is the equivalent model of the cohesion model bond.

void CohesionBond::modelBond(Particle &pi, Particle &pj, ParticleInteraction &pij, const GlobalProperties &properties, Vector &force, Vector &torque_i, Vector &torque_j)

This function uses the same kind of calculations as the previous function. And we will go briefly through the functionality of the model. The history value bond_exists is a scalar that is 1 if a bond exists and 0 otherwise.

if (bond_exists < 0.5)
{
    const bool create_now = (time_step_create_bond_ < 0 || properties.getTimeStep() == time_step_create_bond_);
    const double distance_create = is_wall ? 2.*distance : distance;
    const bool is_close = distance_create < getPairProperty(create_distance_index_, itype, jtype);
    if (create_now && is_close)
        createBond(pij, distance, pi.getPosition(), pij.getNormal());
    else
        return;
}

In case it does not exist it is created either if the current time step is equal to the time_step_create_bond_ which was set using the global scalar registry. Or, if this property was not set, then a bond is create if the two particles are together close enough (using the create_distance pair property). If these conditions are violated the model will simply exit and keep the force at zero (or whatever the Hertz model computed).

else // bond exists
{
    if (is_wall)
    {
        contact_position += pj.getVelocity()*dt;
    }
}

If the bond exists and the particle-wall interaction is computed the velocity of the wall element (which is equal to pj) is used to update the contact position of the interaction.

Next, the length of the bond is computed and if stress_break is off then the bond is broken if the maximum distance is exceeded:

if (!stress_break_ && distance > getPairProperty(break_distance_index_, itype, jtype))
{
    breakBond(pij);
    return;
}

This is followed by the computation of the relative normal and tangential linear velocities and angular velocities:

Vector vr = pij.getRelativeVelocity();
Vector vn = vr.projectOn(normal);
Vector vt = vr - vn;
Vector weighted_wr = is_wall ? pi.getAngularVelocity()*0.5 : (pi.getAngularVelocity()*ri + pj.getAngularVelocity()*rj)*radsuminv;
Vector vtr = vt + normal.cross(weighted_wr)*distance;

Vector wr = pij.getRelativeAngularVelocity();
Vector wn = wr.projectOn(normal);
Vector wt = wr - wn;

With all the established parameters we can then proceed to compute the normal and tangential forces and torques, e.g. here the code for the normal force:

Vector normal_force_damped;
Vector normal_force;
if (displacement < -1e-15 || displacement > 1e-15)
{
    normal_force = normal*(kn*A*displacement);
    normal_force_damped = normal_force;
    if (damping_)
        normal_force_damped -= normal_force.abs().Hadamard(vn.sgn())*getPairProperty(normal_force_damping_index_, itype, jtype);
}

Note the heavy use of Vector functions such as abs, Hadamard and sgn which give the vector of the absolute value of the components, execute a component-wise product and return the vector of the sign of each component.

Now that we have computed the forces we can compute the stresses and in case stress_break is on we can determine whether the bond should break:

if (stress_break_)
{
    const double current_sigma = normal_force.length()/A + tangential_torque.length()*rb/I;
    const bool normal_stress_exceeded = getPairProperty(limit_sigma_index_, itype, jtype) < current_sigma;
    const double current_tau = tangential_force.length()/A + normal_torque.length()*rb/J;
    const bool tangential_stress_exceeded = getPairProperty(limit_tau_index_, itype, jtype) < current_tau;
    if (normal_stress_exceeded || tangential_stress_exceeded)
    {
        breakBond(pij);
        return;
    }
}

If the bond is not broken at this stage we set the final forces and torques and update some history values at the end of the function:

Vector torque = tangential_force_damped.cross(normal);

pij.setHistoryVector("tangential_force", tangential_force);
pij.setHistoryVector("normal_torque", normal_torque);
pij.setHistoryVector("tangential_torque", tangential_torque);

force += normal_force_damped + tangential_force_damped;
torque_i += torque*cri + normal_torque_damped + tangential_torque_damped;
torque_j += torque*crj - normal_torque_damped - tangential_torque_damped;

breakBond

The breakBond function simply breaks the bond by resetting the bond_exists and initial_distance history values.

createBond

The createBond function initiates a bond by setting several history values such as bond_exists and initial_distance. The Vector() objects simply represent the 3-D zero vector.

Running the parallel API example

Linux:

After the compilation which is equivalent to the one shown in our previous tutorial case we can verify that the executable is where we expect it to be:

$ ls install/bin/
bond_model

To execute the code we now use mpirun with two processors:

$ mpirun -np 2 install/bin/bond_model

and watch the simulation output:

Aspherix (Version Aspherix 7.1.0, compiled 2026-03-12-13:20:23 by vagrant, git commit 3efa378ae0fe305d4e235e2e5c66d2244884a712)
Checkout of asx_solver OK.
Checkout of asx_solver OK.
Created orthogonal box = (-0.25 -0.5 -0.15) to (0.25 0.5 0.15)
  1 by 2 by 1 MPI processor grid

[... more output from the Aspherix(R) simulation ...]

Total # of neighbors = 12
Ave neighs/atom = 1.5
Neighbor list builds = 945
Dangerous builds = 0

 >> Simulation produced 1 warnings. See the warnings file (default: warnings_aspherix.txt) for details. <<

In the post folder you can find the aspherix-simulation.pvd that can be opened with ParaView. If you use the ExtractBlock filter to visualize the particles and the contact network you should be able to observe the four pairs of bonded particles interacting with each other.