Use bricks for solid meshes (demo_FEA_brick.cpp)

Tutorial that teaches how to use the FEA module to perform FEA dynamics for bricks of ANCF type.

// =============================================================================
// 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 using the brick element
//
// =============================================================================
#include "chrono/physics/ChSystemSMC.h"
#include "chrono/fea/ChElementBrick.h"
#include "chrono/fea/ChElementBar.h"
#include "chrono/fea/ChLinkPointFrame.h"
#include "chrono/fea/ChLinkDirFrame.h"
#include "chrono/fea/ChVisualizationFEAmesh.h"
#include "chrono/solver/ChIterativeSolverLS.h"
#include "chrono_irrlicht/ChIrrApp.h"
using namespace chrono;
using namespace chrono::fea;
using namespace chrono::irrlicht;
using namespace irr;
int main(int argc, char* argv[]) {
GetLog() << "Copyright (c) 2017 projectchrono.org\nChrono version: " << CHRONO_VERSION << "\n\n";
ChSystemSMC my_system;
// Create the Irrlicht visualization (open the Irrlicht device,
// bind a simple user interface, etc. etc.)
ChIrrApp application(&my_system, L"Brick Elements", core::dimension2d<u32>(800, 600));
// Easy shortcuts to add camera, lights, logo and sky in Irrlicht scene:
application.AddTypicalLogo();
application.AddTypicalSky();
application.AddTypicalLights();
application.AddTypicalCamera(core::vector3df(1.2f, 0.6f, 0.3f), // camera location
core::vector3df(0.2f, -0.2f, 0.f)); // "look at" location
GetLog() << "-----------------------------------------------------------\n";
GetLog() << "-----------------------------------------------------------\n";
GetLog() << " Brick Elements demo with implicit integration \n";
GetLog() << "-----------------------------------------------------------\n";
// The physical system: it contains all physical objects.
// Create a mesh, that is a container for groups
// of elements and their referenced nodes.
auto my_mesh = chrono_types::make_shared<ChMesh>();
// Geometry of the plate
double plate_lenght_x = 1;
double plate_lenght_y = 1;
double plate_lenght_z = 0.05; // small thickness
// Specification of the mesh
int numDiv_x = 4;
int numDiv_y = 4;
int numDiv_z = 1;
int N_x = numDiv_x + 1;
// Number of elements in the z direction is considered as 1
int TotalNumElements = numDiv_x * numDiv_y;
int TotalNumNodes = (numDiv_x + 1) * (numDiv_y + 1) * (numDiv_z + 1);
// For uniform mesh
double dx = plate_lenght_x / numDiv_x;
double dy = plate_lenght_y / numDiv_y;
double dz = plate_lenght_z / numDiv_z;
int MaxMNUM = 1;
int MaxLayNum = 1;
ChMatrixDynamic<double> COORDFlex(TotalNumNodes, 3);
ChMatrixDynamic<double> VELCYFlex(TotalNumNodes, 3);
ChMatrixDynamic<int> NumNodes(TotalNumElements, 8);
ChMatrixDynamic<int> LayNum(TotalNumElements, 1);
ChMatrixDynamic<int> NDR(TotalNumNodes, 3);
ChMatrixDynamic<double> ElemLengthXY(TotalNumElements, 3);
for (int i = 0; i < MaxMNUM; i++) {
MPROP(i, 0) = 500; // Density [kg/m3]
MPROP(i, 1) = 2.1E+05; // H(m)
MPROP(i, 2) = 0.3; // nu
}
auto mmaterial = chrono_types::make_shared<ChContinuumElastic>();
mmaterial->Set_RayleighDampingK(0.0);
mmaterial->Set_RayleighDampingM(0.0);
mmaterial->Set_density(MPROP(0, 0));
mmaterial->Set_E(MPROP(0, 1));
mmaterial->Set_G(MPROP(0, 1) / (2 + 2 * MPROP(0, 2)));
mmaterial->Set_v(MPROP(0, 2));
for (int i = 0; i < TotalNumElements; i++) {
// All the elements belong to the same layer, e.g layer number 1.
LayNum(i, 0) = 1;
// Node number of the 4 nodes which creates element i.
// The nodes are distributed this way. First in the x direction for constant
// y when max x is reached go to the
// next level for y by doing the same distribution but for y+1 and keep
// doing until y limit is reached. Node
// number start from 1.
NumNodes(i, 0) = (i / (numDiv_x)) * (N_x) + i % numDiv_x;
NumNodes(i, 1) = (i / (numDiv_x)) * (N_x) + i % numDiv_x + 1;
NumNodes(i, 2) = (i / (numDiv_x)) * (N_x) + i % numDiv_x + 1 + N_x;
NumNodes(i, 3) = (i / (numDiv_x)) * (N_x) + i % numDiv_x + N_x;
NumNodes(i, 4) = (numDiv_x + 1) * (numDiv_y + 1) + NumNodes(i, 0);
NumNodes(i, 5) = (numDiv_x + 1) * (numDiv_y + 1) + NumNodes(i, 1);
NumNodes(i, 6) = (numDiv_x + 1) * (numDiv_y + 1) + NumNodes(i, 2);
NumNodes(i, 7) = (numDiv_x + 1) * (numDiv_y + 1) + NumNodes(i, 3);
// Let's keep the element length a fixed number in both direction. (uniform
// distribution of nodes in both direction)
ElemLengthXY(i, 0) = dx;
ElemLengthXY(i, 1) = dy;
ElemLengthXY(i, 2) = dz;
if (MaxLayNum < LayNum(i, 0)) {
MaxLayNum = LayNum(i, 0);
}
}
for (int i = 0; i < TotalNumNodes; i++) {
// If the node is the first node from the left side fix the x,y,z degree of
// freedom. 1 for constrained 0 for ...
//-The NDR array is used to define the degree of freedoms that are
// constrained in the 6 variable explained above.
NDR(i, 0) = (i % (numDiv_x + 1) == 0) ? 1 : 0;
NDR(i, 1) = (i % (numDiv_x + 1) == 0) ? 1 : 0;
NDR(i, 2) = (i % (numDiv_x + 1) == 0) ? 1 : 0;
//-COORDFlex are the initial coordinates for each node,
// the first three are the position
COORDFlex(i, 0) = (i % (numDiv_x + 1)) * dx;
COORDFlex(i, 1) = (i / (numDiv_x + 1)) % (numDiv_y + 1) * dy;
COORDFlex(i, 2) = (i) / ((numDiv_x + 1) * (numDiv_y + 1)) * dz;
//-VELCYFlex is essentially the same as COORDFlex, but for the initial
// velocity instead of position.
// let's assume zero initial velocity for nodes
VELCYFlex(i, 0) = 0;
VELCYFlex(i, 1) = 0;
VELCYFlex(i, 2) = 0;
}
// Adding the nodes to the mesh
int i = 0;
while (i < TotalNumNodes) {
auto node = chrono_types::make_shared<ChNodeFEAxyz>(ChVector<>(COORDFlex(i, 0), COORDFlex(i, 1), COORDFlex(i, 2)));
node->SetMass(0.0);
my_mesh->AddNode(node);
if (NDR(i, 0) == 1 && NDR(i, 1) == 1 && NDR(i, 2) == 1) {
node->SetFixed(true);
}
i++;
}
auto nodetip = std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(TotalNumNodes - 1));
int elemcount = 0;
while (elemcount < TotalNumElements) {
auto element = chrono_types::make_shared<ChElementBrick>();
ChVector<double> InertFlexVec(ElemLengthXY(elemcount, 0), ElemLengthXY(elemcount, 1),
ElemLengthXY(elemcount, 2)); // read element length, used in ChElementBrick
element->SetInertFlexVec(InertFlexVec);
element->SetNodes(std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 0))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 1))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 2))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 3))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 4))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 5))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 6))),
std::dynamic_pointer_cast<ChNodeFEAxyz>(my_mesh->GetNode(NumNodes(elemcount, 7))));
element->SetMaterial(mmaterial);
element->SetElemNum(elemcount); // for EAS
element->SetGravityOn(true); // turn gravity on/off from within the element
element->SetMooneyRivlin(false); // turn on/off Mooney Rivlin (Linear Isotropic by default)
ChVectorN<double, 9> stock_alpha_EAS; //
stock_alpha_EAS.setZero();
element->SetStockAlpha(stock_alpha_EAS(0), stock_alpha_EAS(1), stock_alpha_EAS(2),
stock_alpha_EAS(3), stock_alpha_EAS(4), stock_alpha_EAS(5),
stock_alpha_EAS(6), stock_alpha_EAS(7), stock_alpha_EAS(8));
my_mesh->AddElement(element);
elemcount++;
}
// Deactivate automatic gravity in mesh
my_mesh->SetAutomaticGravity(false);
my_system.Set_G_acc(ChVector<>(0, 0, -9.81));
// Remember to add the mesh to the system!
my_system.Add(my_mesh);
// Options for visualization in irrlicht
auto mvisualizemesh = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
mvisualizemesh->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_NODE_P);
mvisualizemesh->SetShrinkElements(true, 0.85);
mvisualizemesh->SetSmoothFaces(false);
my_mesh->AddAsset(mvisualizemesh);
auto mvisualizemeshref = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
mvisualizemeshref->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_SURFACE);
mvisualizemeshref->SetWireframe(true);
mvisualizemeshref->SetDrawInUndeformedReference(true);
my_mesh->AddAsset(mvisualizemeshref);
auto mvisualizemeshC = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
mvisualizemeshC->SetFEMglyphType(ChVisualizationFEAmesh::E_GLYPH_NODE_DOT_POS);
mvisualizemeshC->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_SURFACE);
mvisualizemeshC->SetSymbolsThickness(0.015);
my_mesh->AddAsset(mvisualizemeshC);
auto mvisualizemeshD = chrono_types::make_shared<ChVisualizationFEAmesh>(*(my_mesh.get()));
mvisualizemeshD->SetFEMglyphType(ChVisualizationFEAmesh::E_GLYPH_NONE);
mvisualizemeshD->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_SURFACE);
mvisualizemeshD->SetSymbolsScale(1);
mvisualizemeshD->SetColorscaleMinMax(-0.5, 5);
mvisualizemeshD->SetZbufferHide(false);
my_mesh->AddAsset(mvisualizemeshD);
application.AssetBindAll();
application.AssetUpdateAll();
// Perform a dynamic time integration:
auto solver = chrono_types::make_shared<ChSolverMINRES>();
my_system.SetSolver(solver);
solver->SetMaxIterations(1000);
solver->SetTolerance(1e-10);
solver->EnableDiagonalPreconditioner(true);
solver->SetVerbose(false);
my_system.SetTimestepperType(ChTimestepper::Type::HHT);
auto mystepper = std::dynamic_pointer_cast<ChTimestepperHHT>(my_system.GetTimestepper());
mystepper->SetAlpha(-0.2);
mystepper->SetMaxiters(100);
mystepper->SetAbsTolerances(1e-5);
mystepper->SetMode(ChTimestepperHHT::POSITION);
mystepper->SetScaling(true);
application.SetTimestep(0.004);
while (application.GetDevice()->run()) {
application.BeginScene();
application.DrawAll();
application.DoStep();
application.EndScene();
}
return 0;
}
Eigen::Matrix< T, N, 1 > ChVectorN
Column vector with fixed size (known at compile time).
Definition: ChMatrix.h:98
void Add(std::shared_ptr< ChPhysicsItem > item)
Attach an arbitrary ChPhysicsItem (e.g.
Definition: ChSystem.cpp:146
std::shared_ptr< ChTimestepper > GetTimestepper() const
Get the timestepper currently used for time integration.
Definition: ChSystem.h:136
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > ChMatrixDynamic
Dense matrix with dynamic size (i.e., with size a run-time variable, unknown at compile time).
Definition: ChMatrix.h:73
ChLog & GetLog()
Global function to get the current ChLog object.
Definition: ChLog.cpp:39
void SetTimestep(double val)
Set/Get the time step for time integration.
Definition: ChIrrAppInterface.cpp:552
Class to add some GUI to Irrlicht+ChronoEngine applications.
Definition: ChIrrApp.h:29
virtual void EndScene()
Call this to end the scene draw at the end of each animation frame.
Definition: ChIrrAppInterface.cpp:627
Definition of general purpose 3d vector variables, such as points in 3D.
Definition: ChVector.h:35
virtual void DoStep()
Call this function inside a loop such as.
Definition: ChIrrAppInterface.cpp:637
Class for a physical system in which contact is modeled using a smooth (penalty-based) method.
Definition: ChSystemSMC.h:31
virtual void SetSolver(std::shared_ptr< ChSolver > newsolver)
Attach a solver (derived from ChSolver) for use by this system.
Definition: ChSystem.cpp:270
virtual void BeginScene(bool backBuffer=true, bool zBuffer=true, irr::video::SColor color=irr::video::SColor(255, 0, 0, 0))
Call this to clean the canvas at the beginning of each animation frame.
Definition: ChIrrAppInterface.cpp:610
Main namespace for the Chrono package.
Definition: ChAsset.cpp:18
void SetTimestepperType(ChTimestepper::Type type)
Set the method for time integration (time stepper type).
Definition: ChSystem.cpp:348
void Set_G_acc(const ChVector<> &m_acc)
Set the G (gravity) acceleration vector, affecting all the bodies in the system.
Definition: ChSystem.h:224
virtual void DrawAll()
Call this function inside a loop such as.
Definition: ChIrrAppInterface.cpp:750
Namespace for FEA classes.
Definition: ChChrono.h:52
Eigen::Matrix< T, M, N, Eigen::RowMajor > ChMatrixNM
Dense matrix with fixed size (known at compile time).
Definition: ChMatrix.h:78