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.

// =============================================================================
// PROJECT CHRONO - http://projectchrono.org
//
// Copyright (c) 2014 projectchrono.org
// 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
// http://projectchrono.org/license-chrono.txt.
//
// =============================================================================
// Authors: Alessandro Tasora
// =============================================================================
//
// FEA advanced demo:
// - loading an Abaqus tetrahedron mesh
// - apply a load to the mesh using an external tool,
// say CFD or SPH (here simulated as a function in this .cpp file)
// that is perform a cosimulation.
//
// =============================================================================
#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;
using namespace irr;
// 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) {
vert_output_forces.clear();
vert_output_indexes.clear();
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));
vert_output_indexes.push_back(iv);
}
}
}
// 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 projectchrono.org\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::Engine physical system
//
// CREATE A FINITE ELEMENT MESH
//
// Create the surface material, containing information
// about friction etc.
auto mysurfmaterial = chrono_types::make_shared<ChMaterialSurfaceSMC>();
mysurfmaterial->SetYoungModulus(10e4);
mysurfmaterial->SetFriction(0.3f);
mysurfmaterial->SetRestitution(0.2f);
mysurfmaterial->SetAdhesion(0);
// Create a mesh, that is a container for groups
// of FEA elements and their referenced nodes.
auto my_mesh = chrono_types::make_shared<ChMesh>();
sys.Add(my_mesh);
// 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
mmaterial->Set_v(0.4);
mmaterial->Set_RayleighDampingK(0.004);
mmaterial->Set_density(1000);
// 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 {
GetChronoDataFile("models/tractor_wheel/tractor_wheel_coarse.INP").c_str(),
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);
my_mesh->AddContactSurface(mcontactsurf);
mcontactsurf->AddFacesFromBoundary();
// 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>();
sys.Add(mloadcontainer);
auto mmeshload = chrono_types::make_shared<ChLoadContactSurfaceMesh>(mcontactsurf);
mloadcontainer->Add(mmeshload);
//
// Optional... visualization
//
// Visualization of the FEM mesh.
// This will automatically update a triangle mesh (a ChTriangleMeshShape
// 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->SetFEMdataType(ChVisualShapeFEA::DataType::NODE_SPEED_NORM);
mvisualizemesh->SetColorscaleMinMax(0.0, 10);
mvisualizemesh->SetSmoothFaces(true);
my_mesh->AddVisualShapeFEA(mvisualizemesh);
//
// CREATE A RIGID BODY WITH A MESH
//
// 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>();
sys.Add(mrigidbody);
mrigidbody->SetMass(200);
mrigidbody->SetInertiaXX(ChVector<>(20, 20, 20));
mrigidbody->SetPos(tire_center + ChVector<>(-1, 0, 0));
GetChronoDataFile("models/tractor_wheel/tractor_wheel_fine.obj"));
mesh->Transform(VNULL, Q_from_AngAxis(CH_C_PI, VECT_Y));
auto mesh_asset = chrono_types::make_shared<ChTriangleMeshShape>();
mesh_asset->SetMesh(mesh);
mrigidbody->AddVisualShape(mesh_asset);
// this is used to use the mesh in cosimulation!
auto mrigidmeshload = chrono_types::make_shared<ChLoadBodyMesh>(mrigidbody, *mesh);
mloadcontainer->Add(mrigidmeshload);
// Create the Irrlicht visualization system
auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
vis->AttachSystem(&sys);
vis->SetWindowSize(1280, 720);
vis->SetWindowTitle("demo_FEA_cosimulate_load");
vis->Initialize();
vis->AddLogo();
vis->AddSkyBox();
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));
vis->EnableShadows();
//
// THE SOFT-REAL-TIME CYCLE
//
// 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>();
sys.SetSolver(solver);
solver->SetMaxIterations(40);
solver->SetTolerance(1e-10);
solver->EnableDiagonalPreconditioner(true);
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()) {
vis->BeginScene();
vis->Render();
sys.DoStepDynamics(0.005);
// -------------------------------------------------------------------------
// Here do the cosimulation
// A <--> B
// For example, A is this main program, and B can be an external
// program, ex. a CFD or SPH 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 C function
// in this .cpp file.
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):
vert_pos.clear();
vert_vel.clear();
triangles.clear();
vert_forces.clear();
vert_indexes.clear();
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);
vis->EndScene();
}
return 0;
}
void AddTypicalLights()
Simple shortcut to set two point lights in the scene.
Definition: ChVisualSystemIrrlicht.cpp:286
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:179
COORDSYS:
Definition: ChCoordsys.h:38
void AddSkyBox(const std::string &texture_dir=GetChronoDataFile("skybox/"))
Add a sky box in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:299
void SetSolverForceTolerance(double tolerance)
Set a solver tolerance threshold at force level (default: not specified).
Definition: ChSystem.h:195
Class for exceptions for throw() catch() mechanism.
Definition: ChException.h:25
virtual void Initialize()
Initialize the visualization system.
Definition: ChVisualSystemIrrlicht.cpp:160
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:353
ChLog & GetLog()
Global function to get the current ChLog object.
Definition: ChLog.cpp:39
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:231
bool Run()
Run the Irrlicht device.
Definition: ChVisualSystemIrrlicht.cpp:217
void AddCamera(const ChVector<> &pos, ChVector<> targ=VNULL)
Add a camera in an Irrlicht 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:268
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
virtual void BeginScene(bool backBuffer=true, bool zBuffer=true, ChColor color=ChColor(0, 0, 0))
Clean the canvas at the beginning of each animation frame.
Definition: ChVisualSystemIrrlicht.cpp:501
Namespace with classes for the Irrlicht module.
Definition: ChApiIrr.h:48
Definition of a visual color.
Definition: ChColor.h:26
Definition of general purpose 3d vector variables, such as points in 3D.
Definition: ChVector.h:35
int DoStepDynamics(double step_size)
Advances the dynamical simulation for a single step, of length step_size.
Definition: ChSystem.cpp:1422
Class for a physical system in which contact is modeled using a smooth (penalty-based) method.
Definition: ChSystemSMC.h:30
virtual void EndScene()
End the scene draw at the end of each animation frame.
Definition: ChVisualSystemIrrlicht.cpp:546
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:412
void SetWindowTitle(const std::string &win_title)
Set the windoiw title (default "").
Definition: ChVisualSystemIrrlicht.cpp:124
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:144
virtual void SetSolver(std::shared_ptr< ChSolver > newsolver)
Attach a solver (derived from ChSolver) for use by this system.
Definition: ChSystem.cpp:335
virtual void Render()
Draw all 3D shapes and GUI elements at the current frame.
Definition: ChVisualSystemIrrlicht.cpp:556
Main namespace for the Chrono package.
Definition: ChBarrelShape.cpp:17
void SetTimestepperType(ChTimestepper::Type type)
Set the method for time integration (time stepper type).
Definition: ChSystem.cpp:455
void AddLogo(const std::string &logo_filename=GetChronoDataFile("logo_chronoengine_alpha.png"))
Add a logo in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:260
Namespace for FEA classes.
Definition: ChVisualShapeFEA.h:25
void SetWindowSize(unsigned int width, unsigned int height)
Set the window size (default 640x480).
Definition: ChVisualSystemIrrlicht.cpp:120
Irrlicht-based Chrono run-time visualization system.
Definition: ChVisualSystemIrrlicht.h:47