Added chronos IGA simulation model
This commit is contained in:
iasonmanolas
parent
d8a6fadab7
commit
3092809d02
|
|
@ -0,0 +1,298 @@
|
|||
#include "chronosigasimulationmodel.hpp"
|
||||
#include "chrono/physics/ChLoadContainer.h"
|
||||
#include <chrono/fea/ChBeamSectionEuler.h>
|
||||
#include <chrono/fea/ChBuilderBeam.h>
|
||||
#include <chrono/fea/ChLoadsBeam.h>
|
||||
#include <chrono/fea/ChMesh.h>
|
||||
#include <chrono/fea/ChNodeFEAxyzrot.h>
|
||||
#include <chrono/physics/ChBody.h>
|
||||
#include <chrono/physics/ChSystemSMC.h>
|
||||
#include <chrono/solver/ChIterativeSolverLS.h>
|
||||
|
||||
#include <chrono/assets/ChVisualization.h>
|
||||
#include <chrono/fea/ChElementBeamEuler.h>
|
||||
#include <chrono/fea/ChVisualizationFEAmesh.h>
|
||||
|
||||
using namespace chrono;
|
||||
using namespace chrono::fea;
|
||||
std::shared_ptr<ChMesh> ChronosIGASimulationModel::convertToChronosMesh_IGA(
|
||||
const std::shared_ptr<SimulationMesh> &pMesh,
|
||||
std::vector<std::shared_ptr<ChNodeFEAxyzrot>> &edgeMeshVertsToChronosNodes)
|
||||
{
|
||||
auto mesh_chronos = chrono_types::make_shared<ChMesh>();
|
||||
edgeMeshVertsToChronosNodes.clear();
|
||||
edgeMeshVertsToChronosNodes.resize(pMesh->VN(), nullptr);
|
||||
|
||||
//add nodes
|
||||
for (int vi = 0; vi < pMesh->VN(); vi++) {
|
||||
const auto &vertex = pMesh->vert[vi];
|
||||
ChVector<> vertexPos(vertex.cP()[0], vertex.cP()[1], vertex.cP()[2]);
|
||||
edgeMeshVertsToChronosNodes[vi] = chrono_types::make_shared<ChNodeFEAxyzrot>(
|
||||
ChFrame<>(vertexPos));
|
||||
mesh_chronos->AddNode(edgeMeshVertsToChronosNodes[vi]);
|
||||
}
|
||||
//add elements
|
||||
ChBuilderBeamIGA builder;
|
||||
for (int ei = 0; ei < pMesh->EN(); ei++) {
|
||||
const SimulationMesh::EdgeType &edge = pMesh->edge[ei];
|
||||
//define end-points
|
||||
const auto vi0 = pMesh->getIndex(edge.cV(0));
|
||||
const auto vi1 = pMesh->getIndex(edge.cV(1));
|
||||
//define cross section
|
||||
const Element &element = pMesh->elements[ei];
|
||||
const double beam_wz = element.dimensions.getDim1();
|
||||
const double beam_wy = element.dimensions.getDim2();
|
||||
const double E = element.material.youngsModulus;
|
||||
// const double poisson = element.material.poissonsRatio;
|
||||
const double density = 1e0;
|
||||
|
||||
// auto msection = chrono_types::make_shared<ChBeamSectionEulerAdvanced>();
|
||||
auto msection = chrono_types::make_shared<ChBeamSectionCosseratEasyRectangular>(
|
||||
beam_wy, beam_wz, E, element.material.G, density);
|
||||
// msection->SetDensity(density);
|
||||
// msection->SetYoungModulus(E);
|
||||
// msection->SetGwithPoissonRatio(poisson);
|
||||
// // msection->SetBeamRaleyghDamping(0.0);
|
||||
// msection->SetAsRectangularSection(beam_wy, beam_wz);
|
||||
builder.BuildBeam(
|
||||
mesh_chronos, // the mesh where to put the created nodes and elements
|
||||
msection, // the ChBeamSectionEuler to use for the ChElementBeamEuler elements
|
||||
4, // the number of ChElementBeamEuler to create
|
||||
edgeMeshVertsToChronosNodes[vi0]->GetPos(), // the 'A' point in space (beginning of beam)
|
||||
edgeMeshVertsToChronosNodes[vi1]->GetPos(), // the 'B' point in space (end of beam)
|
||||
ChVector<>(0, 0, 1)
|
||||
// ChVector<>(0, cos(rot_rad), sin(rot_rad))
|
||||
); // the 'Y' up direction of the section for the beam
|
||||
const auto lastBeamNodesVector = builder.GetLastBeamNodes();
|
||||
// assert(lastBeamNodesVector.size() == 4);
|
||||
edgeMeshVertsToChronosNodes[vi0] = std::dynamic_pointer_cast<ChNodeFEAxyzrot>(
|
||||
mesh_chronos->GetNode(mesh_chronos->GetNnodes() - lastBeamNodesVector.size()));
|
||||
edgeMeshVertsToChronosNodes[vi1] = std::dynamic_pointer_cast<ChNodeFEAxyzrot>(
|
||||
mesh_chronos->GetNode(mesh_chronos->GetNnodes() - 1));
|
||||
// std::cout << edgeMeshVertsToChronosNodes[vi1]->GetCoord().pos[0] << " "
|
||||
// << edgeMeshVertsToChronosNodes[vi1]->GetCoord().pos[1] << " "
|
||||
// << edgeMeshVertsToChronosNodes[vi1]->GetCoord().pos[2] << std::endl;
|
||||
// std::cout << lastBeamNodesVector.back()->GetCoord().pos[0] << " "
|
||||
// << lastBeamNodesVector.back()->GetCoord().pos[1] << " "
|
||||
// << lastBeamNodesVector.back()->GetCoord().pos[2] << std::endl;
|
||||
// int i = 0;
|
||||
// i++;
|
||||
}
|
||||
|
||||
return mesh_chronos;
|
||||
}
|
||||
|
||||
//SimulationResults ChronosEulerSimulationModel::executeSimulation(
|
||||
// const std::shared_ptr<SimulationJob> &pJob)
|
||||
//{}
|
||||
|
||||
//chrono::ChSystemSMC convertToChronosSystem(const std::shared_ptr<SimulationJob> &pJob)
|
||||
//{
|
||||
// chrono::ChSystemSMC my_system;
|
||||
//}
|
||||
|
||||
void ChronosIGASimulationModel::parseForces(
|
||||
const std::shared_ptr<chrono::fea::ChMesh> &mesh_chronos,
|
||||
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronoNodes,
|
||||
const std::unordered_map<VertexIndex, Vector6d> &nodalExternalForces)
|
||||
{
|
||||
mesh_chronos->SetAutomaticGravity(false);
|
||||
for (const std::pair<VertexIndex, Vector6d> &externalForce : nodalExternalForces) {
|
||||
const int &forceVi = externalForce.first;
|
||||
const Vector6d &force = externalForce.second;
|
||||
edgeMeshVertsToChronoNodes[forceVi]->SetForce(ChVector<>(force[0], force[1], force[2]));
|
||||
edgeMeshVertsToChronoNodes[forceVi]->SetTorque(ChVector<>(force[3], force[4], force[5]));
|
||||
std::cout << "Force on vertex:" << std::endl;
|
||||
std::cout << edgeMeshVertsToChronoNodes[forceVi]->GetCoord().pos[0] << " "
|
||||
<< edgeMeshVertsToChronoNodes[forceVi]->GetCoord().pos[1] << " "
|
||||
<< edgeMeshVertsToChronoNodes[forceVi]->GetCoord().pos[2] << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
void ChronosIGASimulationModel::parseConstrainedVertices(
|
||||
const std::shared_ptr<const SimulationJob> &pJob,
|
||||
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronoNodes,
|
||||
chrono::ChSystemSMC &my_system)
|
||||
{
|
||||
assert(!edgeMeshVertsToChronoNodes.empty());
|
||||
for (const std::pair<VertexIndex, std::unordered_set<int>> &constrainedVertex :
|
||||
pJob->constrainedVertices) {
|
||||
const int &constrainedVi = constrainedVertex.first;
|
||||
const std::unordered_set<int> &constrainedDoF = constrainedVertex.second;
|
||||
// Create a truss
|
||||
auto truss = chrono_types::make_shared<ChBody>();
|
||||
truss->SetBodyFixed(true);
|
||||
my_system.Add(truss);
|
||||
const auto &constrainedChronoNode = edgeMeshVertsToChronoNodes[constrainedVi];
|
||||
std::cout << "Constrained vertex:" << std::endl;
|
||||
std::cout << edgeMeshVertsToChronoNodes[constrainedVi]->GetCoord().pos[0] << " "
|
||||
<< edgeMeshVertsToChronoNodes[constrainedVi]->GetCoord().pos[1] << " "
|
||||
<< edgeMeshVertsToChronoNodes[constrainedVi]->GetCoord().pos[2] << std::endl;
|
||||
// Create a constraint at the end of the beam
|
||||
auto constr_a = chrono_types::make_shared<ChLinkMateGeneric>();
|
||||
constr_a->SetConstrainedCoords(constrainedDoF.contains(0),
|
||||
constrainedDoF.contains(1),
|
||||
constrainedDoF.contains(2),
|
||||
constrainedDoF.contains(3),
|
||||
constrainedDoF.contains(4),
|
||||
constrainedDoF.contains(5));
|
||||
constr_a->Initialize(constrainedChronoNode,
|
||||
truss,
|
||||
false,
|
||||
constrainedChronoNode->Frame(),
|
||||
constrainedChronoNode->Frame());
|
||||
// const auto frameNode = constrainedChronoNode->Frame();
|
||||
my_system.Add(constr_a);
|
||||
|
||||
// edgeMeshVertsToChronoNodes[constrainedVi]->SetFixed(true);
|
||||
// if (vertexIsFullyConstrained) {
|
||||
// } else {
|
||||
// std::cerr << "Currently only rigid vertices are handled." << std::endl;
|
||||
// // SimulationResults simulationResults;
|
||||
// // simulationResults.converged = false;
|
||||
// // assert(false);
|
||||
// // return simulationResults;
|
||||
// }
|
||||
}
|
||||
}
|
||||
|
||||
SimulationResults ChronosIGASimulationModel::executeSimulation(
|
||||
const std::shared_ptr<SimulationJob> &pJob)
|
||||
{
|
||||
// assert(pJob->pMesh->VN() != 0);
|
||||
// const bool structureInitialized = mesh_chronos != nullptr;
|
||||
// const bool wasInitializedWithDifferentStructure = structureInitialized
|
||||
// && (pJob->pMesh->EN()
|
||||
// != mesh_chronos->GetNelements()
|
||||
// || pJob->pMesh->VN()
|
||||
// != mesh_chronos->GetNnodes());
|
||||
chrono::ChSystemSMC my_system;
|
||||
// if (!structureInitialized || wasInitializedWithDifferentStructure) {
|
||||
setStructure(pJob->pMesh);
|
||||
my_system.Add(mesh_chronos);
|
||||
// }
|
||||
// chrono::irrlicht::ChIrrApp application(&my_system,
|
||||
// L"Irrlicht FEM visualization",
|
||||
// irr::core::dimension2d<irr::u32>(800, 600),
|
||||
// chrono::irrlicht::VerticalDir::Z,
|
||||
// false,
|
||||
// true);
|
||||
// const std::string chronoDataFolderPath = "/home/iason/Coding/build/external "
|
||||
// "dependencies/CHRONO-src/data/";
|
||||
// application.AddTypicalLogo(chronoDataFolderPath + "logo_chronoengine_alpha.png");
|
||||
// application.AddTypicalSky(chronoDataFolderPath + "skybox/");
|
||||
// application.AddTypicalLights();
|
||||
// application.AddTypicalCamera(irr::core::vector3df(0, (irr::f32) 0.6, -1));
|
||||
// my_system.SetTimestepperType(chrono::ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED);
|
||||
const auto node_0 = std::dynamic_pointer_cast<ChNodeFEAxyzrot>(mesh_chronos->GetNode(0));
|
||||
//parse forces
|
||||
// std::cout << node_0->GetCoord().pos[0] << " " << node_0->GetCoord().pos[1] << " "
|
||||
// << node_0->GetCoord().pos[2] << std::endl;
|
||||
parseForces(mesh_chronos, edgeMeshVertsToChronoNodes, pJob->nodalExternalForces);
|
||||
// std::cout << node_0->GetCoord().pos[0] << " " << node_0->GetCoord().pos[1] << " "
|
||||
// << node_0->GetCoord().pos[2] << std::endl;
|
||||
|
||||
//parse constrained vertices
|
||||
// std::cout << node_0->GetCoord().pos[0] << " " << node_0->GetCoord().pos[1] << " "
|
||||
// << node_0->GetCoord().pos[2] << std::endl;
|
||||
parseConstrainedVertices(pJob, edgeMeshVertsToChronoNodes, my_system);
|
||||
// std::cout << node_0->GetCoord().pos[0] << " " << node_0->GetCoord().pos[1] << " "
|
||||
// << node_0->GetCoord().pos[2] << std::endl;
|
||||
// std::dynamic_pointer_cast<std::shared_ptr<ChNodeFEAxyzrot>>(mesh_chronos->GetNode(1))
|
||||
// ->SetFixed(true);
|
||||
// and load containers must be added to your system
|
||||
// auto mvisualizemesh = chrono_types::make_shared<ChVisualizationFEAmesh>(*(mesh_chronos.get()));
|
||||
// mvisualizemesh->SetFEMdataType(ChVisualizationFEAmesh::E_PLOT_NODE_DISP_NORM);
|
||||
// mvisualizemesh->SetColorscaleMinMax(0.0, 5.50);
|
||||
// mvisualizemesh->SetShrinkElements(false, 0.85);
|
||||
// mvisualizemesh->SetSmoothFaces(false);
|
||||
// mesh_chronos->AddAsset(mvisualizemesh);
|
||||
|
||||
// application.AssetBindAll();
|
||||
// application.AssetUpdateAll();
|
||||
|
||||
auto solver = chrono_types::make_shared<ChSolverMINRES>();
|
||||
my_system.SetSolver(solver);
|
||||
solver->SetMaxIterations(1e5);
|
||||
// solver->SetTolerance(1e-12);
|
||||
solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED
|
||||
solver->EnableDiagonalPreconditioner(true);
|
||||
// my_system.SetSolverForceTolerance(1e-9);
|
||||
solver->SetVerbose(true);
|
||||
|
||||
std::cout << node_0->GetCoord().pos[0] << " " << node_0->GetCoord().pos[1] << " "
|
||||
<< node_0->GetCoord().pos[2] << std::endl;
|
||||
SimulationResults simulationResults;
|
||||
simulationResults.converged = my_system.DoStaticNonlinear();
|
||||
// simulationResults.converged = my_system.DoStaticLinear();
|
||||
std::cout << node_0->GetCoord().pos[0] << " " << node_0->GetCoord().pos[1] << " "
|
||||
<< node_0->GetCoord().pos[2] << std::endl;
|
||||
if (!simulationResults.converged) {
|
||||
std::cerr << "Simulation failed" << std::endl;
|
||||
assert(false);
|
||||
return simulationResults;
|
||||
}
|
||||
|
||||
// my_system.SetTimestepperType(ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED);
|
||||
// application.SetTimestep(0.001);
|
||||
|
||||
// while (application.GetDevice()->run()) {
|
||||
// application.BeginScene();
|
||||
// application.DrawAll();
|
||||
// application.EndScene();
|
||||
// }
|
||||
simulationResults.pJob = pJob;
|
||||
simulationResults.displacements.resize(pJob->pMesh->VN());
|
||||
simulationResults.rotationalDisplacementQuaternion.resize(pJob->pMesh->VN());
|
||||
for (size_t vi = 0; vi < pJob->pMesh->VN(); vi++) {
|
||||
const auto node_chronos = edgeMeshVertsToChronoNodes[vi];
|
||||
std::cout << node_chronos->GetCoord().pos[0] << " " << node_chronos->GetCoord().pos[1]
|
||||
<< " " << node_chronos->GetCoord().pos[2] << std::endl;
|
||||
int i = 0;
|
||||
i++;
|
||||
assert(node_chronos != nullptr);
|
||||
const auto posDisplacement = node_chronos->Frame().GetPos()
|
||||
- node_chronos->GetX0().GetPos();
|
||||
// std::cout << "Node " << vi << " coordinate x= " << node_chronos->Frame().GetPos().x()
|
||||
// << " y=" << node_chronos->Frame().GetPos().y()
|
||||
// << " z=" << node_chronos->Frame().GetPos().z() << "\n";
|
||||
//Translations
|
||||
simulationResults.displacements[vi][0] = posDisplacement[0];
|
||||
simulationResults.displacements[vi][1] = posDisplacement[1];
|
||||
simulationResults.displacements[vi][2] = posDisplacement[2];
|
||||
//Rotations
|
||||
chrono::ChQuaternion<double> rotQuat = node_chronos->GetRot();
|
||||
simulationResults.rotationalDisplacementQuaternion[vi]
|
||||
= Eigen::Quaternion<double>(rotQuat.e0(), rotQuat.e1(), rotQuat.e2(), rotQuat.e3());
|
||||
const ChVector<double> eulerAngles = rotQuat.Q_to_Euler123();
|
||||
// std::cout << "Euler angles:" << eulerAngles << std::endl;
|
||||
simulationResults.displacements[vi][3] = eulerAngles[0];
|
||||
simulationResults.displacements[vi][4] = eulerAngles[1];
|
||||
simulationResults.displacements[vi][5] = eulerAngles[2];
|
||||
}
|
||||
|
||||
simulationResults.simulationModelUsed = label;
|
||||
return simulationResults;
|
||||
|
||||
// VCGEdgeMesh deformedMesh;
|
||||
// deformedMesh.copy(*(pJob->pMesh));
|
||||
// for (size_t vi = 0; vi < pJob->pMesh->VN(); vi++) {
|
||||
// const std::shared_ptr<ChNodeFEAxyzrot> node_chronos = edgeMeshVertsToChronosNodes[vi];
|
||||
// deformedMesh.vert[vi].P() = CoordType(node_chronos->GetPos()[0],
|
||||
// node_chronos->GetPos()[1],
|
||||
// node_chronos->GetPos()[2]);
|
||||
// }
|
||||
|
||||
// deformedMesh.updateEigenEdgeAndVertices();
|
||||
// deformedMesh.setLabel("deformed");
|
||||
// deformedMesh.registerForDrawing();
|
||||
// polyscope::show();
|
||||
// return simulationResults;
|
||||
}
|
||||
|
||||
void ChronosIGASimulationModel::setStructure(const std::shared_ptr<SimulationMesh> &pMesh)
|
||||
{
|
||||
mesh_chronos = convertToChronosMesh_IGA(pMesh, edgeMeshVertsToChronoNodes);
|
||||
}
|
||||
|
||||
ChronosIGASimulationModel::ChronosIGASimulationModel() {}
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
#ifndef CHRONOSIGASIMULATIONMODEL_HPP
|
||||
#define CHRONOSIGASIMULATIONMODEL_HPP
|
||||
|
||||
#include "simulationmodel.hpp"
|
||||
|
||||
namespace chrono {
|
||||
class ChSystemSMC;
|
||||
namespace fea {
|
||||
class ChMesh;
|
||||
class ChNodeFEAxyzrot;
|
||||
} // namespace fea
|
||||
} // namespace chrono
|
||||
|
||||
class ChronosIGASimulationModel : public SimulationModel
|
||||
{
|
||||
std::shared_ptr<chrono::fea::ChMesh> mesh_chronos;
|
||||
std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> edgeMeshVertsToChronoNodes;
|
||||
|
||||
static void parseConstrainedVertices(
|
||||
const std::shared_ptr<const SimulationJob> &pJob,
|
||||
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronoNodes,
|
||||
chrono::ChSystemSMC &my_system);
|
||||
static void parseForces(
|
||||
const std::shared_ptr<chrono::fea::ChMesh> &mesh_chronos,
|
||||
const std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronoNodes,
|
||||
const std::unordered_map<VertexIndex, Vector6d> &nodalExternalForces);
|
||||
|
||||
public:
|
||||
ChronosIGASimulationModel();
|
||||
SimulationResults executeSimulation(const std::shared_ptr<SimulationJob> &pJob) override;
|
||||
void setStructure(const std::shared_ptr<SimulationMesh> &pMesh) override;
|
||||
static std::shared_ptr<chrono::fea::ChMesh> convertToChronosMesh_IGA(
|
||||
const std::shared_ptr<SimulationMesh> &pMesh,
|
||||
std::vector<std::shared_ptr<chrono::fea::ChNodeFEAxyzrot>> &edgeMeshVertsToChronosNodes);
|
||||
|
||||
inline const static std::string label{"Chronos_linear_IGA"};
|
||||
};
|
||||
|
||||
#endif // CHRONOSIGASIMULATIONMODEL_HPP
|
||||
Loading…
Reference in New Issue