Apply cosimulated loads to amesh (demo_FEA_cosimulate_load.cpp)

Tutorial that teaches how to use the FEA module to import an .INP Abaqus mesh with a 3D tetrahedral mesh and apply loads to the surface coming from an external process.

The external process here is simply simulated using a function in the same .cpp unit, but in a context of cosimulation it could be an external process/program like, say, a CFD software, that receives the 3D mesh from Chrono and gives back the fluid forces to Chrono.

// =============================================================================
// Copyright (c) 2014
// All rights reserved.
// Use of this source code is governed by a BSD-style license that can be found
// in the LICENSE file at the top level of the distribution and at
// =============================================================================
// Authors: Alessandro Tasora
// =============================================================================
// FEA force-displacement co-simulation
// - loading an Abaqus tetrahedron mesh
// - apply a load to the mesh using an external tool, e.g. CFD.
// (here simulated as a function in this .cpp file)
// =============================================================================
#include "chrono/geometry/ChTriangleMeshConnected.h"
#include "chrono/physics/ChLoadBodyMesh.h"
#include "chrono/physics/ChLoadContainer.h"
#include "chrono/physics/ChLoaderUV.h"
#include "chrono/physics/ChSystemSMC.h"
#include "chrono/solver/ChIterativeSolverLS.h"
#include "chrono/fea/ChElementTetraCorot_4.h"
#include "chrono/fea/ChLoadContactSurfaceMesh.h"
#include "chrono/fea/ChMesh.h"
#include "chrono/fea/ChMeshFileLoader.h"
#include "chrono/assets/ChVisualShapeFEA.h"
#include "chrono_irrlicht/ChVisualSystemIrrlicht.h"
using namespace chrono;
using namespace chrono::fea;
using namespace chrono::irrlicht;
// This function simulates the effect of an external program that
// gets a triangle mesh and outputs the forces acting on the nodes.
// In a real cosimulation scenario, this procedure could even reside
// on a different computing node and manage inputs/outputs via MPI or such.
void PerformExternalCosimulation(const std::vector<ChVector<>>& input_vert_pos,
const std::vector<ChVector<>>& input_vert_vel,
const std::vector<ChVector<int>>& input_triangles,
std::vector<ChVector<>>& vert_output_forces,
std::vector<int>& vert_output_indexes) {
double ky = 10000; // upward stiffness
double ry = 20; // upward damping
// simple example: scan through all vertexes in the mesh, see if they sink below zero,
// apply a penalty upward spring force if so.
for (int iv = 0; iv < input_vert_pos.size(); ++iv) {
if (input_vert_pos[iv].y() < 0) {
double yforce = -ky * input_vert_pos[iv].y() - ry * input_vert_vel[iv].y();
if (yforce > 0) {
vert_output_forces.push_back(ChVector<>(0, yforce, 0));
// note that:
// - we avoided adding forces to vert_output_forces when force was zero.
// - vert_output_forces has the same size of vert_output_indexes, maybe smaller than input_vert_pos
// Utility to draw some triangles that are affected by cosimulation.
// Also plot forces as vectors.
// Mostly for debugging.
void draw_affected_triangles(ChVisualSystemIrrlicht& vis,
std::vector<ChVector<>>& vert_pos,
std::vector<ChVector<int>>& triangles,
std::vector<int>& vert_indexes,
std::vector<ChVector<>>& vert_forces,
double forcescale = 0.01) {
for (int it = 0; it < triangles.size(); ++it) {
bool vert_hit = false;
for (int io = 0; io < vert_indexes.size(); ++io) {
if (triangles[it].x() == vert_indexes[io] || triangles[it].y() == vert_indexes[io] ||
triangles[it].z() == vert_indexes[io])
vert_hit = true;
if (vert_hit == true) {
std::vector<chrono::ChVector<>> fourpoints = {vert_pos[triangles[it].x()], vert_pos[triangles[it].y()],
vert_pos[triangles[it].z()], vert_pos[triangles[it].x()]};
tools::drawPolyline(&vis, fourpoints, ChColor(0.94f, 0.78f, 0.00f), true);
if (forcescale > 0)
for (int io = 0; io < vert_indexes.size(); ++io) {
std::vector<chrono::ChVector<>> forceline = {vert_pos[vert_indexes[io]],
vert_pos[vert_indexes[io]] + vert_forces[io] * forcescale};
tools::drawPolyline(&vis, forceline, ChColor(0.94f, 0.00f, 0.00f), true);
int main(int argc, char* argv[]) {
GetLog() << "Copyright (c) 2017\nChrono version: " << CHRONO_VERSION << "\n\n";
// Global parameter for tire:
double tire_rad = 0.8;
ChVector<> tire_center(0, 0.02 + tire_rad, 0.5);
ChMatrix33<> tire_alignment(Q_from_AngAxis(CH_C_PI, VECT_Y)); // create rotated 180� on y
// Create a Chrono physical system and the associated collision system
// Create the surface material
auto mysurfmaterial = chrono_types::make_shared<ChMaterialSurfaceSMC>();
// Create a mesh, that is a container for groups of FEA elements and their referenced nodes.
auto my_mesh = chrono_types::make_shared<ChMesh>();
// Create a material, that must be assigned to each solid element in the mesh,
// and set its parameters
auto mmaterial = chrono_types::make_shared<ChContinuumElastic>();
mmaterial->Set_E(0.003e9); // rubber 0.01e9, steel 200e9
// Load an ABAQUS .INP tetrahedron mesh file from disk, defining a tetrahedron mesh.
// Note that not all features of INP files are supported. Also, quadratic tetrahedrons are promoted to linear.
// This is much easier than creating all nodes and elements via C++ programming.
// Ex. you can generate these .INP files using Abaqus or exporting from the SolidWorks simulation tool.
std::map<std::string, std::vector<std::shared_ptr<ChNodeFEAbase>>> node_sets;
try {
mmaterial, node_sets, tire_center, tire_alignment);
} catch (const ChException& myerr) {
GetLog() << myerr.what();
return 0;
// Create the contact surface(s).
// Use the AddFacesFromBoundary() to select automatically the outer skin of the tetrahedron mesh.
// Note that the contact material specified here is not used, as contacts will be emulated through cosimulation.
auto mcontactsurf = chrono_types::make_shared<ChContactSurfaceMesh>(mysurfmaterial);
// Create a mesh load for cosimulation, acting on the contact surface above
// (forces on nodes will be computed by an external procedure)
auto mloadcontainer = chrono_types::make_shared<ChLoadContainer>();
auto mmeshload = chrono_types::make_shared<ChLoadContactSurfaceMesh>(mcontactsurf);
// Optional... visualization
// Visualization of the FEM mesh.
// This will automatically update a triangle mesh (a ChVisualShapeTriangleMesh
// asset that is internally managed) by setting proper
// coordinates and vertex colors as in the FEM elements.
auto mvisualizemesh = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemesh->SetColorscaleMinMax(0.0, 10);
// Create also a rigid body with a rigid mesh that will be used for the cosimulation,
// this time the ChLoadContactSurfaceMesh cannot be used as in the FEA case, so we
// will use the ChLoadBodyMesh class:
auto mrigidbody = chrono_types::make_shared<ChBody>();
mrigidbody->SetInertiaXX(ChVector<>(20, 20, 20));
mrigidbody->SetPos(tire_center + ChVector<>(-1, 0, 0));
mesh->Transform(VNULL, Q_from_AngAxis(CH_C_PI, VECT_Y));
auto mesh_asset = chrono_types::make_shared<ChVisualShapeTriangleMesh>();
// this is used to use the mesh in cosimulation!
auto mrigidmeshload = chrono_types::make_shared<ChLoadBodyMesh>(mrigidbody, *mesh);
// Create the Irrlicht visualization system
auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
vis->SetWindowSize(1280, 720);
vis->AddLightWithShadow(ChVector<>(1.5, 5.5, -2.5), ChVector<>(0, 0, 0), 3, 2.2, 7.2, 40, 512,
ChColor(0.8f, 0.8f, 1.0f));
vis->AddCamera(ChVector<>(1.0, 1.4, -1.2), ChVector<>(0, tire_rad, 0));
// Change solver to embedded MINRES
// NOTE! it is strongly advised that you compile the optional MKL module
// if you need higher precision, and switch to its MKL solver - see demos for FEA & MKL.
auto solver = chrono_types::make_shared<ChSolverMINRES>();
solver->EnableWarmStart(true); // Enable for better convergence if using Euler implicit linearized
// Change type of integrator:
sys.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED); // fast, less precise
while (vis->Run()) {
// -------------------------------------------------------------------------
// Here do the cosimulation
// A <--> B
// For example, A is this main program, and B can be an external program,
// e.g. a CFD simulation tool.
// The idea is that A --> B communicates the mesh position,
// then A <-- B receives the computed forces to be applied at nodes.
// In this example, to keep things simple, B is just a simple local function
std::vector<ChVector<>> vert_pos;
std::vector<ChVector<>> vert_vel;
std::vector<ChVector<int>> triangles;
std::vector<ChVector<>> vert_forces;
std::vector<int> vert_indexes;
mmeshload->OutputSimpleMesh(vert_pos, vert_vel, triangles);
PerformExternalCosimulation(vert_pos, vert_vel, triangles, vert_forces, vert_indexes);
mmeshload->InputSimpleForces(vert_forces, vert_indexes);
// now, just for debugging and some fun, draw some triangles
// (only those that have a vertex that has a force applied):
draw_affected_triangles(*vis, vert_pos, triangles, vert_indexes, vert_forces, 0.01);
// Other example: call another cosimulation, this time for the rigid body
// mesh (the second tire, the rigid one):
mrigidmeshload->OutputSimpleMesh(vert_pos, vert_vel, triangles);
PerformExternalCosimulation(vert_pos, vert_vel, triangles, vert_forces, vert_indexes);
mrigidmeshload->InputSimpleForces(vert_forces, vert_indexes);
// now, just for debugging and some fun, draw some triangles
// (only those that have a vertex that has a force applied):
draw_affected_triangles(*vis, vert_pos, triangles, vert_indexes, vert_forces, 0.01);
// End of cosimulation block
// -------------------------------------------------------------------------
tools::drawGrid(vis.get(), 0.1, 0.1, 20, 20, ChCoordsys<>(VNULL, CH_C_PI_2, VECT_X),
ChColor(0.40f, 0.40f, 0.40f), true);
return 0;
void AddTypicalLights()
Simple shortcut to set two point lights in the scene.
Definition: ChVisualSystemIrrlicht.cpp:337
std::string GetChronoDataFile(const std::string &filename)
Obtain the complete path to the specified filename, given relative to the Chrono data directory (thre...
Definition: ChGlobal.cpp:95
void Add(std::shared_ptr< ChPhysicsItem > item)
Attach an arbitrary ChPhysicsItem (e.g.
Definition: ChSystem.cpp:210
virtual void Render() override
Draw all 3D shapes and GUI elements at the current frame.
Definition: ChVisualSystemIrrlicht.cpp:623
Definition: ChCoordsys.h:38
void AddSkyBox(const std::string &texture_dir=GetChronoDataFile("skybox/"))
Add a sky box in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:350
void SetSolverForceTolerance(double tolerance)
Set a solver tolerance threshold at force level (default: not specified).
Definition: ChSystem.h:193
Class for exceptions for throw() catch() mechanism.
Definition: ChException.h:24
virtual void Initialize()
Initialize the visualization system.
Definition: ChVisualSystemIrrlicht.cpp:181
ChLog & GetLog()
Global function to get the current ChLog object.
Definition: ChLog.cpp:39
virtual int AddCamera(const ChVector<> &pos, ChVector<> targ=VNULL) override
Add a camera in an Irrlicht 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:287
Definition of a 3x3 fixed size matrix to represent 3D rotations and inertia tensors.
Definition: ChMatrix33.h:31
static std::shared_ptr< ChTriangleMeshConnected > CreateFromWavefrontFile(const std::string &filename, bool load_normals=true, bool load_uv=false)
Create and return a ChTriangleMeshConnected from a Wavefront OBJ file.
Definition: ChTriangleMeshConnected.cpp:259
ChQuaternion< double > Q_from_AngAxis(double angle, const ChVector< double > &axis)
Get the quaternion from an angle of rotation and an axis, defined in abs coords.
Definition: ChQuaternion.cpp:99
Namespace with classes for the Irrlicht module.
Definition: ChApiIrr.h:47
virtual void EndScene() override
End the scene draw at the end of each animation frame.
Definition: ChVisualSystemIrrlicht.cpp:613
Definition of a visual color.
Definition: ChColor.h:30
virtual void BeginScene() override
Perform any necessary operations at the beginning of each rendering frame.
Definition: ChVisualSystemIrrlicht.cpp:564
Definition of general purpose 3d vector variables, such as points in 3D.
Definition: ChVector.h:35
irr::scene::ILightSceneNode * AddLightWithShadow(const ChVector<> &pos, const ChVector<> &aim, double radius, double near_value, double far_value, double angle, unsigned int resolution=512, ChColor color=ChColor(1, 1, 1), bool directional=false, bool clipborder=true)
Add a point light that cast shadow (using soft shadows/shadow maps) Note that the quality of the shad...
Definition: ChVisualSystemIrrlicht.cpp:404
int DoStepDynamics(double step_size)
Advances the dynamical simulation for a single step, of length step_size.
Definition: ChSystem.cpp:1590
Class for a physical system in which contact is modeled using a smooth (penalty-based) method.
Definition: ChSystemSMC.h:30
Bullet-based collision detection system.
void EnableShadows(std::shared_ptr< ChPhysicsItem > item=nullptr)
Enable shadow maps for all visual models in a scene or only for a single physics item.
Definition: ChVisualSystemIrrlicht.cpp:468
void SetWindowTitle(const std::string &win_title)
Set the windoiw title (default "").
Definition: ChVisualSystemIrrlicht.cpp:138
static void FromAbaqusFile(std::shared_ptr< ChMesh > mesh, const char *filename, std::shared_ptr< ChContinuumMaterial > my_material, std::map< std::string, std::vector< std::shared_ptr< ChNodeFEAbase > > > &node_sets, ChVector<> pos_transform=VNULL, ChMatrix33<> rot_transform=ChMatrix33<>(1), bool discard_unused_nodes=true)
Load tetrahedrons, if any, saved in a .inp file for Abaqus.
Definition: ChMeshFileLoader.cpp:196
virtual void AttachSystem(ChSystem *sys) override
Attach another Chrono system to the run-time visualization system.
Definition: ChVisualSystemIrrlicht.cpp:165
virtual void SetSolver(std::shared_ptr< ChSolver > newsolver)
Attach a solver (derived from ChSolver) for use by this system.
Definition: ChSystem.cpp:370
virtual bool Run() override
Run the Irrlicht device.
Definition: ChVisualSystemIrrlicht.cpp:240
Main namespace for the Chrono package.
Definition: ChCamera.cpp:17
void SetTimestepperType(ChTimestepper::Type type)
Set the method for time integration (time stepper type).
Definition: ChSystem.cpp:463
void AddLogo(const std::string &logo_filename=GetChronoDataFile("logo_chronoengine_alpha.png"))
Add a logo in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:329
virtual void SetCollisionSystemType(ChCollisionSystem::Type type)
Set the collision detection system used by this Chrono system to the specified type.
Definition: ChSystem.cpp:375
Namespace for FEA classes.
Definition: ChVisualShapeFEA.h:28
void SetWindowSize(unsigned int width, unsigned int height)
Set the window size (default 640x480).
Definition: ChVisualSystemIrrlicht.cpp:134
Irrlicht-based Chrono run-time visualization system.
Definition: ChVisualSystemIrrlicht.h:51