diff --git a/chronoseulersimulationmodel.cpp b/chronoseulersimulationmodel.cpp index 8c480e9..119c939 100644 --- a/chronoseulersimulationmodel.cpp +++ b/chronoseulersimulationmodel.cpp @@ -1,11 +1,17 @@ #include "chronoseulersimulationmodel.hpp" +#include "chrono/physics/ChLoadContainer.h" #include #include +#include #include #include #include #include +#include +#include +#include + using namespace chrono; using namespace chrono::fea; std::shared_ptr ChronosEulerSimulationModel::convertToChronosMesh_Euler( @@ -19,9 +25,9 @@ std::shared_ptr ChronosEulerSimulationModel::convertToChronosMesh_Euler( //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))); + 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 @@ -36,19 +42,21 @@ std::shared_ptr ChronosEulerSimulationModel::convertToChronosMesh_Euler( 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; + // const double poisson = element.material.poissonsRatio; + const double density = 1e0; - 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); + // 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 - 5, // the number of ChElementBeamEuler to create + 1, // 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) @@ -61,88 +69,153 @@ std::shared_ptr ChronosEulerSimulationModel::convertToChronosMesh_Euler( ChronosEulerSimulationModel::ChronosEulerSimulationModel() {} -SimulationResults ChronosEulerSimulationModel::executeSimulation( - const std::shared_ptr &pJob) +//SimulationResults ChronosEulerSimulationModel::executeSimulation( +// const std::shared_ptr &pJob) +//{} + +//chrono::ChSystemSMC convertToChronosSystem(const std::shared_ptr &pJob) +//{ +// chrono::ChSystemSMC my_system; +//} + +void ChronosEulerSimulationModel::parseForces( + const std::shared_ptr &mesh_chronos, + const std::vector> &edgeMeshVertsToChronoNodes, + const std::unordered_map &nodalExternalForces) { - 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 + 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])); + } +} + +void ChronosEulerSimulationModel::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; - 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; - } + // Create a truss + auto truss = chrono_types::make_shared(); + truss->SetBodyFixed(true); + my_system.Add(truss); + const auto &constrainedChronoNode = edgeMeshVertsToChronoNodes[constrainedVi]; + // 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 ChronosEulerSimulationModel::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()); + if (!structureInitialized || wasInitializedWithDifferentStructure) { + setStructure(pJob->pMesh); + } + chrono::ChSystemSMC my_system; + // 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); + //parse forces + parseForces(mesh_chronos, edgeMeshVertsToChronoNodes, pJob->nodalExternalForces); + //parse constrained vertices + parseConstrainedVertices(pJob, edgeMeshVertsToChronoNodes, my_system); // 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)); + // 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(); 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); + // solver->SetMaxIterations(100); + // solver->SetTolerance(1e-8); + solver->EnableWarmStart(true); // IMPORTANT for convergence when using EULER_IMPLICIT_LINEARIZED + solver->EnableDiagonalPreconditioner(true); + my_system.SetSolverForceTolerance(1e-6); + solver->SetVerbose(false); + SimulationResults simulationResults; simulationResults.converged = my_system.DoStaticNonlinear(); + // simulationResults.converged = my_system.DoStaticLinear(); if (!simulationResults.converged) { std::cerr << "Simulation failed" << std::endl; assert(false); return simulationResults; } - // my_system.DoStaticLinear(); + + // 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++) { - std::shared_ptr node_chronos = edgeMeshVertsToChronosNodes[vi]; - const auto posDisplacement = node_chronos->GetPos() - node_chronos->GetX0().GetPos(); + const auto node_chronos = edgeMeshVertsToChronoNodes[vi]; + 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]; @@ -152,10 +225,13 @@ SimulationResults ChronosEulerSimulationModel::executeSimulation( 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.setLabelPrefix("chronos"); return simulationResults; // VCGEdgeMesh deformedMesh; @@ -173,3 +249,8 @@ SimulationResults ChronosEulerSimulationModel::executeSimulation( // polyscope::show(); // return simulationResults; } + +void ChronosEulerSimulationModel::setStructure(const std::shared_ptr &pMesh) +{ + mesh_chronos = convertToChronosMesh_Euler(pMesh, edgeMeshVertsToChronoNodes); +} diff --git a/chronoseulersimulationmodel.hpp b/chronoseulersimulationmodel.hpp index 488cb17..9496f55 100644 --- a/chronoseulersimulationmodel.hpp +++ b/chronoseulersimulationmodel.hpp @@ -4,6 +4,7 @@ #include "simulationmodel.hpp" namespace chrono { +class ChSystemSMC; namespace fea { class ChMesh; class ChNodeFEAxyzrot; @@ -12,12 +13,26 @@ class ChNodeFEAxyzrot; class ChronosEulerSimulationModel : public SimulationModel { + std::shared_ptr mesh_chronos; + std::vector> edgeMeshVertsToChronoNodes; + public: ChronosEulerSimulationModel(); SimulationResults executeSimulation(const std::shared_ptr &pJob) override; + void setStructure(const std::shared_ptr &pMesh) override; static std::shared_ptr convertToChronosMesh_Euler( const std::shared_ptr &pMesh, std::vector> &edgeMeshVertsToChronosNodes); + +private: + 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); }; #endif // CHRONOSEULERSIMULATIONMODEL_HPP diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 4a39a3b..651588f 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -3006,6 +3006,14 @@ currentSimulationStep > maxDRMIterations*/ #endif return results; } + +void DRMSimulationModel::setStructure(const std::shared_ptr &pMesh) +{ + std::cout << "This function is currently not implemented" << std::endl; + assert(false); + std::terminate(); +} + SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &pJob, const Settings &settings, const SimulationResults &solutionGuess) diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index a49391e..36b3f93 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -251,7 +251,8 @@ private: void applyKineticDamping(const std::shared_ptr &pJob); - virtual SimulationResults executeSimulation(const std::shared_ptr &pJob); + virtual SimulationResults executeSimulation(const std::shared_ptr &pJob) override; + void setStructure(const std::shared_ptr &pMesh) override; void reset(const std::shared_ptr &pJob); diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp index a507611..f309ca5 100644 --- a/reducedmodeloptimizer.cpp +++ b/reducedmodeloptimizer.cpp @@ -222,6 +222,7 @@ double ReducedModelOptimizer::objective(const std::vector &x, // run simulations std::vector simulationErrorsPerScenario(totalNumberOfSimulationScenarios, 0); + // ChronosEulerSimulationModel simulator; LinearSimulationModel simulator; simulator.setStructure(pReducedPatternSimulationMesh); // simulator.initialize(); @@ -389,7 +390,7 @@ void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel, ReducedModelOptimizer::createSimulationMeshes(fullModel, reducedModel, m_pFullPatternSimulationMesh, - m_pReducedPatternSimulationMesh); + m_pReducedModelSimulationMesh); } void ReducedModelOptimizer::computeMaps( @@ -833,7 +834,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv optimizationState.functions_updateReducedPatternParameter[optimizationParameterIndex]( results.optimalXNameValuePairs[optimizationParameterIndex].second, - m_pReducedPatternSimulationMesh); //NOTE:due to this call this function is not const + m_pReducedModelSimulationMesh); //NOTE:due to this call this function is not const #ifdef POLYSCOPE_DEFINED std::cout << results.optimalXNameValuePairs[optimizationParameterIndex].first << ":" << results.optimalXNameValuePairs[optimizationParameterIndex].second << " "; @@ -841,15 +842,14 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv } #ifdef POLYSCOPE_DEFINED std::cout << std::endl; - m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); + m_pReducedModelSimulationMesh->updateEigenEdgeAndVertices(); #endif - m_pReducedPatternSimulationMesh->reset(); + m_pReducedModelSimulationMesh->reset(); results.fullPatternYoungsModulus = youngsModulus; // Compute obj value per simulation scenario and the raw objective value // updateMeshFunction(optimalX); // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); //TODO:use push_back it will make the code more readable - LinearSimulationModel simulator; results.objectiveValue.totalRaw = 0; results.objectiveValue.perSimulationScenario_translational.resize( optimizationState.simulationScenarioIndices.size()); @@ -869,11 +869,12 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv results.perScenario_fullPatternPotentialEnergy.resize( optimizationState.simulationScenarioIndices.size()); + reducedModelSimulator.setStructure(m_pReducedModelSimulationMesh); for (int i = 0; i < optimizationState.simulationScenarioIndices.size(); i++) { const int simulationScenarioIndex = optimizationState.simulationScenarioIndices[i]; std::shared_ptr &pReducedJob = optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults reducedModelResults = simulator.executeSimulation(pReducedJob); + SimulationResults reducedModelResults = reducedModelSimulator.executeSimulation(pReducedJob); #ifdef POLYSCOPE_DEFINED #ifdef USE_SCENARIO_WEIGHTS @@ -973,7 +974,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv reducedModelResults.registerForDrawing(Colors::reducedDeformed); optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing( - Colors::fullDeformed); + Colors::patternDeformed); polyscope::show(); reducedModelResults.unregister(); optimizationState.fullPatternResults[simulationScenarioIndex].unregister(); @@ -1362,6 +1363,7 @@ void ReducedModelOptimizer::constructSSimulationScenario( * forceMagnitude; } } + job.label = "Axial"; } void ReducedModelOptimizer::constructShearSimulationScenario( @@ -1382,6 +1384,7 @@ void ReducedModelOptimizer::constructShearSimulationScenario( job.constrainedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } } + job.label = "Shear"; } void ReducedModelOptimizer::constructBendingSimulationScenario( @@ -1394,6 +1397,7 @@ void ReducedModelOptimizer::constructBendingSimulationScenario( job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude; job.constrainedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } + job.label = "Bending"; } /*NOTE: From the results it seems as if the forced displacements are different in the linear and in the drm @@ -1454,6 +1458,7 @@ void ReducedModelOptimizer::constructDomeSimulationScenario( job.constrainedVertices[viPair.first] = std::unordered_set{0, 1, 2}; job.constrainedVertices[viPair.second] = std::unordered_set{2}; } + job.label = "Dome"; } void ReducedModelOptimizer::constructSaddleSimulationScenario( @@ -1483,6 +1488,7 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario( * forceMagnitude / 2; } } + job.label = "Saddle"; } double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude) @@ -1497,18 +1503,20 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni // optimizationState.fullPatternInterfaceViPairs, // job); - DRMSimulationModel simulator; - DRMSimulationModel::Settings settings; - // settings.totalExternalForcesNormPercentageTermination = 1e-2; - settings.totalTranslationalKineticEnergyThreshold = 1e-10; - // settings.viscousDampingFactor = 1e-2; - // settings.useKineticDamping = true; - settings.shouldDraw = false; - settings.debugModeStep = 200; - // settings.averageResidualForcesCriterionThreshold = 1e-5; - settings.maxDRMIterations = 200000; - SimulationResults results = simulator.executeSimulation(std::make_shared(job), - settings); + ChronosEulerSimulationModel simulator; + SimulationResults results = simulator.executeSimulation(std::make_shared(job)); + // DRMSimulationModel simulator; + // DRMSimulationModel::Settings settings; + // // settings.totalExternalForcesNormPercentageTermination = 1e-2; + // settings.totalTranslationalKineticEnergyThreshold = 1e-10; + // // settings.viscousDampingFactor = 1e-2; + // // settings.useKineticDamping = true; + // settings.shouldDraw = false; + // settings.debugModeStep = 200; + // // settings.averageResidualForcesCriterionThreshold = 1e-5; + // settings.maxDRMIterations = 200000; + // SimulationResults results = simulator.executeSimulation(std::make_shared(job), + // settings); const double &desiredRotationAngle = g_baseScenarioMaxDisplacementHelperData .desiredMaxRotationAngle; //TODO: move from OptimizationState to a new struct @@ -1599,25 +1607,28 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa // optimizationState.fullPatternInterfaceViPairs, // job); - DRMSimulationModel simulator; - DRMSimulationModel::Settings settings; - // settings.totalResidualForcesNormThreshold = 1e-3; - settings.totalExternalForcesNormPercentageTermination = 1e-2; - settings.totalTranslationalKineticEnergyThreshold = 1e-8; - settings.viscousDampingFactor = 5e-3; - settings.useTranslationalKineticEnergyForKineticDamping = true; - // settings.averageResidualForcesCriterionThreshold = 1e-5; - // settings.useAverage = true; - // settings.totalTranslationalKineticEnergyThreshold = 1e-10; - // settings.shouldUseTranslationalKineticEnergyThreshold = true; - // settings.shouldDraw = true; - // settings.isDebugMode = true; - // settings.drawingStep = 200000; - // settings.beVerbose = true; - // settings.debugModeStep = 200; - settings.maxDRMIterations = 200000; - SimulationResults results = simulator.executeSimulation(std::make_shared(job), - settings); + ChronosEulerSimulationModel simulator; + SimulationResults results = simulator.executeSimulation(std::make_shared(job)); + + // DRMSimulationModel simulator; + // DRMSimulationModel::Settings settings; + // // settings.totalResidualForcesNormThreshold = 1e-3; + // settings.totalExternalForcesNormPercentageTermination = 1e-2; + // settings.totalTranslationalKineticEnergyThreshold = 1e-8; + // settings.viscousDampingFactor = 5e-3; + // settings.useTranslationalKineticEnergyForKineticDamping = true; + // // settings.averageResidualForcesCriterionThreshold = 1e-5; + // // settings.useAverage = true; + // // settings.totalTranslationalKineticEnergyThreshold = 1e-10; + // // settings.shouldUseTranslationalKineticEnergyThreshold = true; + // // settings.shouldDraw = true; + // // settings.isDebugMode = true; + // // settings.drawingStep = 200000; + // // settings.beVerbose = true; + // // settings.debugModeStep = 200; + // settings.maxDRMIterations = 200000; + // SimulationResults results = simulator.executeSimulation(std::make_shared(job), + // settings); const double &desiredDisplacementValue = g_baseScenarioMaxDisplacementHelperData .desiredMaxDisplacementValue; //TODO: move from OptimizationState to a new struct @@ -1630,7 +1641,6 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa .append("mag_" + g_baseScenarioMaxDisplacementHelperData.currentScenarioName)); std::filesystem::create_directories(outputPath); job.save(outputPath); - settings.save(outputPath); // std::terminate(); } else { error = std::abs(results @@ -2049,11 +2059,13 @@ void ReducedModelOptimizer::optimize( results.baseTriangle = baseTriangle; - DRMSimulationModel::Settings drmSettings; - drmSettings.totalExternalForcesNormPercentageTermination = 1e-5; - drmSettings.maxDRMIterations = 200000; + // DRMSimulationModel::Settings drmSettings; + // drmSettings.totalExternalForcesNormPercentageTermination = 1e-5; + // drmSettings.maxDRMIterations = 200000; + // drmSettings.totalTranslationalKineticEnergyThreshold = 1e-8; + // drmSettings.totalTranslationalKineticEnergyThreshold = 1e-9; + // drmSettings.totalTranslationalKineticEnergyThreshold = 1e-10; - drmSettings.totalTranslationalKineticEnergyThreshold = 1e-9; // drmSettings.gradualForcedDisplacementSteps = 20; // drmSettings.linearGuessForceScaleFactor = 1; // drmSettings.viscousDampingFactor = 5e-3; @@ -2061,7 +2073,6 @@ void ReducedModelOptimizer::optimize( // drmSettings.shouldDraw = true; // drmSettings.totalExternalForcesNormPercentageTermination = 1e-2; - drmSettings.totalTranslationalKineticEnergyThreshold = 1e-8; // drmSettings.viscousDampingFactor = 5e-3; // simulationSettings.viscousDampingFactor = 5e-3; // simulationSettings.linearGuessForceScaleFactor = 1; @@ -2109,152 +2120,44 @@ void ReducedModelOptimizer::optimize( .append("FullPatternResults") .append(m_pFullPatternSimulationMesh->getLabel()) .append(pFullPatternSimulationJob->getLabel())); - // .append(pFullPatternSimulationJob->getLabel() + ".json") - //#ifdef POLYSCOPE_DEFINED - // constexpr bool recomputeFullPatternResults = true; - //#else - // constexpr bool recomputeFullPatternResults = true; - //#endif - SimulationResults patternDrmResults; - // if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) { - // fullPatternResults - // .load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), - // pFullPatternSimulationJob); - // } else { - // const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( - // patternMaxForceMagnitudesFilePath); - // if (fullPatternScenarioMagnitudesExist) { - // nlohmann::json json; - // std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); - // ifs >> json; - // fullPatternSimulationScenarioMaxMagnitudes - // = static_cast>>( - // json.at("maxMagn")); - // const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() - // != desiredBaseSimulationScenarioIndices.size(); - // if (!shouldRecompute) { - // return fullPatternSimulationScenarioMaxMagnitudes; - // } - // } - //#ifdef POLYSCOPE_DEFINED - // LinearSimulationModel linearSimulator; - // SimulationResults fullPatternResults = linearSimulator.executeSimulation( - // pFullPatternSimulationJob); - //#else - DRMSimulationModel simulator; - patternDrmResults = simulator.executeSimulation(pFullPatternSimulationJob, drmSettings); - // fullPatternResults.save(jobResultsDirectoryPath); - // } - //#endif - // if (!fullPatternResults.converged) { - // DRMSimulationModel::Settings simulationSettings_secondRound; - // simulationSettings_secondRound.viscousDampingFactor = 2e-3; - // simulationSettings_secondRound.useKineticDamping = true; - // simulationSettings.maxDRMIterations = 200000; - // fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, - // simulationSettings_secondRound); + ChronosEulerSimulationModel patternSimulator; + patternSimulator.setStructure(m_pFullPatternSimulationMesh); + SimulationResults results = patternSimulator.executeSimulation( + pFullPatternSimulationJob); - if (!patternDrmResults.converged) { -#ifdef POLYSCOPE_DEFINED + if (!results.converged) { + //#ifdef POLYSCOPE_DEFINED std::filesystem::path outputPath( std::filesystem::path("../nonConvergingJobs") .append(m_pFullPatternSimulationMesh->getLabel()) .append("final_" + pFullPatternSimulationJob->getLabel())); std::filesystem::create_directories(outputPath); pFullPatternSimulationJob->save(outputPath); -// drmSettings_secondTry.save(outputPath); -#endif - //#ifdef POLYSCOPE_DEFINED - // SimulationResults fullPatternResults_linear; - // if (drawFullPatternSimulationResults) { - // LinearSimulationModel linearSimulator; - // fullPatternResults_linear = linearSimulator.executeSimulation( - // pFullPatternSimulationJob); - // fullPatternResults_linear.labelPrefix += "_linear"; - // fullPatternResults_linear.registerForDrawing(std::array{0, 0, 255}, - // true); - // std::cout << pFullPatternSimulationJob->getLabel() - // << " did not converge on the first try. Trying again.." << std::endl; - // } + std::cerr << m_pFullPatternSimulationMesh->getLabel() << " did not converge." + << std::endl; + assert(false); //#endif - DRMSimulationModel::Settings drmSettings_secondTry = drmSettings; - // drmSettings_secondTry.shouldDraw = true; - // drmSettings_secondTry.debugModeStep = 20000; - // drmSettings_secondTry.shouldCreatePlots = true; - // drmSettings_secondTry.beVerbose = true; - drmSettings_secondTry.linearGuessForceScaleFactor = 1; - drmSettings_secondTry.viscousDampingFactor = 4e-8; - *drmSettings_secondTry.maxDRMIterations *= 5; - drmSettings_secondTry.useTotalRotationalKineticEnergyForKineticDamping = true; - // drmSettings.totalTranslationalKineticEnergyThreshold = std::nullopt; - // drmSettings.totalExternalForcesNormPercentageTermination = 1e-3; - drmSettings.totalTranslationalKineticEnergyThreshold = 1e-10; - // drmSettings_secondTry.gamma = 2; - // drmSettings_secondTry.xi = 0.98; - SimulationResults drmResults_secondTry - = simulator.executeSimulation(pFullPatternSimulationJob, drmSettings_secondTry); - - if (drmResults_secondTry.converged) { -#ifdef POLYSCOPE_DEFINED - std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel() - << " converged after a second try." << std::endl; -#endif - patternDrmResults = drmResults_secondTry; - } else { - results.wasSuccessful = false; - std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() - << " of pattern " << pFullPatternSimulationJob->pMesh->getLabel() - << " did not converge." << std::endl; -#ifdef POLYSCOPE_DEFINED - std::filesystem::path outputPath( - std::filesystem::path("../nonConvergingJobs") - .append(m_pFullPatternSimulationMesh->getLabel()) - .append("final_" + pFullPatternSimulationJob->getLabel())); - std::filesystem::create_directories(outputPath); - pFullPatternSimulationJob->save(outputPath); - drmSettings_secondTry.save(outputPath); - - if (drawFullPatternSimulationResults) { - optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( - ReducedModelOptimization::Colors::fullInitial); - pFullPatternSimulationJob->registerForDrawing( - optimizationState.fullPatternSimulationJobs[0]->pMesh->getLabel(), true); - patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); - // SimulationResults fullPatternResults_linear - // = linearSimulator.executeSimulation(pFullPatternSimulationJob); - // patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); - // fullPatternResults_linear.labelPrefix += "_linear"; - // fullPatternResults_linear - // .registerForDrawing(std::array{0, 0, 255}, true); - polyscope::show(); - patternDrmResults.unregister(); - // fullPatternResults_linear.unregister(); - optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); - } -#endif - return; - } } - optimizationState.fullPatternResults[simulationScenarioIndex] = patternDrmResults; + optimizationState.fullPatternResults[simulationScenarioIndex] = results; SimulationJob reducedPatternSimulationJob; - reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; + reducedPatternSimulationJob.pMesh = m_pReducedModelSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, m_fullToReducedInterfaceViMap, reducedPatternSimulationJob); optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); - // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; + std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; #ifdef POLYSCOPE_DEFINED if (drawFullPatternSimulationResults) { optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( - ReducedModelOptimization::Colors::fullInitial); + ReducedModelOptimization::Colors::patternInitial); pFullPatternSimulationJob->registerForDrawing(optimizationState .fullPatternSimulationJobs[0] ->pMesh->getLabel(), true); - patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); + results.registerForDrawing(std::array{0, 255, 0}, true); // SimulationResults fullPatternResults_linear // = linearSimulator.executeSimulation(pFullPatternSimulationJob); // patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); @@ -2262,7 +2165,7 @@ void ReducedModelOptimizer::optimize( // fullPatternResults_linear // .registerForDrawing(std::array{0, 0, 255}, true); polyscope::show(); - patternDrmResults.unregister(); + results.unregister(); // fullPatternResults_linear.unregister(); optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); } @@ -2271,9 +2174,6 @@ void ReducedModelOptimizer::optimize( auto t2 = std::chrono::high_resolution_clock::now(); auto s_int = std::chrono::duration_cast(t2 - t1); - std::cout << s_int.count() << std::endl; - results.wasSuccessful = false; - return; #ifdef POLYSCOPE_DEFINED std::cout << "Computed ground of truth pattern results in:" << s_int.count() << " seconds." << std::endl; diff --git a/reducedmodeloptimizer.hpp b/reducedmodeloptimizer.hpp index 0ba6e50..2161a05 100644 --- a/reducedmodeloptimizer.hpp +++ b/reducedmodeloptimizer.hpp @@ -1,6 +1,7 @@ #ifndef REDUCEDMODELOPTIMIZER_HPP #define REDUCEDMODELOPTIMIZER_HPP +#include "chronoseulersimulationmodel.hpp" #include "csvfile.hpp" #include "drmsimulationmodel.hpp" #include "edgemesh.hpp" @@ -71,7 +72,7 @@ private: const std::vector> &, SimulationJob &)> constructScenarioFunction; - std::shared_ptr m_pReducedPatternSimulationMesh; + std::shared_ptr m_pReducedModelSimulationMesh; std::shared_ptr m_pFullPatternSimulationMesh; std::unordered_map m_fullToReducedInterfaceViMap; @@ -119,7 +120,7 @@ public: constexpr static std::array simulationScenariosResolution = {12, 12, 22, 22, 22, 22}; constexpr static std::array - baseScenarioWeights = {1, 1, 2, 3, 2}; + baseScenarioWeights = {1, 1, 2, 2, 2}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.end(), @@ -303,7 +304,8 @@ private: ReducedModelOptimization::NumberOfOptimizationVariables> &optimizationParamters); - DRMSimulationModel simulator; + ChronosEulerSimulationModel patternSimulator; + LinearSimulationModel reducedModelSimulator; void computeObjectiveValueNormalizationFactors( const ReducedModelOptimization::Settings &optimizationSettings); diff --git a/simulationmodel.hpp b/simulationmodel.hpp index 6842260..2112ecc 100644 --- a/simulationmodel.hpp +++ b/simulationmodel.hpp @@ -8,6 +8,7 @@ class SimulationModel public: virtual SimulationResults executeSimulation(const std::shared_ptr &simulationJob) = 0; + virtual void setStructure(const std::shared_ptr &pMesh) = 0; }; #endif // SIMULATIONMODEL_HPP