diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d01355..b20f154 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -156,6 +156,13 @@ PUBLIC ${ENSMALLEN_SOURCE_DIR}/include) add_compile_definitions(USE_ENSMALLEN) #endif() +#set(Chrono_DIR "/home/iason/Coding/build/external dependencies/CHRONO-src/cmake") +#find_package(Chrono +# CONFIG REQUIRED) +#target_include_directories(${PROJECT_NAME} PRIVATE ${CHRONO_INCLUDE_DIRS} ) +target_link_libraries(${PROJECT_NAME} PRIVATE "/home/iason/Coding/build/external dependencies/CHRONO-src/build/RelWithDebInfo/lib/libChronoEngine.so") +target_include_directories(${PROJECT_NAME} PRIVATE "/home/iason/Coding/build/external dependencies/CHRONO-src/src" PRIVATE "/home/iason/Coding/build/external dependencies/CHRONO-src/build/RelWithDebInfo/" ) + #find_package(OpenMP REQUIRED) #target_link_libraries(${PROJECT_NAME} OpenMP::OpenMP_CXX) #target_include_directories(${PROJECT_NAME} OpenMP::OpenMP_CXX) diff --git a/chronoseulersimulationmodel.cpp b/chronoseulersimulationmodel.cpp new file mode 100644 index 0000000..e63ce02 --- /dev/null +++ b/chronoseulersimulationmodel.cpp @@ -0,0 +1,175 @@ +#include "chronoseulersimulationmodel.hpp" +#include +#include +#include +#include +#include +#include + +using namespace chrono; +using namespace chrono::fea; +std::shared_ptr ChronosEulerSimulationModel::convertToChronosMesh_Euler( + 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]; + edgeMeshVertsToChronosNodes[vi] = std::make_shared( + ChFrame<>(ChVector<>(vertex.cP()[0], vertex.cP()[1], vertex.cP()[2]), + ChVector<>(0, 1, 0))); + mesh_chronos->AddNode(edgeMeshVertsToChronosNodes[vi]); + } + //add elements + ChBuilderBeamEuler 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.b; + const double beam_wy = element.dimensions.h; + const double E = element.material.youngsModulus; + const double poisson = element.material.poissonsRatio; + const double density = 1e4; + + auto msection = chrono_types::make_shared(); + 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 + 5, // the number of ChElementBeamEuler to create + edgeMeshVertsToChronosNodes[vi0], // the 'A' point in space (beginning of beam) + edgeMeshVertsToChronosNodes[vi1], // 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 + } + + return mesh_chronos; +} + +ChronosEulerSimulationModel::ChronosEulerSimulationModel() {} + +SimulationResults ChronosEulerSimulationModel::executeSimulation( + const std::shared_ptr &pJob) +{ + chrono::ChSystemSMC my_system; + my_system.SetTimestepperType(chrono::ChTimestepper::Type::EULER_IMPLICIT_LINEARIZED); + std::vector> edgeMeshVertsToChronosNodes; + auto mesh_chronos = convertToChronosMesh_Euler(pJob->pMesh, edgeMeshVertsToChronosNodes); + // mesh_chronos->SetAutomaticGravity(false); + //parse constrained vertices + for (const std::pair> &constrainedVertex : + pJob->constrainedVertices) { + const int &constrainedVi = constrainedVertex.first; + const std::unordered_set &constrainedDoF = constrainedVertex.second; + const bool vertexIsFullyConstrained = constrainedDoF.size() == 6 + && constrainedDoF.contains(0) + && constrainedDoF.contains(1) + && constrainedDoF.contains(2) + && constrainedDoF.contains(3) + && constrainedDoF.contains(4) + && constrainedDoF.contains(5); + if (vertexIsFullyConstrained) { + edgeMeshVertsToChronosNodes[constrainedVi]->SetFixed(true); + } else { + std::cerr << "Currently only rigid vertices are handled." << std::endl; + SimulationResults simulationResults; + simulationResults.converged = false; + assert(false); + return simulationResults; + } + } + // std::dynamic_pointer_cast>(mesh_chronos->GetNode(1)) + // ->SetFixed(true); + //parse external forces + for (const std::pair &externalForce : pJob->nodalExternalForces) { + const int &forceVi = externalForce.first; + const Vector6d &force = externalForce.second; + if (Eigen::Vector3d(force[3], force[4], force[5]).norm() > 1e-10) { + std::cerr << "Moments are currently not supported." << std::endl; + SimulationResults simulationResults; + simulationResults.converged = false; + assert(false); + return simulationResults; + } else { + edgeMeshVertsToChronosNodes[forceVi]->SetForce(ChVector<>(force[0], force[1], force[2])); + } + } + // edgeMeshVertsToChronosNodes[0]->SetFixed(true); + // edgeMeshVertsToChronosNodes[1]->SetForce(ChVector<>(1e-10, 0, 0)); + // edgeMeshVertsToChronosNodes[0]->SetFixed(true); + // builder.GetLastBeamNodes().front()->SetFixed(true); + // builder.GetLastBeamNodes().back()->SetForce(ChVector<>(0, 0, 1e1)); + + // edgeMeshVertsToChronosNodes[3]->SetFixed(true); + // edgeMeshVertsToChronosNodes[7]->SetFixed(true); + // edgeMeshVertsToChronosNodes[11]->SetFixed(true); + // edgeMeshVertsToChronosNodes[15]->SetForce(ChVector<>(0.0, 0, 1e-1)); + // edgeMeshVertsToChronosNodes[19]->SetForce(ChVector<>(0.0, 0, 1e-1)); + // edgeMeshVertsToChronosNodes[23]->SetForce(ChVector<>(0.0, 0, 1e-1)); + // builder.GetLastBeamNodes().back()->SetForce(ChVector<>(0.0, 1e2, 0.0)); + + my_system.Add(mesh_chronos); + + auto solver = chrono_types::make_shared(); + solver->SetMaxIterations(5000); + solver->SetTolerance(1e-10); + solver->EnableDiagonalPreconditioner(true); + solver->SetVerbose(true); + my_system.SetSolver(solver); + SimulationResults simulationResults; + simulationResults.converged = my_system.DoStaticNonlinear(); + if (!simulationResults.converged) { + std::cerr << "Simulation failed" << std::endl; + assert(false); + return simulationResults; + } + // my_system.DoStaticLinear(); + 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++) { + std::shared_ptr node_chronos = edgeMeshVertsToChronosNodes[vi]; + const auto posDisplacement = node_chronos->GetPos() - node_chronos->GetX0().GetPos(); + //Translations + simulationResults.displacements[vi][0] = posDisplacement[0]; + simulationResults.displacements[vi][1] = posDisplacement[1]; + simulationResults.displacements[vi][2] = posDisplacement[2]; + //Rotations + const 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(); + simulationResults.displacements[vi][3] = eulerAngles[0]; + simulationResults.displacements[vi][4] = eulerAngles[1]; + simulationResults.displacements[vi][5] = eulerAngles[2]; + } + 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; +} diff --git a/chronoseulersimulationmodel.hpp b/chronoseulersimulationmodel.hpp new file mode 100644 index 0000000..488cb17 --- /dev/null +++ b/chronoseulersimulationmodel.hpp @@ -0,0 +1,23 @@ +#ifndef CHRONOSEULERSIMULATIONMODEL_HPP +#define CHRONOSEULERSIMULATIONMODEL_HPP + +#include "simulationmodel.hpp" + +namespace chrono { +namespace fea { +class ChMesh; +class ChNodeFEAxyzrot; +} // namespace fea +} // namespace chrono + +class ChronosEulerSimulationModel : public SimulationModel +{ +public: + ChronosEulerSimulationModel(); + SimulationResults executeSimulation(const std::shared_ptr &pJob) override; + static std::shared_ptr convertToChronosMesh_Euler( + const std::shared_ptr &pMesh, + std::vector> &edgeMeshVertsToChronosNodes); +}; + +#endif // CHRONOSEULERSIMULATIONMODEL_HPP