Tutorial on advanced API usage

Description:

This tutorial explains the demo API case showing many different functionalities.

Introduction:

This tutorial does not actually produce a really sensible simulation. Instead it severes as a demonstration platform for a wide variety of API features. As such we expect the user to be familiar with the more basic API tutorials and we will skim over the features that were already explained in-depth in those.

The demo API example explained

The example consists of an input script demo.asx with an associated mesh to be found in meshes/plate.stl. The input script creates three particles at defined locations and initial velocities and adds the mesh as solid object to the simulation. As mentioned above, we do not really concern ourselves with the physical behavior of this simulation, so we’ll leave it at that.

There are a total of three C++ file sets that are of interest here:

1. demo.cpp: This contains the main function with the simulation control

2. my_contact_model.cpp/h: Implements a basic contact model

3. my_fix_external.cpp/h: Implements a basic Fix

The CMakeLists.txt file is similar to the ones seen in the previous tutorials so we won’t discuss it further.

The demo.cpp file:

This file controls the flow of the simulation. It begins by starting the simulation using the demo.asx case:

asx.executeInputFile("demo.asx");

Next we get the full list of our three particles and print their properties:

auto plist = asx.getParticleList();
for (const auto &p : plist)
{
    std::cout << "Particle with id " << p.getId() << " has type " << p.getType() << " and mass " << p.getMass() << std::endl;
}
std::cout << "Alternative access of id 2 mass: " << plist[2]->getMass() << std::endl;

Note, in the last line we use the []-operator to access a particle directly using its unique id.

Next, some general properties of the simulation are written out:

std::cout << "Timestep of simulation is " << asx.getGlobalProperties().getDt() << " s" << std::endl;
std::cout << "Simulation currently has " << asx.getGlobalProperties().getNumberOfParticles() << " particles" << std::endl;

In the next lines we define our new fix MyFixExternal that is defined in the my_fix_external.cpp/h files, which are described in detail further below. Additionally, we let Aspherix(R) know that this fix has two output values that can be accessed either via the input script or via the API (see below).

MyFixExternal fix(&asx, "my_fix_external");
fix.hasVectorOutput(2);

asx.executeCommand("status_style custom step time atoms ke f_my_fix_external[1] f_my_fix_external[2]");

The last line modifies the status output that is written to the terminal, log file and the simulation data file, with the new output generated by this fix.

asx.executeCommand("fix fpas all property/atom fpas scalar no no no 1.");
asx.executeCommand("fix fpav all property/atom fpav vector no no no 2. 3. 4.");

Two new fixes of style property/atom are added in order to demonstrate some of the API features below. Note, normally you interact with property/atom for example when using the heat conduction models and access the particle temperature (via Temp) and you will not add them in your input script yourself.

// alternative way of accessing output
auto myfix = asx.getFix("my_fix_external");
std::cout << "Fix my_fix_external computes the following two output values: " << myfix.getVector(0) << ", " << myfix.getVector(1) << std::endl;

As an alternative way to access the output of the fix we just created, we can use the getFix function to get a reference to our fix and access its calculated values via getVector.

The API is also capable of interacting with variables of style atom and equal as shown in the following code:

auto var_rho = asx.getVariable("rho");
std::cout << "Variable rho for atom with id 1 = " << var_rho.getAtomValue(1) << std::endl;
auto var_time = asx.getVariable("time");
std::cout << "Time from variable = " << var_time.getValue() << std::endl;

Next, it is demonstrated how to access the current values of the fixes of style property/atom defined above. Note, the access of fpav with getFixPropertyScalar is designed to fail as it is a vector and not a scalar, so the message will be displayed.

std::cout << "Access fix property/atom fpas, scalar of id 1: " << plist[1]->getFixPropertyScalar("fpas") << std::endl;
try
{
    double scalar = plist[2]->getFixPropertyScalar("fpav");
}
catch(std::exception &e)
{
    std::cout << "fpav fix is not a scalar" << std::endl;
}
std::vector<double> vector = plist[2]->getFixPropertyVector("fpav");
std::cout << "Access fix property/atom fpav, vector of id 2: ";
std::for_each(vector.begin(), vector.end(), [](double v){ std::cout << v << " "; });
std::cout << std::endl;

The subsequent code block demonstrates the access of a mesh object using the API.

Mesh mesh;
try
{
    auto force = mesh.getForce();
}
catch(std::exception &e)
{
    std::cout << "Uninitialized mesh as expected" << std::endl;
}
try
{
    mesh = asx.getMesh("mesh_wrong_id");
}
catch(std::exception &e)
{
    std::cout << "Could not find mesh with mesh_wrong_id as expected" << std::endl;
}
try
{
    mesh = asx.getMesh("mesh1");
}
catch (std::exception &e)
{
    std::cout << "This should not be triggered" << std::endl;
}
auto force = mesh.getForce();
std::cout << "Force acting on 'mesh1' is: " << force.as_string() << std::endl;

The getForce function returns the total force acting on the mesh as a Vector object, which is displayed using the as_string function.

We then inspect the Mesh object by listing its scalar and vector properties and displaying one of each.

auto mesh_property_list_scalar = mesh.getListOfScalarProperties();
std::cout << "List of scalar properties: ";
for (auto prop : mesh_property_list_scalar)
    std::cout << prop << ", ";
std::cout << std::endl;
double mesh_area = mesh.getScalarProperty("areaMesh");
std::cout << "Area of mesh is: " << mesh_area << std::endl;
auto mesh_property_list_vector = mesh.getListOfVectorProperties();
std::cout << "List of vector properties: ";
for (auto prop : mesh_property_list_vector)
    std::cout << prop << ", ";
std::cout << std::endl;
Vector ref_point = mesh.getVectorProperty("p_ref");
std::cout << "Reference point of mesh is: " << ref_point.as_string() << std::endl;

The final objective of the Mesh demo object is the access of properties of single mesh elements which is started with a loop over all elements. We then display the id and z-force of each element

auto meshElementList = mesh.getMeshElementList();
for (auto meshElement = meshElementList.begin(); meshElement != meshElementList.end(); meshElement++)
{
    std::cout << "Mesh element with id " << meshElement->getId() << " has z-force: " << meshElement->getForce()[2] << std::endl;
    [...]
}

We then use the first element and, similar to the global mesh properties, inspect the per-element mesh properties. For this first element we also demonstrate how to alter these properties, once for a scalar and once for a vector.

if (meshElement == meshElementList.begin())
{
    auto property_list_scalar = meshElement->getListOfScalarProperties();
    std::cout << "List of scalar properties: ";
    for (auto prop : property_list_scalar)
        std::cout << prop << ", ";
    std::cout << std::endl;
    auto property_list_vector = meshElement->getListOfVectorProperties();
    std::cout << "List of vector properties: ";
    for (auto prop : property_list_vector)
        std::cout << prop << ", ";
    std::cout << std::endl;
    double area = meshElement->getScalarProperty("area");
    std::cout << "Area of element 0 is: " << area << std::endl;
    meshElement->setScalarProperty("area", area*2.);
    area = meshElement->getScalarProperty("area");
    std::cout << "Area of element 0 is now twice as large: " << area << std::endl;
    Vector normal = meshElement->getVectorProperty("surfaceNorm");
    std::cout << "Surface normal of element 0 is: " << normal.as_string() << std::endl;
    normal *= -1.;
    meshElement->setVectorProperty("surfaceNorm", normal);
    normal = meshElement->getVectorProperty("surfaceNorm");
    std::cout << "Surface normal of element 0 is now flipped: " << normal.as_string() << std::endl;
}

Finally, for all mesh elements we move all three nodes of each triangle by a small distance in x-direction, similar to what a linear mesh movement would do. Just a small word of caution here, if you wish to implement more complicated movements here you might need to alter also additional properties such as the reference point.

for (int i = 0; i < 3; i++)
    std::cout << "Node " << i << " has position: " << meshElement->getNode(i).as_string() << std::endl;
meshElement->move(Vector(0.01, 0., 0.));
for (int i = 0; i < 3; i++)
    std::cout << "Node " << i << " has new position: " << meshElement->getNode(i).as_string() << std::endl;

In the last step of the main function we instantiate the MyContactModel class that will be described further below and use it to model the wall interaction.

MyContactModel model(&asx, "my_contact_model");

asx.executeCommand("disable_command id wall_contact_model");
asx.executeCommand("wall_contact_model normal external settings external_model my_contact_model my_string hooke_mass my_double 0.9");
asx.executeCommand("fix sp all property/global myScalarProperty scalar 2.0");
asx.executeCommand("fix tp all property/global myTypeProperty peratomtype 0.5");
asx.executeCommand("fix pp all property/global myPairProperty peratomtypepair 1 1e5");

The last three lines here add new material properties to the simulation as they are required by our new contact model.

The my_fix_external.h file:

This file defines a class MyFixExternal that derives from Aspherix_API::FixExternal and which has only one task, which is to compute an output vector that has two components. External fixes can actually do more than that as documented here.

The my_fix_external.cpp file:

This file implements the compute_vector function of the MyFixExternalClass, using a pointer to the Aspherix object and getting some global properties

if (n == 0)
    return aspherix()->getGlobalProperties().getTime();
else if (n == 1)
    return static_cast<double>(aspherix()->getGlobalProperties().getTimeStep());

the variable n represents which component of the output vector is requested. Note, internally the counting starts with zero, but outside you would access the the first element with id_myFixExternalID[1].

The my_contact_model.h file:

This file defines the class for our new contact model. It defines a few public functions that are overridden from our parent class Aspherix_API::ContactModelExternal. For details about these different function templates please refer to the tutorial with the bond model which explains these in great detail.

The my_contact_model.cpp file:

This file implements a very basic hooke model as can be seen in the following line

force = pij.getNormal() * pij.getDelta() * factor;

or an alternate version in case the my_string setting is set:

else if (my_string.compare("hooke_mass") == 0)
    force = pij.getNormal() * pij.getDelta() * factor * (pi.getMass() + pj.getMass())*0.5;

as mentioned before, this is not necessarily a physical example and mainly used to demonstrate possible ways of accessing different types of data. For further details and a proper physically sensible example please refer to the tutorial about the bond case.

The stiffness is here clobbered together using different ways of accessing scalar, material and material interaction properties:

double factor = getScalarProperty("myScalarProperty");
factor *= (pi.getTypeProperty("myTypeProperty")+pj.getTypeProperty("myTypeProperty"))*0.5;
factor *= pij.getPairProperty("myPairProperty");

as well as accessing the settings part of the contact model definitions

force *= getRegisterDouble("my_double");

and finally to demonstrate the use of history values, a scalar is defined, retrieved and set:

// definition:
addHistoryValue("oldForceMag", 0);
[...]
// retrieval:
if (getRegisterOnOff("my_flag") && pij.getHistoryValue("oldForceMag") > 1.0)
//                                 ^^^^^^^^^^^^^^^^^^^
[...]
// setting:
pij.setHistoryValue("oldForceMag", force.length());

Running the parallel API example

Linux:

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

$ ls install/bin/
demo

To execute the code we now type:

$ install/bin/demo

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.
Created orthogonal box = (-0.1 -0.1 -1) to (1 1 1)
  1 by 1 by 1 MPI processor grid

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

Particle with id 1 has type 1 and mass 1.30116
Particle with id 2 has type 1 and mass 1.30116
Particle with id 3 has type 1 and mass 1.30116
Alternative access of id 2 mass: 1.30116
Timestep of simulation is 1e-05 s
Simulation currently has 3 particles

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

Fix my_fix_external computes the following two output values: 0.0002, 20
Variable rho for atom with id 1 = 2500
Time from variable = 0.0002
Access fix property/atom fpas, scalar of id 1: 1
fpav fix is not a scalar
Access fix property/atom fpav, vector of id 2: 2 3 4
Uninitialized mesh as expected
Could not find mesh with mesh_wrong_id as expected
Force acting on 'mesh1' is: (0.000000, 0.000000, -9.189875)
List of scalar properties: areaMesh,
Area of mesh is: 0.25
List of vector properties: p_ref,
Reference point of mesh is: (0.000000, 0.000000, 0.000000)
Mesh element with id 0 has z-force: -9.18987
List of scalar properties: sigma_t, area, areaAcc, sigma_n,
List of vector properties: edgeLen, surfaceNorm, f,
Area of element 0 is: 0.125
Area of element 0 is now twice as large: 0.25
Surface normal of element 0 is: (-0.000000, 0.000000, 1.000000)
Surface normal of element 0 is now flipped: (0.000000, -0.000000, -1.000000)
Node 0 has position: (0.000000, 0.500000, -0.050000)
Node 1 has position: (0.500000, 0.000000, -0.050000)
Node 2 has position: (0.500000, 0.500000, -0.050000)
Node 0 has new position: (0.010000, 0.500000, -0.050000)
Node 1 has new position: (0.510000, 0.000000, -0.050000)
Node 2 has new position: (0.510000, 0.500000, -0.050000)
Mesh element with id 1 has z-force: 0
Node 0 has position: (0.000000, 0.500000, -0.050000)
Node 1 has position: (0.500000, 0.000000, -0.050000)
Node 2 has position: (0.000000, 0.000000, -0.050000)
Node 0 has new position: (0.010000, 0.500000, -0.050000)
Node 1 has new position: (0.510000, 0.000000, -0.050000)
Node 2 has new position: (0.010000, 0.000000, -0.050000)

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

Ave neighs/atom = 0.333333
Neighbor list builds = 12
Dangerous builds = 0

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

Observe the output that comes from our little demo case in between the simulations and verify that they produce the correct changes, such as the movement of the mesh.