Refactoring.Addded multiple simulation scenarios

This commit is contained in:
iasonmanolas 2021-04-08 20:55:56 +03:00
parent 49494ccef8
commit b54662be9e
3 changed files with 387 additions and 252 deletions

View File

@ -40,21 +40,21 @@ int main(int argc, char *argv[]) {
reducedPattern.scale(0.03, interfaceNodeIndex);
// Set the optization settings
ReducedModelOptimization::xRange beamE{"E", 0.001, 1000};
ReducedModelOptimization::xRange beamA{"A", 0.001, 1000};
ReducedModelOptimization::xRange beamI2{"I2", 0.001, 1000};
ReducedModelOptimization::xRange beamI3{"I3", 0.001, 1000};
ReducedModelOptimization::xRange beamJ{"J", 0.001, 1000};
ReducedModelOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95};
ReducedModelOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
ReducedModelOptimization::Settings settings_optimization;
ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000};
ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000};
ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95};
ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0};
ReducedPatternOptimization::Settings settings_optimization;
settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3,
innerHexagonSize, innerHexagonAngle};
const bool input_numberOfFunctionCallsDefined = argc >= 4;
settings_optimization.numberOfFunctionCalls =
input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
settings_optimization.normalizationStrategy
= ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon;
= ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
settings_optimization.normalizationParameter = 0.0003;
settings_optimization.solutionAccuracy = 0.001;
@ -66,7 +66,8 @@ int main(int argc, char *argv[]) {
assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]);
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size());
ReducedModelOptimization::Results optimizationResults = optimizer.optimize(settings_optimization);
ReducedPatternOptimization::Results optimizationResults = optimizer.optimize(
settings_optimization);
// Export results
const bool input_resultDirectoryDefined = argc >= 5;
@ -105,7 +106,9 @@ int main(int argc, char *argv[]) {
csv_results << endrow;
#ifdef POLYSCOPE_DEFINED
optimizationResults.saveMeshFiles();
optimizationResults.draw();
#endif
return 0;
}

View File

@ -4,10 +4,8 @@
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
#include <dlib/global_optimization.h>
#include <dlib/optimization.h>
using namespace ReducedModelOptimization;
using namespace ReducedPatternOptimization;
struct GlobalOptimizationVariables {
// std::vector<std::vector<Vector6d>> fullPatternDisplacements;
@ -22,15 +20,12 @@ struct GlobalOptimizationVariables {
std::vector<double> objectiveValueHistory;
Eigen::VectorXd initialParameters;
std::vector<int> simulationScenarioIndices;
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
double innerHexagonInitialRotationAngle{30};
double innerHexagonInitialPos{0};
double minY{DBL_MAX};
std::vector<double> minX;
int numOfSimulationCrashes{false};
int numberOfFunctionCalls{0};
int numberOfOptimizationParameters{5};
ReducedModelOptimization::Settings optimizationSettings;
ReducedPatternOptimization::Settings optimizationSettings;
vcg::Triangle3<double> baseTriangle;
} global;
@ -41,13 +36,13 @@ double ReducedModelOptimizer::computeDisplacementError(
&reducedToFullInterfaceViMap,
const double &normalizationFactor)
{
const double rawError = computeRawDisplacementError(fullPatternDisplacements,
reducedPatternDisplacements,
reducedToFullInterfaceViMap);
const double rawError = computeRawTranslationalError(fullPatternDisplacements,
reducedPatternDisplacements,
reducedToFullInterfaceViMap);
return rawError / normalizationFactor;
}
double ReducedModelOptimizer::computeRawDisplacementError(
double ReducedModelOptimizer::computeRawTranslationalError(
const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@ -56,18 +51,14 @@ double ReducedModelOptimizer::computeRawDisplacementError(
double error = 0;
for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
const VertexIndex reducedModelVi = reducedFullViPair.first;
Eigen::Vector3d reducedVertexDisplacement(
reducedPatternDisplacements[reducedModelVi][0],
reducedPatternDisplacements[reducedModelVi][1],
reducedPatternDisplacements[reducedModelVi][2]);
const VertexIndex fullModelVi = reducedFullViPair.second;
Eigen::Vector3d fullVertexDisplacement(
fullPatternDisplacements[fullModelVi][0],
fullPatternDisplacements[fullModelVi][1],
fullPatternDisplacements[fullModelVi][2]);
Eigen::Vector3d errorVector =
fullVertexDisplacement - reducedVertexDisplacement;
error += errorVector.norm();
const Eigen::Vector3d fullPatternVertexDiplacement = fullPatternDisplacements[fullModelVi]
.getTranslation();
const Eigen::Vector3d reducedPatternVertexDiplacement
= reducedPatternDisplacements[reducedModelVi].getTranslation();
const double vertexError = (fullPatternVertexDiplacement - reducedPatternVertexDiplacement)
.norm();
error += vertexError;
}
return error;
@ -121,7 +112,7 @@ double ReducedModelOptimizer::computeError(
simulationResults_reducedPattern.rotationalDisplacementQuaternion,
reducedToFullInterfaceViMap,
normalizationFactor_rotationalDisplacement);
return translationalError + rotationalError;
return 1.5 * translationalError + 0.5 * rotationalError;
}
double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
@ -226,10 +217,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
global.minX.assign(x, x + n);
}
#ifdef POLYSCOPE_DEFINED
// if (++global.numberOfFunctionCalls % 100 == 0) {
// std::cout << "Number of function calls:" << global.numberOfFunctionCalls
// << std::endl;
// }
++global.numberOfFunctionCalls;
if (global.numberOfFunctionCalls % 100 == 0) {
std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
}
#endif
// compute error and return it
// global.objectiveValueHistory.push_back(totalError);
@ -517,9 +508,6 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
global.initialParameters(2) = mesh->elements[0].J;
global.initialParameters(3) = mesh->elements[0].I2;
global.initialParameters(4) = mesh->elements[0].I3;
global.initialParameters(optimizationParamters-2) = global.innerHexagonInitialPos;
global.innerHexagonInitialRotationAngle = 30;
global.initialParameters(optimizationParamters-1) = global.innerHexagonInitialRotationAngle;
}
void ReducedModelOptimizer::computeReducedModelSimulationJob(
@ -674,89 +662,105 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
}
}
ReducedModelOptimization::Results
ReducedModelOptimizer::runOptimization(const Settings &settings) {
ReducedPatternOptimization::Results ReducedModelOptimizer::getResults(
const dlib::function_evaluation &optimizationResult_dlib, const Settings &settings)
{
ReducedPatternOptimization::Results results;
global.objectiveValueHistory.clear();
dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
xMin(i) = settings.xRanges[i].min;
xMax(i) = settings.xRanges[i].max;
}
auto start = std::chrono::system_clock::now();
dlib::function_evaluation result;
double (*objF)(double, double, double, double, double,double,double) = &objective;
result = dlib::find_min_global(
objF, xMin, xMax,
dlib::max_function_calls(settings.numberOfFunctionCalls),
std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy);
auto end = std::chrono::system_clock::now();
auto elapsed =
std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
ReducedModelOptimization::Results results;
results.baseTriangle = global.baseTriangle;
//Number of crashes
results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
//Value of optimal objective Y
results.objectiveValue.total = optimizationResult_dlib.y;
//Optimal X values
results.optimalXNameValuePairs.resize(settings.xRanges.size());
for(int xVariableIndex=0;xVariableIndex<settings.xRanges.size();xVariableIndex++){
results.optimalXNameValuePairs[xVariableIndex]
= std::make_pair(settings.xRanges[xVariableIndex].label, global.minX[xVariableIndex]);
std::vector<double> optimalX(settings.xRanges.size());
for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) {
if (xVariableIndex < 5) {
results.optimalXNameValuePairs[xVariableIndex]
= std::make_pair(settings.xRanges[xVariableIndex].label,
global.minX[xVariableIndex]
* global.initialParameters(xVariableIndex));
} else {
//Hex size and angle are pure values (not multipliers of the initial values)
results.optimalXNameValuePairs[xVariableIndex]
= std::make_pair(settings.xRanges[xVariableIndex].label,
global.minX[xVariableIndex]);
}
assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex));
optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex);
}
// Compute obj value per simulation scenario and the raw objective value
results.rawObjectiveValue=0;
std::vector<double> optimalX(results.optimalXNameValuePairs.size());
for (int xVariableIndex = 0; xVariableIndex < global.numberOfOptimizationParameters;
xVariableIndex++) {
optimalX[xVariableIndex] = global.minX[xVariableIndex];
assert(global.minX[xVariableIndex] == result.x(xVariableIndex));
}
updateMesh(optimalX.size(), optimalX.data());
results.objectiveValuePerSimulationScenario.resize(totalNumberOfSimulationScenarios);
// results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios);
LinearSimulationModel simulator;
double totalObjectiveValue = 0;
double totalRawObjectiveValue = 0;
results.objectiveValue.totalRaw = 0;
results.objectiveValue.perSimulationScenario_translational.resize(
totalNumberOfSimulationScenarios);
results.objectiveValue.perSimulationScenario_rawTranslational.resize(
totalNumberOfSimulationScenarios);
results.objectiveValue.perSimulationScenario_rotational.resize(totalNumberOfSimulationScenarios);
results.objectiveValue.perSimulationScenario_rawRotational.resize(
totalNumberOfSimulationScenarios);
results.objectiveValue.perSimulationScenario_total.resize(totalNumberOfSimulationScenarios);
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
totalObjectiveValue += computeError(
results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex] = computeError(
global.fullPatternResults[simulationScenarioIndex],
reducedModelResults,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
const double rawTranslationalError = computeRawDisplacementError(
//Raw translational
const double rawTranslationalError = computeRawTranslationalError(
global.fullPatternResults[simulationScenarioIndex].displacements,
reducedModelResults.displacements,
global.reducedToFullInterfaceViMap);
results.objectiveValue.perSimulationScenario_rawTranslational[simulationScenarioIndex]
= rawTranslationalError;
//Raw rotational
const double rawRotationalError = computeRawRotationalError(
global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
reducedModelResults.rotationalDisplacementQuaternion,
global.reducedToFullInterfaceViMap);
results.objectiveValue.perSimulationScenario_rawRotational[simulationScenarioIndex]
= rawRotationalError;
results.rawObjectiveValue += rawTranslationalError + rawRotationalError;
const double normalizedTranslationalError
= rawTranslationalError
/ global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
const double test_normalizedTranslationError = computeDisplacementError(
//Normalized translational
const double normalizedTranslationalError = computeDisplacementError(
global.fullPatternResults[simulationScenarioIndex].displacements,
reducedModelResults.displacements,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
const double normalizedRotationalError
= rawRotationalError
/ global.rotationalDisplacementNormalizationValues[simulationScenarioIndex];
const double test_normalizedRotationalError = computeRotationalError(
results.objectiveValue.perSimulationScenario_translational[simulationScenarioIndex]
= normalizedTranslationalError;
// const double test_normalizedTranslationError = computeDisplacementError(
// global.fullPatternResults[simulationScenarioIndex].displacements,
// reducedModelResults.displacements,
// global.reducedToFullInterfaceViMap,
// global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
//Normalized rotational
const double normalizedRotationalError = computeRotationalError(
global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
reducedModelResults.rotationalDisplacementQuaternion,
global.reducedToFullInterfaceViMap,
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
assert(test_normalizedTranslationError == normalizedTranslationalError);
assert(test_normalizedRotationalError == normalizedRotationalError);
std::cout << "Simulation scenario:" << baseSimulationScenarioNames[simulationScenarioIndex]
results.objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex]
= normalizedRotationalError;
// const double test_normalizedRotationalError = computeRotationalError(
// global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
// reducedModelResults.rotationalDisplacementQuaternion,
// global.reducedToFullInterfaceViMap,
// global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
// assert(test_normalizedTranslationError == normalizedTranslationalError);
// assert(test_normalizedRotationalError == normalizedRotationalError);
std::cout << "Simulation scenario:"
<< global.reducedPatternSimulationJobs[simulationScenarioIndex]->getLabel()
<< std::endl;
std::cout << "raw translational error:" << rawTranslationalError << std::endl;
std::cout << "translation normalization value:"
@ -768,180 +772,307 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
<< global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
<< std::endl;
std::cout << "Rotational error:" << normalizedRotationalError << std::endl;
results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
= normalizedTranslationalError + normalizedRotationalError;
// results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
// = normalizedTranslationalError + normalizedRotationalError;
std::cout << "Objective value:"
<< results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
<< results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex]
<< std::endl;
// totalObjectiveValue += results.objectiveValuePerSimulationScenario[simulationScenarioIndex];
totalRawObjectiveValue += rawTranslationalError + rawRotationalError;
results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError;
std::cout << std::endl;
}
assert(result.y == global.minY);
results.objectiveValue = result.y;
auto functionReturnValue = objective(optimalX.size(), optimalX.data());
std::cout << "Total (functionReturn)" << functionReturnValue << std::endl;
std::cout << "Total:" << totalObjectiveValue << std::endl;
#ifdef POLYSCOPE_DEFINED
std::cout << "Total optimal objective value:" << results.objectiveValue << std::endl;
// std::cout << "Total raw objective value:" << global.minY << std::endl;
if (global.minY != result.y) {
std::cerr << "Global min objective is not equal to result objective" << std::endl;
}
#endif
results.time = elapsed.count() / 1000.0;
const bool printDebugInfo = false;
if (printDebugInfo) {
std::cout << "Finished optimizing." << endl;
// std::cout << "Solution x:" << endl;
// std::cout << result.x << endl;
std::cout << "Objective value:" << global.minY << endl;
std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl;
assert(global.minY == optimizationResult_dlib.y);
}
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
results.fullPatternSimulationJobs.push_back(
global.fullPatternSimulationJobs[simulationScenarioIndex]);
results.reducedPatternSimulationJobs.push_back(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
// const std::string temp = global.reducedPatternSimulationJobs[simulationScenarioIndex]
// ->pMesh->getLabel();
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel("temp");
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing();
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp);
}
return results;
}
ReducedPatternOptimization::Results ReducedModelOptimizer::runOptimization(const Settings &settings)
{
global.objectiveValueHistory.clear();
dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
xMin(i) = settings.xRanges[i].min;
xMax(i) = settings.xRanges[i].max;
}
auto start = std::chrono::system_clock::now();
dlib::function_evaluation result_dlib;
double (*objF)(double, double, double, double, double,double,double) = &objective;
result_dlib = dlib::find_min_global(objF,
xMin,
xMax,
dlib::max_function_calls(settings.numberOfFunctionCalls),
std::chrono::hours(24 * 365 * 290),
settings.solutionAccuracy);
auto end = std::chrono::system_clock::now();
auto elapsed =
std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
ReducedPatternOptimization::Results results = getResults(result_dlib, settings);
results.time = elapsed.count() / 1000.0;
return results;
}
//TODO: create a function that takes as arguments the magnitude range, a lambda that generates the simulation job and the base sim scenario index
std::vector<std::shared_ptr<SimulationJob>>
ReducedModelOptimizer::createScenarios(
const std::shared_ptr<SimulationMesh> &pMesh) {
ReducedModelOptimizer::createFullPatternSimulationScenarios(
const std::shared_ptr<SimulationMesh> &pMesh)
{
std::vector<std::shared_ptr<SimulationJob>> scenarios;
scenarios.resize(totalNumberOfSimulationScenarios);
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalForces;
//// Axial
int simulationScenarioIndex = BaseSimulationScenario::Axial;
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType forceDirection(1, 0, 0);
const auto viPair = *viPairIt;
nodalForces[viPair.first]
= Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
* 500;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
const double maxForceMagnitude_axial = 500;
const double minForceMagnitude_axial = -500;
const int numberOfSimulationScenarios_axial
= simulationScenariosResolution[BaseSimulationScenario::Axial];
const double forceMagnitudeStep_axial = numberOfSimulationScenarios_axial == 1
? maxForceMagnitude_axial
: (maxForceMagnitude_axial - minForceMagnitude_axial)
/ (numberOfSimulationScenarios_axial - 1);
const int baseSimulationScenarioIndexOffset_axial
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Axial,
0);
for (int axialSimulationScenarioIndex = 0;
axialSimulationScenarioIndex < numberOfSimulationScenarios_axial;
axialSimulationScenarioIndex++) {
const double forceMagnitude = (forceMagnitudeStep_axial * axialSimulationScenarioIndex
+ minForceMagnitude_axial);
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType forceDirection(1, 0, 0);
const auto viPair = *viPairIt;
nodalForces[viPair.first]
= Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
* forceMagnitude;
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
scenarios[baseSimulationScenarioIndexOffset_axial + axialSimulationScenarioIndex]
= std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_"
+ std::to_string(axialSimulationScenarioIndex),
fixedVertices,
nodalForces,
{}));
}
scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh,
"Axial_" + std::to_string(simulationScenarioIndex),
fixedVertices,
nodalForces,
{}));
//// Shear
simulationScenarioIndex = BaseSimulationScenario::Shear;
fixedVertices.clear();
nodalForces.clear();
// NewMethod
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType forceDirection(0, 1, 0);
const auto viPair = *viPairIt;
nodalForces[viPair.first]
= Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * 40;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
const double maxForceMagnitude_shear = 50;
const double minForceMagnitude_shear = -50;
const int numberOfSimulationScenarios_shear
= simulationScenariosResolution[BaseSimulationScenario::Shear];
const double forceMagnitudeStep_shear = numberOfSimulationScenarios_shear == 1
? maxForceMagnitude_shear
: (maxForceMagnitude_shear - minForceMagnitude_shear)
/ (numberOfSimulationScenarios_shear - 1);
const int baseSimulationScenarioIndexOffset_shear
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Shear,
0);
for (int shearSimulationScenarioIndex = 0;
shearSimulationScenarioIndex < numberOfSimulationScenarios_shear;
shearSimulationScenarioIndex++) {
fixedVertices.clear();
nodalForces.clear();
const double forceMagnitude = (forceMagnitudeStep_shear * shearSimulationScenarioIndex
+ minForceMagnitude_shear);
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType forceDirection(0, 1, 0);
const auto viPair = *viPairIt;
nodalForces[viPair.first]
= Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
* forceMagnitude;
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
scenarios[baseSimulationScenarioIndexOffset_shear + shearSimulationScenarioIndex]
= std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_"
+ std::to_string(shearSimulationScenarioIndex),
fixedVertices,
nodalForces,
{}));
}
scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
{}));
//// Bending
simulationScenarioIndex = BaseSimulationScenario::Bending;
fixedVertices.clear();
nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViPairs) {
nodalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * 0.005;
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
const double maxForceMagnitude_bending = 0.005;
const double minForceMagnitude_bending = -0.005;
const int numberOfSimulationScenarios_bending
= simulationScenariosResolution[BaseSimulationScenario::Bending];
const double forceMagnitudeStep_bending = numberOfSimulationScenarios_bending == 1
? maxForceMagnitude_bending
: (maxForceMagnitude_bending
- minForceMagnitude_bending)
/ (numberOfSimulationScenarios_bending - 1);
const int baseSimulationScenarioIndexOffset_bending
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Bending,
0);
for (int bendingSimulationScenarioIndex = 0;
bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending;
bendingSimulationScenarioIndex++) {
fixedVertices.clear();
nodalForces.clear();
const double forceMagnitude = (forceMagnitudeStep_bending * bendingSimulationScenarioIndex
+ minForceMagnitude_bending);
for (const auto &viPair : m_fullPatternOppositeInterfaceViPairs) {
nodalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude;
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
scenarios[baseSimulationScenarioIndexOffset_bending + bendingSimulationScenarioIndex]
= std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_"
+ std::to_string(bendingSimulationScenarioIndex),
fixedVertices,
nodalForces,
{}));
}
scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
{}));
//// Dome
simulationScenarioIndex = BaseSimulationScenario::Dome;
fixedVertices.clear();
nodalForces.clear();
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
const auto viPair = *viPairIt;
if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
CoordType interfaceV = (pMesh->vert[viPair.first].cP()
- pMesh->vert[viPair.second].cP());
nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceV[0],
-interfaceV[1],
0)
* 0.025;
nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceV[0],
interfaceV[1],
0)
* 0.025;
// 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 {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
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}) * 0.005;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.005;
const double maxForceMagnitude_dome = 0.025;
const double minForceMagnitude_dome = -0.025;
const int numberOfSimulationScenarios_dome
= simulationScenariosResolution[BaseSimulationScenario::Dome];
const double forceMagnitudeStep_dome = numberOfSimulationScenarios_dome == 1
? maxForceMagnitude_dome
: (maxForceMagnitude_dome - minForceMagnitude_dome)
/ (numberOfSimulationScenarios_dome - 1);
const int baseSimulationScenarioIndexOffset_dome
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Dome,
0);
for (int domeSimulationScenarioIndex = 0;
domeSimulationScenarioIndex < numberOfSimulationScenarios_dome;
domeSimulationScenarioIndex++) {
fixedVertices.clear();
nodalForces.clear();
std::unordered_map<VertexIndex, Eigen::Vector3d> nodalForcedDisplacements;
const double forceMagnitude = (forceMagnitudeStep_dome * domeSimulationScenarioIndex
+ minForceMagnitude_dome);
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
const auto viPair = *viPairIt;
CoordType interfaceVector = (pMesh->vert[viPair.first].cP()
- pMesh->vert[viPair.second].cP());
VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0))
* interfaceVector.Normalize();
if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceVector[0],
-interfaceVector[1],
0)
* std::abs(forceMagnitude);
nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
interfaceVector[1],
0)
* std::abs(forceMagnitude);
// 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 {
nodalForces[viPair.first]
= Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]})
* forceMagnitude / 5;
nodalForces[viPair.second]
= Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
* forceMagnitude / 5;
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
}
}
scenarios[baseSimulationScenarioIndexOffset_dome + domeSimulationScenarioIndex]
= std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_"
+ std::to_string(domeSimulationScenarioIndex),
fixedVertices,
nodalForces,
nodalForcedDisplacements));
}
scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
nodalForcedDisplacements));
//// Saddle
simulationScenarioIndex = BaseSimulationScenario::Saddle;
fixedVertices.clear();
nodalForces.clear();
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
const auto &viPair = *viPairIt;
CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^
CoordType(0, 0, -1).Normalize();
if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.002;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.002;
} else {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
const double maxForceMagnitude_saddle = 0.005;
const double minForceMagnitude_saddle = -0.005;
const int numberOfSimulationScenarios_saddle
= simulationScenariosResolution[BaseSimulationScenario::Saddle];
const double forceMagnitudeStep_saddle = numberOfSimulationScenarios_saddle == 1
? maxForceMagnitude_saddle
: (maxForceMagnitude_saddle
- minForceMagnitude_saddle)
/ (numberOfSimulationScenarios_saddle - 1);
const int baseSimulationScenarioIndexOffset_saddle
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Saddle,
0);
for (int saddleSimulationScenarioIndex = 0;
saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle;
saddleSimulationScenarioIndex++) {
fixedVertices.clear();
nodalForces.clear();
const double forceMagnitude = (forceMagnitudeStep_saddle * saddleSimulationScenarioIndex
+ minForceMagnitude_saddle);
for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin();
viPairIt != m_fullPatternOppositeInterfaceViPairs.end();
viPairIt++) {
const auto &viPair = *viPairIt;
CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
^ CoordType(0, 0, -1).Normalize();
if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) {
nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude;
nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude;
} else {
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.001;
nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.001;
nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
/ 2;
nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
/ 2;
}
}
scenarios[baseSimulationScenarioIndexOffset_saddle + saddleSimulationScenarioIndex]
= std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_"
+ std::to_string(saddleSimulationScenarioIndex),
fixedVertices,
nodalForces,
{}));
}
scenarios[simulationScenarioIndex] = std::make_shared<SimulationJob>(
SimulationJob(pMesh,
baseSimulationScenarioNames[simulationScenarioIndex],
fixedVertices,
nodalForces,
{}));
return scenarios;
}
@ -1050,31 +1181,32 @@ Results ReducedModelOptimizer::optimize(
global.numOfSimulationCrashes = 0;
global.numberOfFunctionCalls = 0;
global.optimizationSettings = optimizationSettings;
global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh);
global.fullPatternSimulationJobs = createFullPatternSimulationScenarios(
m_pFullPatternSimulationMesh);
// polyscope::removeAllStructures();
DRMSimulationModel::Settings simulationSettings;
simulationSettings.shouldDraw = false;
// global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
// ReducedModelOptimization::Colors::fullInitial);
// ReducedPatternOptimization::Colors::fullInitial);
// LinearSimulationModel linearSimulator;
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings);
// SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
// pFullPatternSimulationJob);
// fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
// pFullPatternSimulationJob);
// fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed,
// true,
// true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// true,
// true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// true,
// true);
// polyscope::show();
// fullPatternResults.unregister();
// fullPatternResults_linear.unregister();
// fullPatternResults_linear.unregister();
global.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
@ -1083,6 +1215,7 @@ Results ReducedModelOptimizer::optimize(
reducedPatternSimulationJob);
global.reducedPatternSimulationJobs[simulationScenarioIndex]
= std::make_shared<SimulationJob>(reducedPatternSimulationJob);
// std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl;
}
// global.fullPatternSimulationJobs[0]->pMesh->unregister();
@ -1091,12 +1224,6 @@ Results ReducedModelOptimizer::optimize(
computeObjectiveValueNormalizationFactors();
// }
Results optResults = runOptimization(optimizationSettings);
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
optResults.fullPatternSimulationJobs.push_back(
global.fullPatternSimulationJobs[simulationScenarioIndex]);
optResults.reducedPatternSimulationJobs.push_back(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
}
return optResults;
}

View File

@ -9,6 +9,8 @@
#include "reducedmodeloptimizer_structs.hpp"
#include "simulationmesh.hpp"
#include <Eigen/Dense>
#include <dlib/global_optimization.h>
#include <dlib/optimization.h>
#ifdef POLYSCOPE_DEFINED
#include "polyscope/color_management.h"
@ -28,17 +30,17 @@ class ReducedModelOptimizer
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
public:
constexpr static std::array<int, 5> simulationScenariosResolution = {1, 1, 1, 1, 1};
constexpr static std::array<int, 5> simulationScenariosResolution = {2, 2, 6, 6, 6};
inline static int totalNumberOfSimulationScenarios
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.end(),
0);
inline static int fanSize{6};
inline static VectorType patternPlaneNormal{0, 0, 1};
ReducedModelOptimization::Results optimize(
const ReducedModelOptimization::Settings &xRanges,
const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios
= std::vector<ReducedModelOptimization::BaseSimulationScenario>());
ReducedPatternOptimization::Results optimize(
const ReducedPatternOptimization::Settings &xRanges,
const std::vector<ReducedPatternOptimization::BaseSimulationScenario> &simulationScenarios
= std::vector<ReducedPatternOptimization::BaseSimulationScenario>());
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
@ -64,7 +66,7 @@ public:
double innerHexagonRotationAngle);
static double objective(double b, double r, double E);
static std::vector<std::shared_ptr<SimulationJob>> createScenarios(
static std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationScenarios(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::unordered_map<size_t, size_t> &fullPatternOppositeInterfaceViMap);
@ -85,7 +87,7 @@ public:
static void visualizeResults(
const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
const std::vector<ReducedModelOptimization::BaseSimulationScenario> &simulationScenarios,
const std::vector<ReducedPatternOptimization::BaseSimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static void registerResultsForDrawing(
@ -94,7 +96,7 @@ public:
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap);
static double computeRawDisplacementError(
static double computeRawTranslationalError(
const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
@ -136,9 +138,9 @@ private:
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static ReducedModelOptimization::Results runOptimization(
const ReducedModelOptimization::Settings &settings);
std::vector<std::shared_ptr<SimulationJob>> createScenarios(
static ReducedPatternOptimization::Results runOptimization(
const ReducedPatternOptimization::Settings &settings);
std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationScenarios(
const std::shared_ptr<SimulationMesh> &pMesh);
void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern);
void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
@ -148,6 +150,9 @@ private:
static double objective(long n, const double *x);
DRMSimulationModel simulator;
void computeObjectiveValueNormalizationFactors();
static ReducedPatternOptimization::Results getResults(
const dlib::function_evaluation &optimizationResult_dlib,
const ReducedPatternOptimization::Settings &settings);
};
void updateMesh(long n, const double *x);
#endif // REDUCEDMODELOPTIMIZER_HPP