Refactoring. Replaced PatternGeometry with ReducedModel everywhere

This commit is contained in:
iasonmanolas 2022-05-06 16:38:50 +03:00
parent 7f11ea8a47
commit 90dc2ac081
2 changed files with 298 additions and 179 deletions

View File

@ -4,14 +4,17 @@
#include "simulationhistoryplotter.hpp"
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <armadillo>
#ifdef USE_ENSMALLEN
#include "ensmallen.hpp"
#endif
#include <atomic>
#include <chrono>
#include <execution>
#include <experimental/numeric>
#include <functional>
#include <parallel/parallel.h>
//#include <parallel/parallel.h>
//#define USE_SCENARIO_WEIGHTS
#include <armadillo>
using namespace ReducedModelOptimization;
@ -233,6 +236,7 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x,
// for (const int simulationScenarioIndex : optimizationState.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &reducedJob
= optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex];
//#ifdef POLYSCOPE_DEFINED
// std::cout << reducedJob->getLabel() << ":" << std::endl;
//#endif
@ -243,7 +247,19 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x,
= std::numeric_limits<double>::max();
optimizationState.numOfSimulationCrashes++;
#ifdef POLYSCOPE_DEFINED
std::cout << "Failed simulation" << std::endl;
std::cout << "Failed simulation with x:" << std::endl;
std::cout.precision(17);
for (int i = 0; i < x.size(); i++) {
std::cout << x[i] << " ";
}
std::cout << std::endl;
// pReducedPatternSimulationMesh->registerForDrawing();
// polyscope::show();
// pReducedPatternSimulationMesh->unregister();
// std::filesystem::create_directories("failedJob");
// reducedJob->save("failedJob");
// Results debugRes;
#endif
} else {
// const double simulationScenarioError_PE = std::abs(
@ -380,7 +396,7 @@ void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel,
void ReducedModelOptimizer::computeMaps(
const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
ReducedModel &reducedPattern,
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
@ -427,11 +443,11 @@ void ReducedModelOptimizer::computeMaps(
const size_t reducedModelBaseTriangleInterfaceVi = pu_reducedModel.remap[baseTriangleInterfaceVi];
const size_t reducedModelInterfaceVertexOffset
= reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/;
Results::applyOptimizationResults_innerHexagon(initialHexagonSize,
0,
baseTriangle,
reducedPattern);
reducedPattern.createFan({0}); //TODO: should be an input parameter from main
Results::applyOptimizationResults_reducedModel_nonFanned(initialHexagonSize,
0,
baseTriangle,
reducedPattern);
reducedPattern.createFan(); //TODO: should be an input parameter from main
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex
+ reducedModelBaseTriangleInterfaceVi]
@ -503,12 +519,11 @@ void ReducedModelOptimizer::computeMaps(
#endif
}
void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern,
PatternGeometry &reducedPattern)
void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, ReducedModel &reducedModel)
{
ReducedModelOptimizer::computeMaps(slotToNode,
fullPattern,
reducedPattern,
reducedModel,
optimizationState.reducedToFullInterfaceViMap,
m_fullToReducedInterfaceViMap,
m_fullPatternOppositeInterfaceViPairs);
@ -529,38 +544,42 @@ ReducedModelOptimizer::ReducedModelOptimizer()
scenarioIsSymmetrical[BaseSimulationScenario::Shear] = false;
constructBaseScenarioFunctions[BaseSimulationScenario::Bending]
= &constructBendingSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Bending] = false;
scenarioIsSymmetrical[BaseSimulationScenario::Bending] = true;
constructBaseScenarioFunctions[BaseSimulationScenario::Dome] = &constructDomeSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Dome] = false;
scenarioIsSymmetrical[BaseSimulationScenario::Dome] = true;
constructBaseScenarioFunctions[BaseSimulationScenario::Saddle]
= &constructSaddleSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = false;
scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true;
constructBaseScenarioFunctions[BaseSimulationScenario::S] = &constructSSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::S] = false;
scenarioIsSymmetrical[BaseSimulationScenario::S] = true;
}
void ReducedModelOptimizer::initializePatterns(
PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
ReducedModel &reducedModel,
const std::array<xRange, NumberOfOptimizationVariables> &optimizationParameters)
{
assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN());
assert(fullPattern.VN() == reducedModel.VN() && fullPattern.EN() >= reducedModel.EN());
#if POLYSCOPE_DEFINED
#ifdef POLYSCOPE_DEFINED
polyscope::removeAllStructures();
#endif
// Create copies of the input models
PatternGeometry copyFullPattern;
PatternGeometry copyReducedPattern;
ReducedModel copyReducedModel;
copyFullPattern.copy(fullPattern);
copyReducedPattern.copy(reducedPattern);
copyReducedModel.copy(reducedModel);
#ifdef POLYSCOPE_DEFINED
copyFullPattern.updateEigenEdgeAndVertices();
copyReducedPattern.updateEigenEdgeAndVertices();
baseTriangle = copyReducedPattern.getBaseTriangle();
#endif
copyReducedModel.updateEigenEdgeAndVertices();
baseTriangle = copyReducedModel.getBaseTriangle();
computeMaps(copyFullPattern, copyReducedPattern);
computeMaps(copyFullPattern, copyReducedModel);
optimizationState.fullPatternOppositeInterfaceViPairs = m_fullPatternOppositeInterfaceViPairs;
createSimulationMeshes(copyFullPattern, copyReducedPattern);
g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs
= m_fullPatternOppositeInterfaceViPairs;
createSimulationMeshes(copyFullPattern, copyReducedModel);
initializeUpdateReducedPatternFunctions();
initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
}
@ -570,9 +589,9 @@ void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern,
ReducedModelOptimization::Results &optimizationResults)
{
//scale pattern and reduced model
fullPattern.scale(optimizationSettings.targetBaseTriangleSize, interfaceNodeIndex);
fullPattern.scale(optimizationSettings.targetBaseTriangleSize);
ReducedModel reducedModel;
reducedModel.scale(optimizationSettings.targetBaseTriangleSize, interfaceNodeIndex);
reducedModel.scale(optimizationSettings.targetBaseTriangleSize);
auto start = std::chrono::system_clock::now();
optimizationResults.label = fullPattern.getLabel();
optimizationResults.baseTriangleFullPattern.copy(fullPattern);
@ -589,18 +608,23 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
{
const vcg::Triangle3<double> &baseTriangle = this->baseTriangle;
optimizationState.functions_updateReducedPatternParameter[R] =
[=](const double &newR, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
[=](const double &newR, std::shared_ptr<SimulationMesh> &pReducedModelSimulationMesh) {
// std::shared_ptr<VCGEdgeMesh> pReducedModel_edgeMesh = std::dynamic_pointer_cast<VCGEdgeMesh>(
// pReducedModelSimulationMesh);
// ReducedModel::updateFannedGeometry_R(newR, pReducedModel_edgeMesh);
const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR,
newR / 2,
newR / 2);
// const vcg::Triangle3<double> &baseTriangle = getBaseTriangle();
const CoordType hexagonBaseTriangleVertexPosition
= baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0]
+ baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1]
+ baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2];
// constexpr int fanSize = 6;
for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
pReducedModelSimulationMesh->vert[2 * rotationCounter].P()
= vcg::RotationMatrix(PatternGeometry::DefaultNormal,
vcg::math::ToRad(60.0 * rotationCounter))
* hexagonBaseTriangleVertexPosition;
@ -608,17 +632,21 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
};
optimizationState.functions_updateReducedPatternParameter[Theta] =
[](const double &newTheta, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
const CoordType baseTriangleHexagonVertexPosition
= pReducedPatternSimulationMesh->vert[0].cP();
[](const double &newTheta_degrees,
std::shared_ptr<SimulationMesh> &pReducedModelSimulationMesh) {
// std::shared_ptr<VCGEdgeMesh> pReducedModel_edgeMesh
// = std::dynamic_pointer_cast<VCGEdgeMesh>(pReducedModelSimulationMesh);
// ReducedModel::updateFannedGeometry_theta(newTheta_degrees, pReducedModel_edgeMesh);
const CoordType baseTriangleHexagonVertexPosition = pReducedModelSimulationMesh->vert[0]
.cP();
const CoordType thetaRotatedHexagonBaseTriangleVertexPosition
= vcg::RotationMatrix(PatternGeometry::DefaultNormal, vcg::math::ToRad(newTheta))
= vcg::RotationMatrix(PatternGeometry::DefaultNormal,
vcg::math::ToRad(newTheta_degrees))
* baseTriangleHexagonVertexPosition;
for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
pReducedModelSimulationMesh->vert[2 * rotationCounter].P()
= vcg::RotationMatrix(PatternGeometry::DefaultNormal,
vcg::math::ToRad(60.0 * rotationCounter))
* thetaRotatedHexagonBaseTriangleVertexPosition;
@ -814,9 +842,9 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
}
#ifdef POLYSCOPE_DEFINED
std::cout << std::endl;
m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
#endif
m_pReducedPatternSimulationMesh->reset();
results.fullPatternYoungsModulus = youngsModulus;
// Compute obj value per simulation scenario and the raw objective value
// updateMeshFunction(optimalX);
@ -844,8 +872,9 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
optimizationState.simulationScenarioIndices.size());
for (int i = 0; i < optimizationState.simulationScenarioIndices.size(); i++) {
const int simulationScenarioIndex = optimizationState.simulationScenarioIndices[i];
SimulationResults reducedModelResults = simulator.executeSimulation(
optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]);
std::shared_ptr<SimulationJob> &pReducedJob
= optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex];
SimulationResults reducedModelResults = simulator.executeSimulation(pReducedJob);
#ifdef POLYSCOPE_DEFINED
#ifdef USE_SCENARIO_WEIGHTS
@ -943,11 +972,12 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
<< std::endl;
std::cout << std::endl;
// reducedModelResults.registerForDrawing(Colors::reducedDeformed);
// optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed);
// polyscope::show();
// reducedModelResults.unregister();
// optimizationState.fullPatternResults[simulationScenarioIndex].unregister();
reducedModelResults.registerForDrawing(Colors::reducedDeformed);
optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing(
Colors::fullDeformed);
polyscope::show();
reducedModelResults.unregister();
optimizationState.fullPatternResults[simulationScenarioIndex].unregister();
#endif
}
@ -968,6 +998,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
// results.draw();
}
#ifdef DLIB_DEFINED
std::array<double, NumberOfBaseSimulationScenarios>
ReducedModelOptimizer::computeFullPatternMaxSimulationForces(
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenario) const
@ -1028,6 +1059,7 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces(
return fullPatternSimulationScenarioMaxMagnitudes;
}
#endif
void ReducedModelOptimizer::runOptimization(const Settings &settings,
ReducedModelOptimization::Results &results)
@ -1087,8 +1119,8 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
= optimizationState.parametersInitialValue[parameterIndex];
return x[xIndex] * parameterInitialValue;
}();
// std::cout << "Optimization parameter:" << parameterIndex << std::endl;
// std::cout << "New value:" << parameterNewValue << std::endl;
// std::cout << "Optimization parameter:" << parameterIndex << std::endl;
// std::cout << "New value:" << parameterNewValue << std::endl;
optimizationState
.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue,
pMesh);
@ -1141,7 +1173,7 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
#else
ens::SA optimizer;
#endif
// ens::LBestPSO optimizer;
// ens::LBestPSO optimizer;
parameterGroup_optimalResult.y = optimizer.Optimize(optimizationFunction, x);
// for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
// if (x[xIndex] > optimizationState.xMax[xIndex]) {
@ -1155,7 +1187,6 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
std::cout << "optimal x:"
<< "\n"
<< x << std::endl;
std::cout << "Minima:" << minima << std::endl;
std::cout << "optimal objective:" << optimizationFunction.Evaluate(x) << std::endl;
#endif
parameterGroup_optimalResult.x.clear();
@ -1382,36 +1413,47 @@ void ReducedModelOptimizer::constructDomeSimulationScenario(
- job.pMesh->vert[viPair.second].cP());
VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0))
* interfaceVector.Normalize();
if (viPairIt == oppositeInterfaceViPairs.begin()) {
job.nodalForcedDisplacements[viPair.first]
= Eigen::Vector3d(-interfaceVector[0],
-interfaceVector[1],
0)
* 0.001
* std::abs(
forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude?
// * std::abs(forceMagnitude / maxForceMagnitude_dome);
job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
interfaceVector[1],
0)
* 0.001 * std::abs(forceMagnitude);
// * std::abs(forceMagnitude / maxForceMagnitude_dome);
// CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// ^ CoordType(0, 0, -1).Normalize();
// nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
// * 0.0001;
// nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
// * 0.0001;
} else {
job.nodalExternalForces[viPair.first]
= Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude
/ 3;
job.nodalExternalForces[viPair.second]
= Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
* forceMagnitude / 3;
job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{2};
}
// if (viPairIt == oppositeInterfaceViPairs.begin()) {
// // job.nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceVector[0],
// // -interfaceVector[1],
// // 0)
// // * 0.01
// /** std::abs(
// forceMagnitude)*/
// // ; //NOTE:Should the forced displacement change relatively to the magnitude?
// // * std::abs(forceMagnitude / maxForceMagnitude_dome);
// // job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
// // interfaceVector[1],
// // 0)
// // * 0.01 /** std::abs(forceMagnitude)*/;
// // * std::abs(forceMagnitude / maxForceMagnitude_dome);
// // CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// // ^ CoordType(0, 0, -1).Normalize();
// // nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
// // * 0.0001;
// // nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
// // * 0.0001;
// job.constrainedVertices[viPair.first] = {0, 1, 2};
// job.constrainedVertices[viPair.second] = {2};
// job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
// job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 2};
// } else {
// job.nodalExternalForces[viPair.first]
// = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude
// / 3;
// job.nodalExternalForces[viPair.second]
// = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
// * forceMagnitude / 3;
// job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
// job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{2};
// }
job.nodalExternalForces[viPair.first]
= Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude / 3;
job.nodalExternalForces[viPair.second]
= Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) * forceMagnitude
/ 3;
job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{2};
}
}
@ -1564,7 +1606,7 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-8;
settings.viscousDampingFactor = 5e-3;
settings.useKineticDamping = true;
settings.useTranslationalKineticEnergyForKineticDamping = true;
// settings.averageResidualForcesCriterionThreshold = 1e-5;
// settings.useAverage = true;
// settings.totalTranslationalKineticEnergyThreshold = 1e-10;
@ -1616,6 +1658,7 @@ struct ForceMagnitudeOptimization
};
#endif
#ifdef DLIB_DEFINED
double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
const BaseSimulationScenario &scenario) const
{
@ -1733,6 +1776,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
return forceMagnitude;
}
#endif
void ReducedModelOptimizer::computeScenarioWeights(
const std::vector<BaseSimulationScenario> &baseSimulationScenarios,
@ -1818,7 +1862,7 @@ std::vector<std::shared_ptr<SimulationJob>> ReducedModelOptimizer::createFullPat
const double forceMagnitudeStep = numberOfSimulationScenarios == 1
? maxForceMagnitude
: (maxForceMagnitude - minForceMagnitude)
/ (numberOfSimulationScenarios);
/ (numberOfSimulationScenarios - 1);
const int baseSimulationScenarioIndexOffset
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + baseScenario,
@ -1949,11 +1993,15 @@ void ReducedModelOptimizer::optimize(
continue;
}
}
std::set<BaseSimulationScenario> set_desiredScenarios(desiredBaseSimulationScenarios.begin(),
desiredBaseSimulationScenarios.end());
std::set<BaseSimulationScenario> set_zeroScenarios(zeroMagnitudeScenarios.begin(),
zeroMagnitudeScenarios.end());
std::vector<BaseSimulationScenario> usedBaseSimulationScenarios;
std::set_difference(desiredBaseSimulationScenarios.begin(),
desiredBaseSimulationScenarios.end(),
zeroMagnitudeScenarios.begin(),
zeroMagnitudeScenarios.end(),
std::set_difference(set_desiredScenarios.begin(),
set_desiredScenarios.end(),
set_zeroScenarios.begin(),
set_zeroScenarios.end(),
std::back_inserter(usedBaseSimulationScenarios));
assert(!optimizationSettings.optimizationStrategy.empty());
@ -1979,11 +2027,16 @@ void ReducedModelOptimizer::optimize(
0));
}
// constexpr bool recomputeForceMagnitudes = true;
// constexpr bool recomputeForceMagnitudes = false;
// std::array<double, NumberOfBaseSimulationScenarios> fullPatternSimulationScenarioMaxMagnitudes
// = getFullPatternMaxSimulationForces(usedBaseSimulationScenarios,
// optimizationSettings.intermediateResultsDirectoryPath,
// recomputeForceMagnitudes);
// for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) {
// fullPatternSimulationScenarioMaxMagnitudes[baseScenario]
// = std::min(fullPatternSimulationScenarioMaxMagnitudes[baseScenario],
// optimizationSettings.baseScenarioMaxMagnitudes[baseScenario]);
// }
// optimizationSettings.baseScenarioMaxMagnitudes = fullPatternSimulationScenarioMaxMagnitudes;
std::array<double, NumberOfBaseSimulationScenarios> fullPatternSimulationScenarioMaxMagnitudes
= optimizationSettings.baseScenarioMaxMagnitudes;
@ -1997,47 +2050,60 @@ void ReducedModelOptimizer::optimize(
results.baseTriangle = baseTriangle;
DRMSimulationModel::Settings simulationSettings;
simulationSettings.totalExternalForcesNormPercentageTermination = 1e-5;
simulationSettings.maxDRMIterations = 200000;
simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-15;
// simulationSettings.viscousDampingFactor = 5e-3;
// simulationSettings.debugModeStep = 100000;
// simulationSettings.shouldDraw = true;
// simulationSettings.beVerbose = true;
// simulationSettings.useKineticDamping = true;
// simulationSettings.save(std::filesystem::path(std::string(__FILE__))
// .parent_path()
// .parent_path()
// .append("DefaultSettings")
// .append("DRMSettings")
// .append("defaultDRMSimulationSettings.json"));
DRMSimulationModel::Settings drmSettings;
drmSettings.totalExternalForcesNormPercentageTermination = 1e-5;
drmSettings.maxDRMIterations = 200000;
// drmSettings.totalTranslationalKineticEnergyThreshold = 1e-10;
drmSettings.totalTranslationalKineticEnergyThreshold = 1e-9;
// drmSettings.gradualForcedDisplacementSteps = 20;
// drmSettings.linearGuessForceScaleFactor = 1;
// drmSettings.viscousDampingFactor = 5e-3;
// drmSettings.useTotalRotationalKineticEnergyForKineticDamping = true;
// drmSettings.shouldDraw = true;
// simulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
// simulationSettings.viscousDampingFactor = 1e-3;
// simulationSettings.beVerbose = true;
// simulationSettings.shouldDraw = true;
// simulationSettings.isDebugMode = true;
// simulationSettings.debugModeStep = 100000;
#ifdef POLYSCOPE_DEFINED
// drmSettings.totalExternalForcesNormPercentageTermination = 1e-2;
drmSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
// drmSettings.viscousDampingFactor = 5e-3;
// simulationSettings.viscousDampingFactor = 5e-3;
// simulationSettings.linearGuessForceScaleFactor = 1;
// simulationSettings.gamma = 0.3;
// simulationSettings.xi = 0.9999;
// drmSettings.debugModeStep = 100;
// drmSettings.shouldDraw = true;
// drmSettings.shouldCreatePlots = true;
// drmSettings.beVerbose = true;
// simulationSettings.useKineticDamping = true;
// simulationSettings.save(std::filesystem::path(std::string(__FILE__))
// .parent_path()
// .parent_path()
// .append("DefaultSettings")
// .append("DRMSettings")
// .append("defaultDRMSimulationSettings.json"));
// simulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
// simulationSettings.viscousDampingFactor = 1e-3;
// simulationSettings.beVerbose = true;
// simulationSettings.shouldDraw = true;
// simulationSettings.isDebugMode = true;
// simulationSettings.debugModeStep = 100000;
//#ifdef POLYSCOPE_DEFINED
constexpr bool drawFullPatternSimulationResults = false;
if (drawFullPatternSimulationResults) {
optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
ReducedModelOptimization::Colors::fullInitial);
}
#endif
auto t1 = std::chrono::high_resolution_clock::now();
//#endif
results.wasSuccessful = true;
optimizationState.fullPatternResults.resize(totalNumberOfSimulationScenarios);
optimizationState.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios);
std::for_each(
//#ifndef POLYSCOPE_DEFINED
std::execution::par_unseq,
//#endif
optimizationState.simulationScenarioIndices.begin(),
optimizationState.simulationScenarioIndices.end(),
[&](const int &simulationScenarioIndex) {
// for (int simulationScenarioIndex : optimizationState.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= optimizationState.fullPatternSimulationJobs[simulationScenarioIndex];
// std::cout << "Executing " << pFullPatternSimulationJob->getLabel() << std::endl;
std::filesystem::path jobResultsDirectoryPath(
std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath)
@ -2050,7 +2116,7 @@ void ReducedModelOptimizer::optimize(
//#else
// constexpr bool recomputeFullPatternResults = true;
//#endif
SimulationResults fullPatternResults;
SimulationResults patternDrmResults;
// if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) {
// fullPatternResults
// .load(std::filesystem::path(jobResultsDirectoryPath).append("Results"),
@ -2078,76 +2144,101 @@ void ReducedModelOptimizer::optimize(
//#else
DRMSimulationModel simulator;
fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings);
// 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);
#ifdef POLYSCOPE_DEFINED
// std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
// << " used viscous damping." << std::endl;
#endif
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);
if (!fullPatternResults.converged) {
results.wasSuccessful = false;
std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel()
<< " did not converge." << std::endl;
if (!patternDrmResults.converged) {
#ifdef POLYSCOPE_DEFINED
// DRMSimulationModel::Settings debugSimulationSettings;
// debugSimulationSettings.debugModeStep = 50;
// // // debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.shouldDraw = true;
// // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// debugSimulationSettings.shouldCreatePlots = true;
// // // debugSimulationSettings.Dtini = 0.06;
// debugSimulationSettings.beVerbose = true;
// debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
// debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
// debugSimulationSettings.viscousDampingFactor = 1e-2;
// // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob,
// debugSimulationSettings);
// debugResults.setLabelPrefix("debugResults");
// debugResults.registerForDrawing();
// polyscope::show();
// debugResults.unregister();
std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs")
.append(m_pFullPatternSimulationMesh->getLabel())
.append("final_" + pFullPatternSimulationJob->getLabel()));
std::filesystem::create_directories(outputPath);
pFullPatternSimulationJob->save(outputPath);
simulationSettings.save(outputPath);
std::terminate();
return;
// 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<double, 3>{0, 0, 255},
// true);
// std::cout << pFullPatternSimulationJob->getLabel()
// << " did not converge on the first try. Trying again.." << std::endl;
// }
//#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
if (drawFullPatternSimulationResults) {
// SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
// pFullPatternSimulationJob);
fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// true,
// true);
polyscope::show();
fullPatternResults.unregister();
// fullPatternResults_linear.unregister();
}
std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
<< " converged after a second try." << std::endl;
#endif
optimizationState.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
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<double, 3>{0, 255, 0}, true);
// SimulationResults fullPatternResults_linear
// = linearSimulator.executeSimulation(pFullPatternSimulationJob);
// patternDrmResults.registerForDrawing(std::array<double, 3>{0, 255, 0}, true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear
// .registerForDrawing(std::array<double, 3>{0, 0, 255}, true);
polyscope::show();
patternDrmResults.unregister();
// fullPatternResults_linear.unregister();
optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister();
}
#endif
return;
}
}
optimizationState.fullPatternResults[simulationScenarioIndex] = patternDrmResults;
SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob,
@ -2155,13 +2246,42 @@ void ReducedModelOptimizer::optimize(
reducedPatternSimulationJob);
optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]
= std::make_shared<SimulationJob>(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);
pFullPatternSimulationJob->registerForDrawing(optimizationState
.fullPatternSimulationJobs[0]
->pMesh->getLabel(),
true);
patternDrmResults.registerForDrawing(std::array<double, 3>{0, 255, 0}, true);
// SimulationResults fullPatternResults_linear
// = linearSimulator.executeSimulation(pFullPatternSimulationJob);
// patternDrmResults.registerForDrawing(std::array<double, 3>{0, 255, 0}, true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear
// .registerForDrawing(std::array<double, 3>{0, 0, 255}, true);
polyscope::show();
patternDrmResults.unregister();
// fullPatternResults_linear.unregister();
optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister();
}
#endif
});
auto t2 = std::chrono::high_resolution_clock::now();
auto s_int = std::chrono::duration_cast<std::chrono::seconds>(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;
if (drawFullPatternSimulationResults) {
optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister();
}
optimizationState.reducedPatternSimulationJobs[0]->pMesh->updateEigenEdgeAndVertices();
#endif
// if (optimizationState.optimizationSettings.normalizationStrategy
// == Settings::NormalizationStrategy::Epsilon) {
@ -2171,5 +2291,4 @@ void ReducedModelOptimizer::optimize(
std::cout << "Running reduced model optimization" << std::endl;
#endif
runOptimization(optimizationSettings, results);
results.notes = optimizationNotes;
}

View File

@ -117,7 +117,7 @@ public:
inline static std::array<std::string, ReducedModelOptimization::NumberOfOptimizationVariables>
parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"};
constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
simulationScenariosResolution = {11, 11, 21, 21, 21, 21};
simulationScenariosResolution = {12, 12, 22, 22, 22, 22};
constexpr static std::array<double, ReducedModelOptimization::NumberOfBaseSimulationScenarios>
baseScenarioWeights = {1, 1, 2, 3, 2};
inline static int totalNumberOfSimulationScenarios
@ -150,7 +150,7 @@ public:
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh);
void computeMaps(const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
ReducedModel &reducedPattern,
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
@ -263,7 +263,7 @@ public:
ReducedModelOptimization::Results optimize(
ConstPatternGeometry &fullPattern,
const ReducedModelOptimization::Settings &optimizatifullyonSettings);
const ReducedModelOptimization::Settings &optimizationSettings);
void optimize(ConstPatternGeometry &fullPattern,
ReducedModelOptimization::Settings &optimizationSettings,
@ -285,7 +285,7 @@ private:
ReducedModelOptimization::S}));
void initializePatterns(PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
ReducedModel &reducedPattern,
const std::array<ReducedModelOptimization::xRange,
ReducedModelOptimization::NumberOfOptimizationVariables>
&optimizationParameters);
@ -295,7 +295,7 @@ private:
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
void runOptimization(const ReducedModelOptimization::Settings &settings,
ReducedModelOptimization::Results &results);
void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel);
void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
void initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh,