diff --git a/chronosigasimulationmodel.cpp b/chronosigasimulationmodel.cpp new file mode 100644 index 0000000..76b4bb9 --- /dev/null +++ b/chronosigasimulationmodel.cpp @@ -0,0 +1,298 @@ +#include "chronosigasimulationmodel.hpp" +#include "chrono/physics/ChLoadContainer.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace chrono; +using namespace chrono::fea; +std::shared_ptr ChronosIGASimulationModel::convertToChronosMesh_IGA( + const std::shared_ptr &pMesh, + std::vector> &edgeMeshVertsToChronosNodes) +{ + auto mesh_chronos = chrono_types::make_shared(); + 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( + 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(); + auto msection = chrono_types::make_shared( + 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( + mesh_chronos->GetNode(mesh_chronos->GetNnodes() - lastBeamNodesVector.size())); + edgeMeshVertsToChronosNodes[vi1] = std::dynamic_pointer_cast( + 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 &pJob) +//{} + +//chrono::ChSystemSMC convertToChronosSystem(const std::shared_ptr &pJob) +//{ +// chrono::ChSystemSMC my_system; +//} + +void ChronosIGASimulationModel::parseForces( + const std::shared_ptr &mesh_chronos, + const std::vector> &edgeMeshVertsToChronoNodes, + const std::unordered_map &nodalExternalForces) +{ + mesh_chronos->SetAutomaticGravity(false); + for (const std::pair &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 &pJob, + const std::vector> &edgeMeshVertsToChronoNodes, + chrono::ChSystemSMC &my_system) +{ + assert(!edgeMeshVertsToChronoNodes.empty()); + for (const std::pair> &constrainedVertex : + pJob->constrainedVertices) { + const int &constrainedVi = constrainedVertex.first; + const std::unordered_set &constrainedDoF = constrainedVertex.second; + // Create a truss + auto truss = chrono_types::make_shared(); + 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(); + 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 &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(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(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>(mesh_chronos->GetNode(1)) + // ->SetFixed(true); + // and load containers must be added to your system + // auto mvisualizemesh = chrono_types::make_shared(*(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(); + 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 rotQuat = node_chronos->GetRot(); + simulationResults.rotationalDisplacementQuaternion[vi] + = Eigen::Quaternion(rotQuat.e0(), rotQuat.e1(), rotQuat.e2(), rotQuat.e3()); + const ChVector 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 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 &pMesh) +{ + mesh_chronos = convertToChronosMesh_IGA(pMesh, edgeMeshVertsToChronoNodes); +} + +ChronosIGASimulationModel::ChronosIGASimulationModel() {} diff --git a/chronosigasimulationmodel.hpp b/chronosigasimulationmodel.hpp new file mode 100644 index 0000000..b30c6cc --- /dev/null +++ b/chronosigasimulationmodel.hpp @@ -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 mesh_chronos; + std::vector> edgeMeshVertsToChronoNodes; + + static void parseConstrainedVertices( + const std::shared_ptr &pJob, + const std::vector> &edgeMeshVertsToChronoNodes, + chrono::ChSystemSMC &my_system); + static void parseForces( + const std::shared_ptr &mesh_chronos, + const std::vector> &edgeMeshVertsToChronoNodes, + const std::unordered_map &nodalExternalForces); + +public: + ChronosIGASimulationModel(); + SimulationResults executeSimulation(const std::shared_ptr &pJob) override; + void setStructure(const std::shared_ptr &pMesh) override; + static std::shared_ptr convertToChronosMesh_IGA( + const std::shared_ptr &pMesh, + std::vector> &edgeMeshVertsToChronosNodes); + + inline const static std::string label{"Chronos_linear_IGA"}; +}; + +#endif // CHRONOSIGASIMULATIONMODEL_HPP