Added simulation model factory class. refactoring

This commit is contained in:
iasonmanolas 2022-07-22 12:15:52 +03:00
parent 3092809d02
commit 6762c266c0
4 changed files with 212 additions and 158 deletions

View File

@ -1,8 +1,8 @@
#include "reducedmodelevaluator.hpp" #include "reducedmodelevaluator.hpp"
#include "chronoseulersimulationmodel.hpp"
#include "hexagonremesher.hpp" #include "hexagonremesher.hpp"
#include "reducedmodel.hpp" #include "reducedmodel.hpp"
#include "reducedmodeloptimizer.hpp" #include "reducedmodeloptimizer.hpp"
#include "simulationmodelfactory.hpp"
#include "trianglepatterngeometry.hpp" #include "trianglepatterngeometry.hpp"
#include <execution> #include <execution>
#include <filesystem> #include <filesystem>
@ -531,15 +531,18 @@ ReducedModelEvaluator::Results ReducedModelEvaluator::evaluateReducedModel(
<< tilledPatternResultsFolderPath << std::endl; << tilledPatternResultsFolderPath << std::endl;
//Full //Full
std::cout << "Executing:" << jobLabel << std::endl; std::cout << "Executing:" << jobLabel << std::endl;
SimulationModelFactory factory;
std::unique_ptr<SimulationModel> 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) { if (simulationModelLabel == DRMSimulationModel::label) {
DRMSimulationModel drmSimulationModel; simulationResults_tilledPattern = static_cast<DRMSimulationModel *>(
simulationResults_tilledPattern = drmSimulationModel pTilledPatternSimulationModel.get())
.executeSimulation(pJob_tiledFullPattern, ->executeSimulation(pJob_tiledFullPattern,
drmSimulationSettings); drmSimulationSettings);
} else if (simulationModelLabel == ChronosEulerSimulationModel::label) { } else if (simulationModelLabel == ChronosEulerSimulationModel::label) {
ChronosEulerSimulationModel chronosSimulationModel; simulationResults_tilledPattern
simulationResults_tilledPattern = chronosSimulationModel.executeSimulation( = pTilledPatternSimulationModel->executeSimulation(pJob_tiledFullPattern);
pJob_tiledFullPattern);
} else { } else {
std::cerr << "Simulation model used for computing the optimization results was " std::cerr << "Simulation model used for computing the optimization results was "
"not recognized" "not recognized"

View File

@ -13,6 +13,8 @@
#include <experimental/numeric> #include <experimental/numeric>
#include <functional> #include <functional>
//#define USE_SCENARIO_WEIGHTS //#define USE_SCENARIO_WEIGHTS
#include "chronosigasimulationmodel.hpp"
#include "simulationmodelfactory.hpp"
#include <armadillo> #include <armadillo>
using namespace ReducedModelOptimization; using namespace ReducedModelOptimization;
@ -1499,81 +1501,81 @@ void ReducedModelOptimizer::constructSaddleSimulationScenario(
job.label = "Saddle"; job.label = "Saddle";
} }
double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude) //double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude)
{ //{
SimulationJob job; // SimulationJob job;
job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; // job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh;
g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( // g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction(
forceMagnitude, // forceMagnitude,
g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, // g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs,
job);
// ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
// optimizationState.fullPatternInterfaceViPairs,
// job); // job);
// // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
// // optimizationState.fullPatternInterfaceViPairs,
// // job);
ChronosEulerSimulationModel simulator; // ChronosEulerSimulationModel simulator;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job)); // SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job));
// DRMSimulationModel simulator; // // DRMSimulationModel simulator;
// DRMSimulationModel::Settings settings; // // DRMSimulationModel::Settings settings;
// // settings.totalExternalForcesNormPercentageTermination = 1e-2; // // // settings.totalExternalForcesNormPercentageTermination = 1e-2;
// settings.totalTranslationalKineticEnergyThreshold = 1e-10; // // settings.totalTranslationalKineticEnergyThreshold = 1e-10;
// // settings.viscousDampingFactor = 1e-2; // // // settings.viscousDampingFactor = 1e-2;
// // settings.useKineticDamping = true; // // // settings.useKineticDamping = true;
// settings.shouldDraw = false; // // settings.shouldDraw = false;
// settings.debugModeStep = 200; // // settings.debugModeStep = 200;
// // settings.averageResidualForcesCriterionThreshold = 1e-5; // // // settings.averageResidualForcesCriterionThreshold = 1e-5;
// settings.maxDRMIterations = 200000; // // settings.maxDRMIterations = 200000;
// SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job), // // SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
// settings); // // settings);
const double &desiredRotationAngle // const double &desiredRotationAngle
= g_baseScenarioMaxDisplacementHelperData // = g_baseScenarioMaxDisplacementHelperData
.desiredMaxRotationAngle; //TODO: move from OptimizationState to a new struct // .desiredMaxRotationAngle; //TODO: move from OptimizationState to a new struct
double error; // double error;
// job.pMesh->setLabel("initial"); // // job.pMesh->setLabel("initial");
// job.pMesh->registerForDrawing(); // // job.pMesh->registerForDrawing();
// results.registerForDrawing(); // // results.registerForDrawing();
// polyscope::show(); // // polyscope::show();
std::string saveJobToPath; // std::string saveJobToPath;
if (!results.converged) { // if (!results.converged) {
// std::cout << "Force used:" << forceMagnitude << std::endl; // // std::cout << "Force used:" << forceMagnitude << std::endl;
error = std::numeric_limits<double>::max(); // error = std::numeric_limits<double>::max();
// DRMSimulationModel::Settings debugSimulationSettings; // // DRMSimulationModel::Settings debugSimulationSettings;
// debugSimulationSettings.isDebugMode = true; // // debugSimulationSettings.isDebugMode = true;
// debugSimulationSettings.debugModeStep = 1000; // // debugSimulationSettings.debugModeStep = 1000;
// // debugSimulationSettings.maxDRMIterations = 100000; // // // debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.shouldDraw = true; // // debugSimulationSettings.shouldDraw = true;
// debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; // // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// debugSimulationSettings.shouldCreatePlots = true; // // debugSimulationSettings.shouldCreatePlots = true;
// // debugSimulationSettings.Dtini = 0.06; // // // debugSimulationSettings.Dtini = 0.06;
// debugSimulationSettings.beVerbose = true; // // debugSimulationSettings.beVerbose = true;
// debugSimulationSettings.useAverage = true; // // debugSimulationSettings.useAverage = true;
// // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; // // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(std::make_shared<SimulationJob>(job), // // auto debugResults = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
// debugSimulationSettings); // // debugSimulationSettings);
// std::terminate(); // // std::terminate();
saveJobToPath = "../nonConvergingJobs"; // saveJobToPath = "../nonConvergingJobs";
} else { // } else {
error = std::abs( // error = std::abs(
results // results
.rotationalDisplacementQuaternion[g_baseScenarioMaxDisplacementHelperData // .rotationalDisplacementQuaternion[g_baseScenarioMaxDisplacementHelperData
.interfaceViForComputingScenarioError] // .interfaceViForComputingScenarioError]
.angularDistance(Eigen::Quaterniond::Identity()) // .angularDistance(Eigen::Quaterniond::Identity())
- desiredRotationAngle); // - desiredRotationAngle);
saveJobToPath = "../convergingJobs"; // saveJobToPath = "../convergingJobs";
} // }
// std::filesystem::path outputPath(std::filesystem::path(saveJobToPath) // // std::filesystem::path outputPath(std::filesystem::path(saveJobToPath)
// .append(job.pMesh->getLabel()) // // .append(job.pMesh->getLabel())
// .append("mag_" + optimizationState.currentScenarioName)); // // .append("mag_" + optimizationState.currentScenarioName));
// std::filesystem::create_directories(outputPath); // // std::filesystem::create_directories(outputPath);
// job.save(outputPath); // // job.save(outputPath);
// settings.save(outputPath); // // settings.save(outputPath);
#ifdef POLYSCOPE_DEFINED //#ifdef POLYSCOPE_DEFINED
std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl; // std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl;
#endif //#endif
return error; // return error;
} //}
void search(const std::function<double(const double &)> objectiveFunction, void search(const std::function<double(const double &)> objectiveFunction,
const double &targetY, const double &targetY,
@ -1603,68 +1605,68 @@ void search(const std::function<double(const double &)> objectiveFunction,
} }
} }
double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude) //double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude)
{ //{
SimulationJob job; // SimulationJob job;
job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh; // job.pMesh = g_baseScenarioMaxDisplacementHelperData.pFullPatternSimulationMesh;
g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction( // g_baseScenarioMaxDisplacementHelperData.constructScenarioFunction(
forceMagnitude, // forceMagnitude,
g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs, // g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs,
job);
// ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
// optimizationState.fullPatternInterfaceViPairs,
// job); // job);
// // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
// // optimizationState.fullPatternInterfaceViPairs,
// // job);
ChronosEulerSimulationModel simulator; // ChronosEulerSimulationModel simulator;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job)); // SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job));
// DRMSimulationModel simulator; // // DRMSimulationModel simulator;
// DRMSimulationModel::Settings settings; // // DRMSimulationModel::Settings settings;
// // settings.totalResidualForcesNormThreshold = 1e-3; // // // settings.totalResidualForcesNormThreshold = 1e-3;
// settings.totalExternalForcesNormPercentageTermination = 1e-2; // // settings.totalExternalForcesNormPercentageTermination = 1e-2;
// settings.totalTranslationalKineticEnergyThreshold = 1e-8; // // settings.totalTranslationalKineticEnergyThreshold = 1e-8;
// settings.viscousDampingFactor = 5e-3; // // settings.viscousDampingFactor = 5e-3;
// settings.useTranslationalKineticEnergyForKineticDamping = true; // // settings.useTranslationalKineticEnergyForKineticDamping = true;
// // settings.averageResidualForcesCriterionThreshold = 1e-5; // // // settings.averageResidualForcesCriterionThreshold = 1e-5;
// // settings.useAverage = true; // // // settings.useAverage = true;
// // settings.totalTranslationalKineticEnergyThreshold = 1e-10; // // // settings.totalTranslationalKineticEnergyThreshold = 1e-10;
// // settings.shouldUseTranslationalKineticEnergyThreshold = true; // // // settings.shouldUseTranslationalKineticEnergyThreshold = true;
// // settings.shouldDraw = true; // // // settings.shouldDraw = true;
// // settings.isDebugMode = true; // // // settings.isDebugMode = true;
// // settings.drawingStep = 200000; // // // settings.drawingStep = 200000;
// // settings.beVerbose = true; // // // settings.beVerbose = true;
// // settings.debugModeStep = 200; // // // settings.debugModeStep = 200;
// settings.maxDRMIterations = 200000; // // settings.maxDRMIterations = 200000;
// SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job), // // SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
// settings); // // settings);
const double &desiredDisplacementValue // const double &desiredDisplacementValue
= g_baseScenarioMaxDisplacementHelperData // = g_baseScenarioMaxDisplacementHelperData
.desiredMaxDisplacementValue; //TODO: move from OptimizationState to a new struct // .desiredMaxDisplacementValue; //TODO: move from OptimizationState to a new struct
double error; // double error;
if (!results.converged) { // if (!results.converged) {
error = std::numeric_limits<double>::max(); // error = std::numeric_limits<double>::max();
std::filesystem::path outputPath( // std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs") // std::filesystem::path("../nonConvergingJobs")
.append(job.pMesh->getLabel()) // .append(job.pMesh->getLabel())
.append("mag_" + g_baseScenarioMaxDisplacementHelperData.currentScenarioName)); // .append("mag_" + g_baseScenarioMaxDisplacementHelperData.currentScenarioName));
std::filesystem::create_directories(outputPath); // std::filesystem::create_directories(outputPath);
job.save(outputPath); // job.save(outputPath);
// std::terminate(); // // std::terminate();
} else { // } else {
error = std::abs(results // error = std::abs(results
.displacements[g_baseScenarioMaxDisplacementHelperData // .displacements[g_baseScenarioMaxDisplacementHelperData
.interfaceViForComputingScenarioError] // .interfaceViForComputingScenarioError]
.getTranslation() // .getTranslation()
.norm() // .norm()
- desiredDisplacementValue); // - desiredDisplacementValue);
} // }
#ifdef POLYSCOPE_DEFINED //#ifdef POLYSCOPE_DEFINED
std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl; // std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl;
#endif //#endif
return error; // return error;
} //}
#ifdef USE_ENSMALLEN #ifdef USE_ENSMALLEN
struct ForceMagnitudeOptimization struct ForceMagnitudeOptimization
@ -2106,7 +2108,7 @@ void ReducedModelOptimizer::optimize(
// simulationSettings.isDebugMode = true; // simulationSettings.isDebugMode = true;
// simulationSettings.debugModeStep = 100000; // simulationSettings.debugModeStep = 100000;
//#ifdef POLYSCOPE_DEFINED //#ifdef POLYSCOPE_DEFINED
constexpr bool drawFullPatternSimulationResults = false; constexpr bool drawFullPatternSimulationResults = true;
auto t1 = std::chrono::high_resolution_clock::now(); auto t1 = std::chrono::high_resolution_clock::now();
//#endif //#endif
results.wasSuccessful = true; results.wasSuccessful = true;
@ -2115,12 +2117,12 @@ void ReducedModelOptimizer::optimize(
std::for_each( std::for_each(
//#ifndef POLYSCOPE_DEFINED //#ifndef POLYSCOPE_DEFINED
std::execution::par_unseq, // std::execution::par_unseq,
//#endif //#endif
optimizationState.simulationScenarioIndices.begin(), optimizationState.simulationScenarioIndices.begin(),
optimizationState.simulationScenarioIndices.end(), optimizationState.simulationScenarioIndices.end(),
[&](const int &simulationScenarioIndex) { [&](const int &simulationScenarioIndex) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob const std::shared_ptr<SimulationJob> &pPatternSimulationJob
= optimizationState.fullPatternSimulationJobs[simulationScenarioIndex]; = optimizationState.fullPatternSimulationJobs[simulationScenarioIndex];
// std::cout << "Executing " << pFullPatternSimulationJob->getLabel() << std::endl; // std::cout << "Executing " << pFullPatternSimulationJob->getLabel() << std::endl;
@ -2128,31 +2130,45 @@ void ReducedModelOptimizer::optimize(
std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath)
.append("FullPatternResults") .append("FullPatternResults")
.append(m_pFullPatternSimulationMesh->getLabel()) .append(m_pFullPatternSimulationMesh->getLabel())
.append(pFullPatternSimulationJob->getLabel())); .append(pPatternSimulationJob->getLabel()));
std::unique_ptr<SimulationModel> pPatternSimulator = SimulationModelFactory::create(
optimizationSettings.simulationModelLabel);
pPatternSimulator->setStructure(m_pFullPatternSimulationMesh);
SimulationResults patternSimulationResults = pPatternSimulator->executeSimulation(
pPatternSimulationJob);
ChronosEulerSimulationModel patternSimulator; ChronosIGASimulationModel debug_igaSimulationModel;
patternSimulator.setStructure(m_pFullPatternSimulationMesh); const SimulationResults debug_igaPatternSimulationResults
SimulationResults results = patternSimulator.executeSimulation( = debug_igaSimulationModel.executeSimulation(pPatternSimulationJob);
pFullPatternSimulationJob);
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 //#ifdef POLYSCOPE_DEFINED
std::filesystem::path outputPath( std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs") std::filesystem::path("../nonConvergingJobs")
.append(m_pFullPatternSimulationMesh->getLabel()) .append(m_pFullPatternSimulationMesh->getLabel())
.append("final_" + pFullPatternSimulationJob->getLabel())); .append("final_" + pPatternSimulationJob->getLabel()));
std::filesystem::create_directories(outputPath); std::filesystem::create_directories(outputPath);
pFullPatternSimulationJob->save(outputPath); pPatternSimulationJob->save(outputPath);
std::cerr << m_pFullPatternSimulationMesh->getLabel() << " did not converge." std::cerr << m_pFullPatternSimulationMesh->getLabel() << " did not converge."
<< std::endl; << std::endl;
assert(false); assert(false);
//#endif //#endif
} }
optimizationState.fullPatternResults[simulationScenarioIndex] = results; optimizationState.fullPatternResults[simulationScenarioIndex] = patternSimulationResults;
SimulationJob reducedPatternSimulationJob; SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedModelSimulationMesh; reducedPatternSimulationJob.pMesh = m_pReducedModelSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob, computeReducedModelSimulationJob(*pPatternSimulationJob,
m_fullToReducedInterfaceViMap, m_fullToReducedInterfaceViMap,
reducedPatternSimulationJob); reducedPatternSimulationJob);
optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]
@ -2162,11 +2178,11 @@ void ReducedModelOptimizer::optimize(
if (drawFullPatternSimulationResults) { if (drawFullPatternSimulationResults) {
optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
ReducedModelOptimization::Colors::patternInitial); ReducedModelOptimization::Colors::patternInitial);
pFullPatternSimulationJob->registerForDrawing(optimizationState pPatternSimulationJob->registerForDrawing(optimizationState
.fullPatternSimulationJobs[0] .fullPatternSimulationJobs[0]
->pMesh->getLabel(), ->pMesh->getLabel(),
true); true);
results.registerForDrawing(std::array<double, 3>{0, 255, 0}, true); patternSimulationResults.registerForDrawing(std::array<double, 3>{0, 255, 0}, true);
// SimulationResults fullPatternResults_linear // SimulationResults fullPatternResults_linear
// = linearSimulator.executeSimulation(pFullPatternSimulationJob); // = linearSimulator.executeSimulation(pFullPatternSimulationJob);
// patternDrmResults.registerForDrawing(std::array<double, 3>{0, 255, 0}, true); // patternDrmResults.registerForDrawing(std::array<double, 3>{0, 255, 0}, true);
@ -2174,7 +2190,7 @@ void ReducedModelOptimizer::optimize(
// fullPatternResults_linear // fullPatternResults_linear
// .registerForDrawing(std::array<double, 3>{0, 0, 255}, true); // .registerForDrawing(std::array<double, 3>{0, 0, 255}, true);
polyscope::show(); polyscope::show();
results.unregister(); patternSimulationResults.unregister();
// fullPatternResults_linear.unregister(); // fullPatternResults_linear.unregister();
optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister();
} }

View File

@ -0,0 +1,18 @@
#include "simulationmodelfactory.hpp"
SimulationModelFactory::SimulationModelFactory() {}
std::unique_ptr<SimulationModel> SimulationModelFactory::create(
const std::string &simulationModelLabel)
{
if (simulationModelLabel == DRMSimulationModel::label) {
return std::make_unique<DRMSimulationModel>();
} else if (simulationModelLabel == ChronosEulerSimulationModel::label) {
return std::make_unique<ChronosEulerSimulationModel>();
}
std::cerr << "Simulation model used for computing the optimization results was "
"not recognized"
<< std::endl;
assert(false);
return nullptr;
}

View File

@ -0,0 +1,17 @@
#ifndef SIMULATIONMODELFACTORY_HPP
#define SIMULATIONMODELFACTORY_HPP
#include "chronoseulersimulationmodel.hpp"
#include "der_leimer.hpp"
#include "drmsimulationmodel.hpp"
#include "linearsimulationmodel.hpp"
#include <string>
class SimulationModelFactory
{
public:
SimulationModelFactory();
static std::unique_ptr<SimulationModel> create(const std::string &simulationModelLabel);
};
#endif // SIMULATIONMODELFACTORY_HPP