Perform thermal analysis (demo_FEA_thermal.cpp)

Tutorial that teaches how to use the FEA module to compute the scalar temperature field in a solid (Poisson problem) given boundary conditions.

// =============================================================================
// 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 visualization using Irrlicht
// =============================================================================
#include "chrono/physics/ChSystemSMC.h"
#include "chrono/solver/ChIterativeSolverLS.h"
#include "chrono/fea/ChElementBar.h"
#include "chrono/fea/ChElementTetraCorot_4.h"
#include "chrono/fea/ChElementTetraCorot_10.h"
#include "chrono/fea/ChElementHexaCorot_8.h"
#include "chrono/fea/ChElementHexaCorot_20.h"
#include "chrono/fea/ChContinuumThermal.h"
#include "chrono/fea/ChContinuumElectrostatics.h"
#include "chrono/fea/ChNodeFEAxyzP.h"
#include "chrono/fea/ChMesh.h"
#include "chrono/fea/ChMeshFileLoader.h"
#include "chrono/fea/ChLinkPointFrame.h"
#include "chrono/assets/ChVisualShapeFEA.h"
#include "chrono_irrlicht/ChVisualSystemIrrlicht.h"
// Remember to use the namespace 'chrono' because all classes
// of Chrono::Engine belong to this namespace and its children...
using namespace chrono;
using namespace chrono::fea;
using namespace chrono::irrlicht;
int main(int argc, char* argv[]) {
GetLog() << "Copyright (c) 2017\nChrono version: " << CHRONO_VERSION << "\n\n";
// Create a Chrono::Engine physical system
// Create a mesh, that is a container for groups
// of elements and their referenced nodes.
auto my_mesh = chrono_types::make_shared<ChMesh>();
// Create a material, that must be assigned to each element,
// and set its parameters
auto mmaterial = chrono_types::make_shared<ChContinuumThermal>();
// Add some TETAHEDRONS:
// Load a .node file and a .ele file from disk, defining a complicate tetrahedron mesh.
// This is much easier than creating all nodes and elements via C++ programming.
// You can generate these files using the TetGen tool.
try {
ChMeshFileLoader::FromTetGenFile(my_mesh, GetChronoDataFile("fea/beam.node").c_str(),
GetChronoDataFile("fea/beam.ele").c_str(), mmaterial);
} catch (const ChException& myerr) {
GetLog() << myerr.what();
return 0;
for (unsigned int inode = 0; inode < my_mesh->GetNnodes(); ++inode) {
if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(inode))) {
mnode->SetPos(mnode->GetPos() * ChVector<>(3, 1, 3));
// Set some BOUNDARY CONDITIONS on nodes:
// Impose load on the 180th node
auto mnode3 = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(180));
mnode3->SetF(20); // thermal load: heat flux [W] into node
// Impose field on two top nodes (remember the SetFixed(true); )
auto mnode1 = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(my_mesh->GetNnodes() - 1));
mnode1->SetP(0.5); // field: temperature [K]
auto mnode2 = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(my_mesh->GetNnodes() - 2));
mnode2->SetP(0.5); // field: temperature [K]
// Impose field on the base points:
for (unsigned int inode = 0; inode < my_mesh->GetNnodes(); ++inode) {
if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(inode))) {
if (mnode->GetPos().y() < 0.01) {
mnode->SetP(10); // field: temperature [K]
// Remember to add the mesh to the system!
// 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.
// Such triangle mesh can be rendered by Irrlicht or POVray or whatever
// postprocessor that can handle a colored ChVisualShapeTriangleMesh).
// Do not forget AddVisualShapeFEA() at the end!
// This will paint the colored mesh with temperature scale (NODE_P is the scalar field of the Poisson
// problem)
auto mvisualizemesh = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemesh->SetColorscaleMinMax(-1, 12);
mvisualizemesh->SetShrinkElements(false, 0.85);
// This will paint the wireframe
auto mvisualizemeshB = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
// This will paint the heat flux as line vectors
auto mvisualizemeshC = chrono_types::make_shared<ChVisualShapeFEA>(my_mesh);
mvisualizemeshC->SetDefaultSymbolsColor(ChColor(0.1f, 0.2f, 0.2f));
// Create the Irrlicht visualization system
auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
vis->SetWindowSize(800, 600);
vis->SetWindowTitle("FEM thermal");
vis->AddLight(ChVector<>(+20, 20, +20), 90, ChColor(0.5, 0.5, 0.5));
vis->AddLight(ChVector<>(-20, 20, -20), 90, ChColor(0.7f, 0.8f, 0.8f));
vis->AddCamera(ChVector<>(0, 0.7, -1), ChVector<>(0, 0.4, 0));
// Use MINRES solver to handle stiffness matrices.
auto solver = chrono_types::make_shared<ChSolverMINRES>();
solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED
sys.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED); // fast, less precise
// Note: if you are interested only in a single LINEAR STATIC solution
// (not a transient thermal solution, but rather the steady-state solution),
// at this point you can uncomment the following line:
// sys.DoStaticLinear();
// Also, in the following while() loop, remove application.DoStep();
// so you can spin the 3D view and look at the solution.
while (vis->Run()) {
if (sys.GetChTime() > 5)
// Print some node temperatures..
for (unsigned int inode = 0; inode < my_mesh->GetNnodes(); ++inode) {
if (auto mnode = std::dynamic_pointer_cast<ChNodeFEAxyzP>(my_mesh->GetNode(inode))) {
if (mnode->GetPos().x() < 0.01) {
GetLog() << "Node at y=" << mnode->GetPos().y() << " has T=" << mnode->GetP() << "\n";
return 0;
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
void AddSkyBox(const std::string &texture_dir=GetChronoDataFile("skybox/"))
Add a sky box in a 3D scene.
Definition: ChVisualSystemIrrlicht.cpp:350
Class for exceptions for throw() catch() mechanism.
Definition: ChException.h:24
virtual void Initialize()
Initialize the visualization system.
Definition: ChVisualSystemIrrlicht.cpp:181
irr::scene::ILightSceneNode * AddLight(const ChVector<> &pos, double radius, ChColor color=ChColor(0.7f, 0.7f, 0.7f))
Add a point light to the scene.
Definition: ChVisualSystemIrrlicht.cpp:395
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
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
static void FromTetGenFile(std::shared_ptr< ChMesh > mesh, const char *filename_node, const char *filename_ele, std::shared_ptr< ChContinuumMaterial > my_material, ChVector<> pos_transform=VNULL, ChMatrix33<> rot_transform=ChMatrix33<>(1))
Load tetrahedrons from .node and .ele files as saved by TetGen.
Definition: ChMeshFileLoader.cpp:39
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
void SetWindowTitle(const std::string &win_title)
Set the windoiw title (default "").
Definition: ChVisualSystemIrrlicht.cpp:138
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
double GetChTime() const
Get the simulation time of this system.
Definition: ChSystem.h:212
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
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