diff --git a/reducedmodelevaluator.cpp b/reducedmodelevaluator.cpp index c6512d1..a1b31f2 100644 --- a/reducedmodelevaluator.cpp +++ b/reducedmodelevaluator.cpp @@ -1,8 +1,8 @@ #include "reducedmodelevaluator.hpp" -#include "chronoseulersimulationmodel.hpp" #include "hexagonremesher.hpp" #include "reducedmodel.hpp" #include "reducedmodeloptimizer.hpp" +#include "simulationmodelfactory.hpp" #include "trianglepatterngeometry.hpp" #include #include @@ -531,15 +531,18 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel( << tilledPatternResultsFolderPath << std::endl; //Full std::cout << "Executing:" << jobLabel << std::endl; + SimulationModelFactory factory; + std::unique_ptr pTilledPatternSimulationModel = factory.create( + simulationModelLabel); + //TODO: since the drm simulation model does not have a common interface with the rest of simulation models I need to cast it in order to pass simulation settings. Fix it by removing the SimulationSettings argument if (simulationModelLabel == DRMSimulationModel::label) { - DRMSimulationModel drmSimulationModel; - simulationResults_tilledPattern = drmSimulationModel - .executeSimulation(pJob_tiledFullPattern, - drmSimulationSettings); + simulationResults_tilledPattern = static_cast( + pTilledPatternSimulationModel.get()) + ->executeSimulation(pJob_tiledFullPattern, + drmSimulationSettings); } else if (simulationModelLabel == ChronosEulerSimulationModel::label) { - ChronosEulerSimulationModel chronosSimulationModel; - simulationResults_tilledPattern = chronosSimulationModel.executeSimulation( - pJob_tiledFullPattern); + simulationResults_tilledPattern + = pTilledPatternSimulationModel->executeSimulation(pJob_tiledFullPattern); } else { std::cerr << "Simulation model used for computing the optimization results was " "not recognized" diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp index 19402d2..3774873 100644 --- a/reducedmodeloptimizer.cpp +++ b/reducedmodeloptimizer.cpp @@ -13,6 +13,8 @@ #include #include //#define USE_SCENARIO_WEIGHTS +#include "chronosigasimulationmodel.hpp" +#include "simulationmodelfactory.hpp" #include using namespace ReducedModelOptimization; @@ -1499,81 +1501,81 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario( job.label = "Saddle"; } -double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude) -{ - SimulationJob job; - job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; - g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( - forceMagnitude, - g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, - job); - // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, - // optimizationState.fullPatternInterfaceViPairs, - // job); +//double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude) +//{ +// SimulationJob job; +// job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; +// g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( +// forceMagnitude, +// g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, +// job); +// // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, +// // optimizationState.fullPatternInterfaceViPairs, +// // job); - 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 - double error; - // job.pMesh->setLabel("initial"); - // job.pMesh->registerForDrawing(); - // results.registerForDrawing(); - // polyscope::show(); - std::string saveJobToPath; - if (!results.converged) { - // std::cout << "Force used:" << forceMagnitude << std::endl; - error = std::numeric_limits::max(); - // DRMSimulationModel::Settings debugSimulationSettings; - // debugSimulationSettings.isDebugMode = true; - // debugSimulationSettings.debugModeStep = 1000; - // // debugSimulationSettings.maxDRMIterations = 100000; - // debugSimulationSettings.shouldDraw = true; - // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; - // debugSimulationSettings.shouldCreatePlots = true; - // // debugSimulationSettings.Dtini = 0.06; - // debugSimulationSettings.beVerbose = true; - // debugSimulationSettings.useAverage = true; - // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; - // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; - // auto debugResults = simulator.executeSimulation(std::make_shared(job), - // debugSimulationSettings); - // std::terminate(); - saveJobToPath = "../nonConvergingJobs"; - } else { - error = std::abs( - results - .rotationalDisplacementQuaternion[g_baseScenarioMaxDisplacementHelperData - .interfaceViForComputingScenarioError] - .angularDistance(Eigen::Quaterniond::Identity()) - - desiredRotationAngle); - saveJobToPath = "../convergingJobs"; - } - // std::filesystem::path outputPath(std::filesystem::path(saveJobToPath) - // .append(job.pMesh->getLabel()) - // .append("mag_" + optimizationState.currentScenarioName)); - // std::filesystem::create_directories(outputPath); - // job.save(outputPath); - // settings.save(outputPath); +// 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 +// double error; +// // job.pMesh->setLabel("initial"); +// // job.pMesh->registerForDrawing(); +// // results.registerForDrawing(); +// // polyscope::show(); +// std::string saveJobToPath; +// if (!results.converged) { +// // std::cout << "Force used:" << forceMagnitude << std::endl; +// error = std::numeric_limits::max(); +// // DRMSimulationModel::Settings debugSimulationSettings; +// // debugSimulationSettings.isDebugMode = true; +// // debugSimulationSettings.debugModeStep = 1000; +// // // debugSimulationSettings.maxDRMIterations = 100000; +// // debugSimulationSettings.shouldDraw = true; +// // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; +// // debugSimulationSettings.shouldCreatePlots = true; +// // // debugSimulationSettings.Dtini = 0.06; +// // debugSimulationSettings.beVerbose = true; +// // debugSimulationSettings.useAverage = true; +// // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; +// // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; +// // auto debugResults = simulator.executeSimulation(std::make_shared(job), +// // debugSimulationSettings); +// // std::terminate(); +// saveJobToPath = "../nonConvergingJobs"; +// } else { +// error = std::abs( +// results +// .rotationalDisplacementQuaternion[g_baseScenarioMaxDisplacementHelperData +// .interfaceViForComputingScenarioError] +// .angularDistance(Eigen::Quaterniond::Identity()) +// - desiredRotationAngle); +// saveJobToPath = "../convergingJobs"; +// } +// // std::filesystem::path outputPath(std::filesystem::path(saveJobToPath) +// // .append(job.pMesh->getLabel()) +// // .append("mag_" + optimizationState.currentScenarioName)); +// // std::filesystem::create_directories(outputPath); +// // job.save(outputPath); +// // settings.save(outputPath); -#ifdef POLYSCOPE_DEFINED - std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl; -#endif - return error; -} +//#ifdef POLYSCOPE_DEFINED +// std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl; +//#endif +// return error; +//} void search(const std::function objectiveFunction, const double &targetY, @@ -1603,68 +1605,68 @@ void search(const std::function objectiveFunction, } } -double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude) -{ - SimulationJob job; - job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; - g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( - forceMagnitude, - g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, - job); - // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, - // optimizationState.fullPatternInterfaceViPairs, - // job); +//double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude) +//{ +// SimulationJob job; +// job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; +// g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( +// forceMagnitude, +// g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, +// job); +// // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, +// // optimizationState.fullPatternInterfaceViPairs, +// // job); - ChronosEulerSimulationModel simulator; - SimulationResults results = simulator.executeSimulation(std::make_shared(job)); +// 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 - double error; - if (!results.converged) { - error = std::numeric_limits::max(); - std::filesystem::path outputPath( - std::filesystem::path("../nonConvergingJobs") - .append(job.pMesh->getLabel()) - .append("mag_" + g_baseScenarioMaxDisplacementHelperData.currentScenarioName)); - std::filesystem::create_directories(outputPath); - job.save(outputPath); - // std::terminate(); - } else { - error = std::abs(results - .displacements[g_baseScenarioMaxDisplacementHelperData - .interfaceViForComputingScenarioError] - .getTranslation() - .norm() - - desiredDisplacementValue); - } +// // 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 +// double error; +// if (!results.converged) { +// error = std::numeric_limits::max(); +// std::filesystem::path outputPath( +// std::filesystem::path("../nonConvergingJobs") +// .append(job.pMesh->getLabel()) +// .append("mag_" + g_baseScenarioMaxDisplacementHelperData.currentScenarioName)); +// std::filesystem::create_directories(outputPath); +// job.save(outputPath); +// // std::terminate(); +// } else { +// error = std::abs(results +// .displacements[g_baseScenarioMaxDisplacementHelperData +// .interfaceViForComputingScenarioError] +// .getTranslation() +// .norm() +// - desiredDisplacementValue); +// } -#ifdef POLYSCOPE_DEFINED - std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl; -#endif +//#ifdef POLYSCOPE_DEFINED +// std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl; +//#endif - return error; -} +// return error; +//} #ifdef USE_ENSMALLEN struct ForceMagnitudeOptimization @@ -2106,7 +2108,7 @@ void ReducedModelOptimizer::optimize( // simulationSettings.isDebugMode = true; // simulationSettings.debugModeStep = 100000; //#ifdef POLYSCOPE_DEFINED - constexpr bool drawFullPatternSimulationResults = false; + constexpr bool drawFullPatternSimulationResults = true; auto t1 = std::chrono::high_resolution_clock::now(); //#endif results.wasSuccessful = true; @@ -2115,12 +2117,12 @@ void ReducedModelOptimizer::optimize( std::for_each( //#ifndef POLYSCOPE_DEFINED - std::execution::par_unseq, + // std::execution::par_unseq, //#endif optimizationState.simulationScenarioIndices.begin(), optimizationState.simulationScenarioIndices.end(), [&](const int &simulationScenarioIndex) { - const std::shared_ptr &pFullPatternSimulationJob + const std::shared_ptr &pPatternSimulationJob = optimizationState.fullPatternSimulationJobs[simulationScenarioIndex]; // std::cout << "Executing " << pFullPatternSimulationJob->getLabel() << std::endl; @@ -2128,31 +2130,45 @@ void ReducedModelOptimizer::optimize( std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) .append("FullPatternResults") .append(m_pFullPatternSimulationMesh->getLabel()) - .append(pFullPatternSimulationJob->getLabel())); + .append(pPatternSimulationJob->getLabel())); + std::unique_ptr pPatternSimulator = SimulationModelFactory::create( + optimizationSettings.simulationModelLabel); + pPatternSimulator->setStructure(m_pFullPatternSimulationMesh); + SimulationResults patternSimulationResults = pPatternSimulator->executeSimulation( + pPatternSimulationJob); - ChronosEulerSimulationModel patternSimulator; - patternSimulator.setStructure(m_pFullPatternSimulationMesh); - SimulationResults results = patternSimulator.executeSimulation( - pFullPatternSimulationJob); + ChronosIGASimulationModel debug_igaSimulationModel; + const SimulationResults debug_igaPatternSimulationResults + = debug_igaSimulationModel.executeSimulation(pPatternSimulationJob); - if (!results.converged) { + double error; + if (!patternSimulationResults + .isEqual(Utilities::toEigenMatrix( + debug_igaPatternSimulationResults.displacements), + error)) { + std::cerr + << "The computed ground truth does not match the IGA model results. Error is:" + << error << std::endl; + } + + if (!patternSimulationResults.converged) { //#ifdef POLYSCOPE_DEFINED std::filesystem::path outputPath( std::filesystem::path("../nonConvergingJobs") .append(m_pFullPatternSimulationMesh->getLabel()) - .append("final_" + pFullPatternSimulationJob->getLabel())); + .append("final_" + pPatternSimulationJob->getLabel())); std::filesystem::create_directories(outputPath); - pFullPatternSimulationJob->save(outputPath); + pPatternSimulationJob->save(outputPath); std::cerr << m_pFullPatternSimulationMesh->getLabel() << " did not converge." << std::endl; assert(false); //#endif } - optimizationState.fullPatternResults[simulationScenarioIndex] = results; + optimizationState.fullPatternResults[simulationScenarioIndex] = patternSimulationResults; SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedModelSimulationMesh; - computeReducedModelSimulationJob(*pFullPatternSimulationJob, + computeReducedModelSimulationJob(*pPatternSimulationJob, m_fullToReducedInterfaceViMap, reducedPatternSimulationJob); optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] @@ -2162,11 +2178,11 @@ void ReducedModelOptimizer::optimize( if (drawFullPatternSimulationResults) { optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( ReducedModelOptimization::Colors::patternInitial); - pFullPatternSimulationJob->registerForDrawing(optimizationState - .fullPatternSimulationJobs[0] - ->pMesh->getLabel(), - true); - results.registerForDrawing(std::array{0, 255, 0}, true); + pPatternSimulationJob->registerForDrawing(optimizationState + .fullPatternSimulationJobs[0] + ->pMesh->getLabel(), + true); + patternSimulationResults.registerForDrawing(std::array{0, 255, 0}, true); // SimulationResults fullPatternResults_linear // = linearSimulator.executeSimulation(pFullPatternSimulationJob); // patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); @@ -2174,7 +2190,7 @@ void ReducedModelOptimizer::optimize( // fullPatternResults_linear // .registerForDrawing(std::array{0, 0, 255}, true); polyscope::show(); - results.unregister(); + patternSimulationResults.unregister(); // fullPatternResults_linear.unregister(); optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); } diff --git a/simulationmodelfactory.cpp b/simulationmodelfactory.cpp new file mode 100644 index 0000000..b0a2ec9 --- /dev/null +++ b/simulationmodelfactory.cpp @@ -0,0 +1,18 @@ +#include "simulationmodelfactory.hpp" + +SimulationModelFactory::SimulationModelFactory() {} + +std::unique_ptr SimulationModelFactory::create( + const std::string &simulationModelLabel) +{ + if (simulationModelLabel == DRMSimulationModel::label) { + return std::make_unique(); + } else if (simulationModelLabel == ChronosEulerSimulationModel::label) { + return std::make_unique(); + } + std::cerr << "Simulation model used for computing the optimization results was " + "not recognized" + << std::endl; + assert(false); + return nullptr; +} diff --git a/simulationmodelfactory.hpp b/simulationmodelfactory.hpp new file mode 100644 index 0000000..543a011 --- /dev/null +++ b/simulationmodelfactory.hpp @@ -0,0 +1,17 @@ +#ifndef SIMULATIONMODELFACTORY_HPP +#define SIMULATIONMODELFACTORY_HPP + +#include "chronoseulersimulationmodel.hpp" +#include "der_leimer.hpp" +#include "drmsimulationmodel.hpp" +#include "linearsimulationmodel.hpp" +#include + +class SimulationModelFactory +{ +public: + SimulationModelFactory(); + static std::unique_ptr create(const std::string &simulationModelLabel); +}; + +#endif // SIMULATIONMODELFACTORY_HPP