manolas
/
ReducedModelOptimization
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
ReducedModelOptimization/src/reducedmodeloptimizer.cpp

1808 lines
95 KiB
C++
Raw Blame History

#include "reducedmodeloptimizer.hpp"
#include "linearsimulationmodel.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
#include <functional>
using namespace ReducedPatternOptimization;
struct GlobalOptimizationVariables {
// std::vector<std::vector<Vector6d>> fullPatternDisplacements;
std::vector<SimulationResults> fullPatternResults;
std::vector<double> translationalDisplacementNormalizationValues;
std::vector<double> rotationalDisplacementNormalizationValues;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
reducedToFullInterfaceViMap;
std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
fullPatternInterfaceViPairs;
matplot::line_handle gPlotHandle;
std::vector<double> objectiveValueHistoryY;
std::vector<double> objectiveValueHistoryX;
std::vector<double> plotColors;
size_t iterationOfMinima{0};
Eigen::VectorXd initialParameters;
std::vector<int> simulationScenarioIndices;
double minY{DBL_MAX};
std::vector<double> minX;
int numOfSimulationCrashes{false};
int numberOfFunctionCalls{0};
int numberOfOptimizationParameters{5};
ReducedPatternOptimization::Settings optimizationSettings;
vcg::Triangle3<double> baseTriangle;
//Variables for finding the full pattern simulation forces
std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh;
std::function<void(const double &,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>> &,
SimulationJob &)>
constructScenarioFunction;
FullPatternVertexIndex interfaceViForComputingScenarioError;
double desiredMaxDisplacementValue;
double desiredMaxRotationAngle;
std::string currentScenarioName;
std::vector<std::function<void(std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh,
const double &newValue)>>
updateReducedPatternFunctions_material;
} global;
double ReducedModelOptimizer::computeDisplacementError(
const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor)
{
const double rawError = computeRawTranslationalError(fullPatternDisplacements,
reducedPatternDisplacements,
reducedToFullInterfaceViMap);
// std::cout << "raw trans error:" << rawError << std::endl;
// std::cout << "raw trans error:" << normalizationFactor << std::endl;
return rawError / normalizationFactor;
}
double ReducedModelOptimizer::computeRawTranslationalError(
const std::vector<Vector6d> &fullPatternDisplacements,
const std::vector<Vector6d> &reducedPatternDisplacements,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap)
{
double error = 0;
for (const auto reducedFullViPair : reducedToFullInterfaceViMap) {
const VertexIndex reducedModelVi = reducedFullViPair.first;
const VertexIndex fullModelVi = reducedFullViPair.second;
const Eigen::Vector3d fullPatternVertexDiplacement = fullPatternDisplacements[fullModelVi]
.getTranslation();
const Eigen::Vector3d reducedPatternVertexDiplacement
= reducedPatternDisplacements[reducedModelVi].getTranslation();
const double vertexError = (fullPatternVertexDiplacement - reducedPatternVertexDiplacement)
.norm();
error += vertexError;
}
return error;
}
double ReducedModelOptimizer::computeRawRotationalError(
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap)
{
double rawRotationalError = 0;
for (const auto &reducedToFullInterfaceViPair : reducedToFullInterfaceViMap) {
const double vertexRotationalError
= rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance(
rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first]);
rawRotationalError += vertexRotationalError;
}
return rawRotationalError;
}
double ReducedModelOptimizer::computeRotationalError(
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_fullPattern,
const std::vector<Eigen::Quaternion<double>> &rotatedQuaternion_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor)
{
const double rawRotationalError = computeRawRotationalError(rotatedQuaternion_fullPattern,
rotatedQuaternion_reducedPattern,
reducedToFullInterfaceViMap);
return rawRotationalError / normalizationFactor;
}
double ReducedModelOptimizer::computeError(
const SimulationResults &simulationResults_fullPattern,
const SimulationResults &simulationResults_reducedPattern,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
const double &normalizationFactor_translationalDisplacement,
const double &normalizationFactor_rotationalDisplacement)
{
const double translationalError
= computeDisplacementError(simulationResults_fullPattern.displacements,
simulationResults_reducedPattern.displacements,
reducedToFullInterfaceViMap,
normalizationFactor_translationalDisplacement);
// std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl;
// std::cout << "trans error:" << translationalError << std::endl;
const double rotationalError
= computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
simulationResults_reducedPattern.rotationalDisplacementQuaternion,
reducedToFullInterfaceViMap,
normalizationFactor_rotationalDisplacement);
return global.optimizationSettings.objectiveWeights.translational * translationalError
+ global.optimizationSettings.objectiveWeights.rotational * rotationalError;
}
//double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
// double innerHexagonSize,
// double innerHexagonRotationAngle) {
// std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle};
// return ReducedModelOptimizer::objective(x.size(), x.data());
//}
double ReducedModelOptimizer::objective(const dlib::matrix<double, 0, 1> &x)
{
// std::cout.precision(17);
// for (int i = 0; i < x.size(); i++) {
// std::cout << x(i) << " ";
// }
// std::cout << std::endl;
// std::cout << x[n - 2] << " " << x[n - 1] << std::endl;
// const Element &e =
// global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
// e.axialConstFactor << " " << e.torsionConstFactor << " "
// << e.firstBendingConstFactor << " " <<
// e.secondBendingConstFactor
// << std::endl;
const int n = x.size();
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh
= global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
function_updateReducedPattern(x, pReducedPatternSimulationMesh);
// global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
// global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
// polyscope::show();
// global.reducedPatternSimulationJobs[0]->pMesh->unregister();
// std::cout << e.axialConstFact|A,I2,I3,J,r,thetaor << " " << e.torsionConstFactor << " "
// << e.firstBendingConstFactor << " " <<
// e.secondBendingConstFactor
// << std::endl;
// run simulations
double totalError = 0;
LinearSimulationModel simulator;
// FormFinder simulator;
for (const int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &reducedJob
= global.reducedPatternSimulationJobs[simulationScenarioIndex];
//#ifdef POLYSCOPE_DEFINED
// std::cout << reducedJob->getLabel() << ":" << std::endl;
//#endif
SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob);
// std::string filename;
if (!reducedModelResults.converged) {
totalError += std::numeric_limits<double>::max();
global.numOfSimulationCrashes++;
#ifdef POLYSCOPE_DEFINED
std::cout << "Failed simulation" << std::endl;
#endif
} else {
constexpr bool usePotentialEnergy = false;
double simulationScenarioError;
if (usePotentialEnergy) {
simulationScenarioError = std::abs(
reducedModelResults.internalPotentialEnergy
- global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
} else {
simulationScenarioError = computeError(
global.fullPatternResults[simulationScenarioIndex],
reducedModelResults,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
}
//#ifdef POLYSCOPE_DEFINED
// reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
// reducedModelResults.registerForDrawing(Colors::reducedDeformed);
// global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
// global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
// Colors::fullDeformed);
// polyscope::show();
// reducedModelResults.unregister();
// global.pFullPatternSimulationMesh->unregister();
// global.fullPatternResults[simulationScenarioIndex].unregister();
//#endif
// if (global.optimizationSettings.normalizationStrategy !=
// NormalizationStrategy::Epsilon &&
// simulationScenarioError > 1) {
// std::cout << "Simulation scenario "
// <<
// simulationScenarioStrings[simulationScenarioIndex]
// << " results in an error bigger than one." <<
// std::endl;
// for (size_t parameterIndex = 0; parameterIndex < n;
// parameterIndex++) {
// std::cout << "x[" + std::to_string(parameterIndex) + "]="
// << x[parameterIndex] << std::endl;
// }
// }
// #ifdef POLYSCOPE_DEFINED
// ReducedModelOptimizer::visualizeResults(
// global.fullPatternSimulationJobs[simulationScenarioIndex],
// global.reducedPatternSimulationJobs[simulationScenarioIndex],
// global.reducedToFullInterfaceViMap, false);
// ReducedModelOptimizer::visualizeResults(
// global.fullPatternSimulationJobs[simulationScenarioIndex],
// std::make_shared<SimulationJob>(
// reducedPatternMaximumDisplacementSimulationJobs
// [simulationScenarioIndex]),
// global.reducedToFullInterfaceViMap, true);
// polyscope::removeAllStructures();
// #endif // POLYSCOPE_DEFINED
totalError += simulationScenarioError * simulationScenarioError;
}
}
// std::cout << error << std::endl;
// global.objectiveValueHistory.push_back(totalError);
// global.plotColors.push_back(10);
++global.numberOfFunctionCalls;
if (totalError < global.minY) {
global.minY = totalError;
global.iterationOfMinima = global.numberOfFunctionCalls;
// std::cout << "New best:" << totalError << std::endl;
// // global.minX.assign(x.begin(), x.begin() + n);
// std::cout.precision(17);
// for (int i = 0; i < x.size(); i++) {
// std::cout << x(i) << " ";
// }
// std::cout << std::endl;
// global.objectiveValueHistoryY.push_back(std::log(totalError));
// global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls);
// global.plotColors.push_back(0.1);
// auto xPlot = matplot::linspace(0,
// global.objectiveValueHistoryY.size(),
// global.objectiveValueHistoryY.size());
// global.gPlotHandle = matplot::scatter(global.objectiveValueHistoryX,
// global.objectiveValueHistoryY,
// 4,
// global.plotColors);
// matplot::show();
// SimulationResultsReporter::createPlot("Number of Steps",
// "Objective value",
// global.objectiveValueHistoryY);
}
#ifdef POLYSCOPE_DEFINED
if (global.optimizationSettings.numberOfFunctionCalls >= 100
&& global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100)
== 0) {
std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl;
}
#endif
// compute error and return it
return totalError;
}
void ReducedModelOptimizer::createSimulationMeshes(
PatternGeometry &fullModel, PatternGeometry &reducedModel,
std::shared_ptr<SimulationMesh> &pFullPatternSimulationMesh,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
std::cerr << "Error: A rectangular cross section is expected." << std::endl;
terminate();
}
const double ym = 1 * 1e9;
// Full pattern
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
pFullPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pFullPatternSimulationMesh->setBeamMaterial(0.3, ym);
// Reduced pattern
pReducedPatternSimulationMesh =
std::make_shared<SimulationMesh>(reducedModel);
pReducedPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pReducedPatternSimulationMesh->setBeamMaterial(0.3, ym);
}
void ReducedModelOptimizer::createSimulationMeshes(
PatternGeometry &fullModel, PatternGeometry &reducedModel) {
ReducedModelOptimizer::createSimulationMeshes(
fullModel, reducedModel, m_pFullPatternSimulationMesh,
m_pReducedPatternSimulationMesh);
}
void ReducedModelOptimizer::computeMaps(
const std::unordered_map<size_t, std::unordered_set<size_t>> &slotToNode,
PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap,
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
&fullToReducedInterfaceViMap,
std::vector<std::pair<FullPatternVertexIndex, ReducedPatternVertexIndex>>
&fullPatternOppositeInterfaceViPairs)
{
// Compute the offset between the interface nodes
const size_t interfaceSlotIndex = 4; // bottom edge
assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() &&
slotToNode.find(interfaceSlotIndex)->second.size() == 1);
// Assuming that in the bottom edge there is only one vertex which is also the
// interface
const size_t baseTriangleInterfaceVi =
*(slotToNode.find(interfaceSlotIndex)->second.begin());
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<PatternGeometry::VertexPointer> pu_fullModel;
fullPattern.deleteDanglingVertices(pu_fullModel);
const size_t fullModelBaseTriangleInterfaceVi =
pu_fullModel.remap.empty() ? baseTriangleInterfaceVi
: pu_fullModel.remap[baseTriangleInterfaceVi];
const size_t fullModelBaseTriangleVN = fullPattern.VN();
fullPattern.createFan();
const size_t duplicateVerticesPerFanPair =
fullModelBaseTriangleVN - fullPattern.VN() / 6;
const size_t fullPatternInterfaceVertexOffset =
fullModelBaseTriangleVN - duplicateVerticesPerFanPair;
// std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair <<
// std::endl;
// Save excluded edges TODO:this changes the global object. Should this be a
// function parameter?
// global.reducedPatternExludedEdges.clear();
// const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN();
// for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
// for (const size_t ei : reducedModelExcludedEdges) {
// global.reducedPatternExludedEdges.insert(
// fanIndex * reducedBaseTriangleNumberOfEdges + ei);
// }
// }
// Construct reduced->full and full->reduced interface vi map
reducedToFullInterfaceViMap.clear();
vcg::tri::Allocator<VCGEdgeMesh>::PointerUpdater<PatternGeometry::VertexPointer> pu_reducedModel;
reducedPattern.deleteDanglingVertices(pu_reducedModel);
const size_t reducedModelBaseTriangleInterfaceVi =
pu_reducedModel.remap[baseTriangleInterfaceVi];
const size_t reducedModelInterfaceVertexOffset
= reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/;
Results::applyOptimizationResults_innerHexagon(initialHexagonSize,
0,
global.baseTriangle,
reducedPattern);
reducedPattern.createFan({0}); //TODO: should be an input parameter from main
for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) {
reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex +
reducedModelBaseTriangleInterfaceVi] =
fullModelBaseTriangleInterfaceVi +
fanIndex * fullPatternInterfaceVertexOffset;
}
fullToReducedInterfaceViMap.clear();
constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap);
// Create opposite vertex map
fullPatternOppositeInterfaceViPairs.clear();
fullPatternOppositeInterfaceViPairs.reserve(fanSize / 2);
// for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) {
for (int fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) {
const size_t vi0 = fullModelBaseTriangleInterfaceVi
+ fanIndex * fullPatternInterfaceVertexOffset;
const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset;
assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN());
fullPatternOppositeInterfaceViPairs.emplace_back(std::make_pair(vi0, vi1));
}
global.fullPatternInterfaceViPairs = fullPatternOppositeInterfaceViPairs;
#if POLYSCOPE_DEFINED
const bool debugMapping = false;
if (debugMapping) {
reducedPattern.registerForDrawing();
// std::vector<glm::vec3> colors_reducedPatternExcludedEdges(
// reducedPattern.EN(), glm::vec3(0, 0, 0));
// for (const size_t ei : global.reducedPatternExludedEdges) {
// colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0);
// }
// const std::string label = reducedPattern.getLabel();
// polyscope::getCurveNetwork(label)
// ->addEdgeColorQuantity("Excluded edges",
// colors_reducedPatternExcludedEdges)
// ->setEnabled(true);
// polyscope::show();
std::vector<glm::vec3> nodeColorsOpposite(fullPattern.VN(),
glm::vec3(0, 0, 0));
for (const std::pair<size_t, size_t> oppositeVerts : fullPatternOppositeInterfaceViPairs) {
auto color = polyscope::getNextUniqueColor();
nodeColorsOpposite[oppositeVerts.first] = color;
nodeColorsOpposite[oppositeVerts.second] = color;
}
fullPattern.registerForDrawing();
polyscope::getCurveNetwork(fullPattern.getLabel())
->addNodeColorQuantity("oppositeMap", nodeColorsOpposite)
->setEnabled(true);
polyscope::show();
std::vector<glm::vec3> nodeColorsReducedToFull_reduced(reducedPattern.VN(),
glm::vec3(0, 0, 0));
std::vector<glm::vec3> nodeColorsReducedToFull_full(fullPattern.VN(),
glm::vec3(0, 0, 0));
for (size_t vi = 0; vi < reducedPattern.VN(); vi++) {
if (global.reducedToFullInterfaceViMap.contains(vi)) {
auto color = polyscope::getNextUniqueColor();
nodeColorsReducedToFull_reduced[vi] = color;
nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] =
color;
}
}
polyscope::getCurveNetwork(reducedPattern.getLabel())
->addNodeColorQuantity("reducedToFull_reduced",
nodeColorsReducedToFull_reduced)
->setEnabled(true);
polyscope::getCurveNetwork(fullPattern.getLabel())
->addNodeColorQuantity("reducedToFull_full",
nodeColorsReducedToFull_full)
->setEnabled(true);
polyscope::show();
}
#endif
}
void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern,
PatternGeometry &reducedPattern)
{
ReducedModelOptimizer::computeMaps(slotToNode,
fullPattern,
reducedPattern,
global.reducedToFullInterfaceViMap,
m_fullToReducedInterfaceViMap,
m_fullPatternOppositeInterfaceViPairs);
}
ReducedModelOptimizer::ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot)
{
FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot);
FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
constructBaseScenarioFunctions.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios);
scenarioIsSymmetrical.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios);
constructBaseScenarioFunctions[BaseSimulationScenario::Axial] = &constructAxialSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Axial] = false;
constructBaseScenarioFunctions[BaseSimulationScenario::Shear] = &constructShearSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Shear] = false;
constructBaseScenarioFunctions[BaseSimulationScenario::Bending]
= &constructBendingSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Bending] = true;
constructBaseScenarioFunctions[BaseSimulationScenario::Dome] = &constructDomeSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Dome] = true;
constructBaseScenarioFunctions[BaseSimulationScenario::Saddle]
= &constructSaddleSimulationScenario;
scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true;
}
void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
const std::vector<xRange> &optimizationParameters)
{
assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN());
fullPatternNumberOfEdges = fullPattern.EN();
#if POLYSCOPE_DEFINED
polyscope::removeAllStructures();
#endif
// Create copies of the input models
PatternGeometry copyFullPattern;
PatternGeometry copyReducedPattern;
copyFullPattern.copy(fullPattern);
copyReducedPattern.copy(reducedPattern);
global.baseTriangle = copyReducedPattern.getBaseTriangle();
computeMaps(copyFullPattern, copyReducedPattern);
createSimulationMeshes(copyFullPattern, copyReducedPattern);
initializeUpdateReducedPatternFunctions();
initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters);
}
void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
{
function_updateReducedPattern_geometry = [&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh>
&pReducedPatternSimulationMesh) {
const int n = x.size();
assert(n >= 2);
// CoordType center_barycentric(1, 0, 0);
// CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5);
// CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric)
// * x[n - 2]);
CoordType movableVertex_barycentric(1 - x(n - 2), x(n - 2) / 2, x(n - 2) / 2);
CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0)
* movableVertex_barycentric[0]
+ global.baseTriangle.cP(1)
* movableVertex_barycentric[1]
+ global.baseTriangle.cP(2)
* movableVertex_barycentric[2];
baseTriangleMovableVertexPosition
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(x(n - 1)))
* baseTriangleMovableVertexPosition;
for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(60.0 * rotationCounter))
* baseTriangleMovableVertexPosition;
}
};
function_updateReducedPattern_material_E =
[&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newE) {
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE));
}
};
function_updateReducedPattern_material_A =
[&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newA) {
const double beamWidth = std::sqrt(newA);
const double beamHeight = beamWidth;
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight));
}
};
function_updateReducedPattern_material_I =
[&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI) {
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.inertia.I2 = newI;
e.inertia.I3 = newI;
}
};
function_updateReducedPattern_material_I2 =
[&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI2) {
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.inertia.I2 = newI2;
}
};
function_updateReducedPattern_material_I3 =
[&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newI3) {
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.inertia.I3 = newI3;
}
};
function_updateReducedPattern_material_J =
[&](std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh, const double &newJ) {
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.inertia.J = newJ;
}
};
}
void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh, const std::vector<xRange> &optimizationParameters)
{
global.numberOfOptimizationParameters = optimizationParameters.size();
global.initialParameters.resize(global.numberOfOptimizationParameters);
for (int optimizationParameterIndex = 0;
optimizationParameterIndex < optimizationParameters.size();
optimizationParameterIndex++) {
const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex];
if (xOptimizationParameter.label == "E") {
global.initialParameters(optimizationParameterIndex) = mesh->elements[0]
.material.youngsModulus;
global.updateReducedPatternFunctions_material.push_back(
function_updateReducedPattern_material_E);
} else if (xOptimizationParameter.label == "A") {
global.initialParameters(optimizationParameterIndex) = mesh->elements[0].A;
global.updateReducedPatternFunctions_material.push_back(
function_updateReducedPattern_material_A);
} else if (xOptimizationParameter.label == "I") {
global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2;
global.updateReducedPatternFunctions_material.push_back(
function_updateReducedPattern_material_I);
} else if (xOptimizationParameter.label == "I2") {
global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2;
global.updateReducedPatternFunctions_material.push_back(
function_updateReducedPattern_material_I2);
} else if (xOptimizationParameter.label == "I3") {
global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I3;
global.updateReducedPatternFunctions_material.push_back(
function_updateReducedPattern_material_I3);
} else if (xOptimizationParameter.label == "J") {
global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.J;
global.updateReducedPatternFunctions_material.push_back(
function_updateReducedPattern_material_J);
}
//NOTE:Assuming that I2,I3 and J are be passed after E and A because the update of E and A changes I2,I3,J.
// case "HexSize":
// updateReducedPatternFunctions_material.push_back(
// function_updateReducedPattern_material_E);
// break;
// case "HexAngle":
// updateReducedPatternFunctions_material.push_back(
// function_updateReducedPattern_material_E);
// break;
}
if (global.updateReducedPatternFunctions_material.size() != 0) {
function_updateReducedPattern_material =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
for (int optimizationParameterIndex = 0;
optimizationParameterIndex
< global.updateReducedPatternFunctions_material.size();
optimizationParameterIndex++) {
global.updateReducedPatternFunctions_material[optimizationParameterIndex](
pReducedPatternSimulationMesh,
global.initialParameters(optimizationParameterIndex)
* x(optimizationParameterIndex));
}
};
} else {
function_updateReducedPattern_material =
[](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {};
}
}
void ReducedModelOptimizer::computeReducedModelSimulationJob(
const SimulationJob &simulationJobOfFullModel,
const std::unordered_map<size_t, size_t> &fullToReducedMap,
SimulationJob &simulationJobOfReducedModel)
{
assert(simulationJobOfReducedModel.pMesh->VN() != 0);
std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
reducedModelFixedVertices;
for (auto fullModelFixedVertex :
simulationJobOfFullModel.constrainedVertices) {
reducedModelFixedVertices[fullToReducedMap.at(fullModelFixedVertex.first)]
= fullModelFixedVertex.second;
}
std::unordered_map<VertexIndex, Vector6d> reducedModelNodalForces;
for (auto fullModelNodalForce :
simulationJobOfFullModel.nodalExternalForces) {
reducedModelNodalForces[fullToReducedMap.at(fullModelNodalForce.first)]
= fullModelNodalForce.second;
}
std::unordered_map<VertexIndex, Eigen::Vector3d> reducedNodalForcedDisplacements;
for (auto fullForcedDisplacement : simulationJobOfFullModel.nodalForcedDisplacements) {
reducedNodalForcedDisplacements[fullToReducedMap.at(fullForcedDisplacement.first)]
= fullForcedDisplacement.second;
}
simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices;
simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces;
simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel();
simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements;
}
//#if POLYSCOPE_DEFINED
//void ReducedModelOptimizer::visualizeResults(
// const std::vector<std::shared_ptr<SimulationJob>> &fullPatternSimulationJobs,
// const std::vector<std::shared_ptr<SimulationJob>> &reducedPatternSimulationJobs,
// const std::vector<BaseSimulationScenario> &simulationScenarios,
// const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
// &reducedToFullInterfaceViMap)
//{
// DRMSimulationModel simulator;
// std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh = fullPatternSimulationJobs[0]->pMesh;
// pFullPatternSimulationMesh->registerForDrawing();
// pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply");
// reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel()
// + "_undeformed.ply");
// double totalError = 0;
// for (const int simulationScenarioIndex : simulationScenarios) {
// const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
// = fullPatternSimulationJobs[simulationScenarioIndex];
// pFullPatternSimulationJob->registerForDrawing(pFullPatternSimulationMesh->getLabel());
// SimulationResults fullModelResults = simulator.executeSimulation(pFullPatternSimulationJob);
// fullModelResults.registerForDrawing();
// // fullModelResults.saveDeformedModel();
// const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob
// = reducedPatternSimulationJobs[simulationScenarioIndex];
// SimulationResults reducedModelResults = simulator.executeSimulation(
// pReducedPatternSimulationJob);
// double normalizationFactor = 1;
// if (global.optimizationSettings.normalizationStrategy
// != ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
// normalizationFactor
// = global.translationalDisplacementNormalizationValues[simulationScenarioIndex];
// }
// reducedModelResults.saveDeformedModel();
// fullModelResults.saveDeformedModel();
// double error = computeDisplacementError(reducedModelResults.displacements,
// fullModelResults.displacements,
// reducedToFullInterfaceViMap,
// normalizationFactor);
// std::cout << "Error of simulation scenario "
// < baseSimulationScenarioNames[simulationScenarioIndex] << " is " << error << std::endl;
// totalError += error;
// reducedModelResults.registerForDrawing();
// // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing();
// // registerWorldAxes();
// const std::string screenshotFilename
// = "/home/iason/Coding/Projects/Approximating shapes with flat "
// "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
// "Images/"
// + pFullPatternSimulationMesh->getLabel() + "_"
// + baseSimulationScenarioNames[simulationScenarioIndex];
// polyscope::show();
// polyscope::screenshot(screenshotFilename, false);
// fullModelResults.unregister();
// reducedModelResults.unregister();
// // firstOptimizationRoundResults[simulationScenarioIndex].unregister();
// }
// std::cout << "Total error:" << totalError << std::endl;
//}
//void ReducedModelOptimizer::registerResultsForDrawing(
// const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob,
// const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob,
// const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
// &reducedToFullInterfaceViMap) {
// DRMSimulationModel simulator;
// std::shared_ptr<SimulationMesh> pFullPatternSimulationMesh =
// pFullPatternSimulationJob->pMesh;
// pFullPatternSimulationMesh->registerForDrawing();
// // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel()
// // +
// // "_undeformed.ply");
// // reducedPatternSimulationJobs[0]->pMesh->savePly(
// // reducedPatternSimulationJobs[0]->pMesh->getLabel() +
// // "_undeformed.ply");
// pFullPatternSimulationJob->registerForDrawing(
// pFullPatternSimulationMesh->getLabel());
// SimulationResults fullModelResults =
// simulator.executeSimulation(pFullPatternSimulationJob);
// fullModelResults.registerForDrawing();
// // fullModelResults.saveDeformedModel();
// SimulationResults reducedModelResults =
// simulator.executeSimulation(pReducedPatternSimulationJob);
// // reducedModelResults.saveDeformedModel();
// // fullModelResults.saveDeformedModel();
// double error = computeRawDisplacementError(
// reducedModelResults.displacements, fullModelResults.displacements,
// reducedToFullInterfaceViMap/*,
// global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/);
// std::cout << "Error is " << error << std::endl;
// reducedModelResults.registerForDrawing();
//}
//#endif // POLYSCOPE_DEFINED
void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {
assert(optimalDisplacementsOfReducedModel.rows() != 0 &&
optimalDisplacementsOfReducedModel.cols() == 3);
optimalDisplacementsOfReducedModel.setZero(
optimalDisplacementsOfReducedModel.rows(),
optimalDisplacementsOfReducedModel.cols());
for (auto reducedFullViPair : displacementsReducedToFullMap) {
const VertexIndex fullModelVi = reducedFullViPair.second;
const Vector6d fullModelViDisplacements =
fullModelResults.displacements[fullModelVi];
optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) =
Eigen::Vector3d(fullModelViDisplacements[0],
fullModelViDisplacements[1],
fullModelViDisplacements[2]);
}
}
void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimizationResult_dlib,
const Settings &settings,
ReducedPatternOptimization::Results &results)
{
//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());
std::vector<double> optimalX(settings.xRanges.size());
for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) {
if (xVariableIndex < settings.xRanges.size() - 2) {
results.optimalXNameValuePairs[xVariableIndex]
= std::make_pair(settings.xRanges[xVariableIndex].label,
optimizationResult_dlib.x(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,
optimizationResult_dlib.x(xVariableIndex));
}
optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex);
#ifdef POLYSCOPE_DEFINED
std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":"
<< optimalX[xVariableIndex] << " ";
#endif
}
#ifdef POLYSCOPE_DEFINED
std::cout << std::endl;
#endif
// Compute obj value per simulation scenario and the raw objective value
// updateMeshFunction(optimalX);
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh
= global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh;
function_updateReducedPattern(optimizationResult_dlib.x, pReducedPatternSimulationMesh);
// results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios);
//TODO:use push_back it will make the code more readable
LinearSimulationModel simulator;
results.objectiveValue.totalRaw = 0;
results.objectiveValue.perSimulationScenario_translational.resize(
global.simulationScenarioIndices.size());
results.objectiveValue.perSimulationScenario_rawTranslational.resize(
global.simulationScenarioIndices.size());
results.objectiveValue.perSimulationScenario_rotational.resize(
global.simulationScenarioIndices.size());
results.objectiveValue.perSimulationScenario_rawRotational.resize(
global.simulationScenarioIndices.size());
results.objectiveValue.perSimulationScenario_total.resize(
global.simulationScenarioIndices.size());
//#ifdef POLYSCOPE_DEFINED
// global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
//#endif
results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size());
for (int i = 0; i < global.simulationScenarioIndices.size(); i++) {
const int simulationScenarioIndex = global.simulationScenarioIndices[i];
SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
results.objectiveValue.perSimulationScenario_total[i] = computeError(
global.fullPatternResults[simulationScenarioIndex],
reducedModelResults,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]);
//Raw translational
const double rawTranslationalError = computeRawTranslationalError(
global.fullPatternResults[simulationScenarioIndex].displacements,
reducedModelResults.displacements,
global.reducedToFullInterfaceViMap);
results.objectiveValue.perSimulationScenario_rawTranslational[i] = rawTranslationalError;
//Raw rotational
const double rawRotationalError = computeRawRotationalError(
global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion,
reducedModelResults.rotationalDisplacementQuaternion,
global.reducedToFullInterfaceViMap);
results.objectiveValue.perSimulationScenario_rawRotational[i] = rawRotationalError;
//Normalized translational
const double normalizedTranslationalError = computeDisplacementError(
global.fullPatternResults[simulationScenarioIndex].displacements,
reducedModelResults.displacements,
global.reducedToFullInterfaceViMap,
global.translationalDisplacementNormalizationValues[simulationScenarioIndex]);
results.objectiveValue.perSimulationScenario_translational[i] = 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]);
results.objectiveValue.perSimulationScenario_rotational[i] = 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);
results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError;
results.perScenario_fullPatternPotentialEnergy[i]
= global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy;
#ifdef POLYSCOPE_DEFINED
std::cout << "Simulation scenario:"
<< global.reducedPatternSimulationJobs[simulationScenarioIndex]->getLabel()
<< std::endl;
std::cout << "raw translational error:" << rawTranslationalError << std::endl;
std::cout << "translation normalization value:"
<< global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
<< std::endl;
std::cout << "raw Rotational error:" << rawRotationalError << std::endl;
std::cout << "rotational normalization value:"
<< global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
<< std::endl;
std::cout << "Translational error:" << normalizedTranslationalError << std::endl;
std::cout << "Rotational error:" << normalizedRotationalError << std::endl;
// results.objectiveValuePerSimulationScenario[simulationScenarioIndex]
// = normalizedTranslationalError + normalizedRotationalError;
std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i]
<< std::endl;
std::cout << std::endl;
// reducedModelResults.registerForDrawing(Colors::reducedDeformed);
// global.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed);
// polyscope::show();
// reducedModelResults.unregister();
// global.fullPatternResults[simulationScenarioIndex].unregister();
#endif
}
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);
}
results.iterationOfMinima = global.iterationOfMinima;
// results.draw();
}
std::vector<std::pair<BaseSimulationScenario, double>>
ReducedModelOptimizer::computeFullPatternMaxSimulationForces(
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenario)
{
std::vector<std::pair<BaseSimulationScenario, double>> fullPatternMaxSimulationForces;
fullPatternMaxSimulationForces.reserve(desiredBaseSimulationScenario.size());
for (const BaseSimulationScenario &scenario : desiredBaseSimulationScenario) {
const double maxForce = computeFullPatternMaxSimulationForce(scenario);
fullPatternMaxSimulationForces.emplace_back(std::make_pair(scenario, maxForce));
}
return fullPatternMaxSimulationForces;
}
std::vector<std::pair<BaseSimulationScenario, double>>
ReducedModelOptimizer::getFullPatternMaxSimulationForces(
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
{
std::vector<std::pair<BaseSimulationScenario, double>> fullPatternSimulationScenarioMaxMagnitudes;
#ifdef POLYSCOPE_DEFINED
const std::filesystem::path forceMagnitudesDirectoryPath(std::filesystem::current_path()
.parent_path()
.append("IntermediateResults")
.append("ForceMagnitudes"));
std::filesystem::path patternMaxForceMagnitudesFilePath(
std::filesystem::path(forceMagnitudesDirectoryPath)
.append(m_pFullPatternSimulationMesh->getLabel() + ".json"));
const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists(
patternMaxForceMagnitudesFilePath);
if (fullPatternScenarioMagnitudesExist) {
nlohmann::json json;
std::ifstream ifs(patternMaxForceMagnitudesFilePath.string());
ifs >> json;
fullPatternSimulationScenarioMaxMagnitudes
= static_cast<std::vector<std::pair<BaseSimulationScenario, double>>>(
json.at("maxMagn"));
const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size()
!= desiredBaseSimulationScenarioIndices.size();
if (!shouldRecompute) {
return fullPatternSimulationScenarioMaxMagnitudes;
}
}
#endif
fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces(
desiredBaseSimulationScenarioIndices);
#ifdef POLYSCOPE_DEFINED
nlohmann::json json;
json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes;
std::filesystem::create_directories(forceMagnitudesDirectoryPath);
std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string());
jsonFile << json;
#endif
assert(fullPatternSimulationScenarioMaxMagnitudes.size()
== desiredBaseSimulationScenarioIndices.size());
return fullPatternSimulationScenarioMaxMagnitudes;
}
void ReducedModelOptimizer::runOptimization(const Settings &settings,
ReducedPatternOptimization::Results &results)
{
#if POLYSCOPE_DEFINED
global.objectiveValueHistoryY.clear();
global.objectiveValueHistoryY.reserve(settings.numberOfFunctionCalls);
global.objectiveValueHistoryX.reserve(settings.numberOfFunctionCalls);
global.plotColors.reserve(settings.numberOfFunctionCalls);
#endif
double (*objF)(const dlib::matrix<double, 0, 1> &) = &objective;
dlib::function_evaluation optimalResult;
constexpr bool shouldJointlyOptimizeGeometryAndMaterial = false;
const auto hexAngleParameterIt = std::find_if(settings.xRanges.begin(),
settings.xRanges.end(),
[](const xRange &x) {
return x.label == "HexAngle";
});
const bool hasHexAngleParameter = hexAngleParameterIt != settings.xRanges.end();
const auto hexSizeParameterIt = std::find_if(settings.xRanges.begin(),
settings.xRanges.end(),
[](const xRange &x) {
return x.label == "HexSize";
});
const bool hasHexSizeParameter = hexSizeParameterIt != settings.xRanges.end();
const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter;
assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter));
const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0;
const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters
- numberOfGeometryOptimizationParameters;
const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters
&& numberOfMaterialOptimizationParameters != 0;
if (!shouldJointlyOptimizeGeometryAndMaterial && hasGeometryAndMaterialParameters) {
//Geometry optimization of the reduced pattern
dlib::matrix<double, 0, 1> xGeometryMin(numberOfGeometryOptimizationParameters);
dlib::matrix<double, 0, 1> xGeometryMax(numberOfGeometryOptimizationParameters);
// const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(),
// hexAngleParameterIt);
// const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(),
// hexSizeParameterIt);
xGeometryMin(0) = hexSizeParameterIt->min;
xGeometryMax(0) = hexSizeParameterIt->max;
xGeometryMin(1) = hexAngleParameterIt->min;
xGeometryMax(1) = hexAngleParameterIt->max;
//Set reduced pattern update functions to be used during optimization
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
const int optimizationFunctionCalls_geometry
= global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters
? settings.numberOfFunctionCalls
: static_cast<int>(settings.numberOfFunctionCalls / 3.0);
dlib::function_evaluation result_dlib_geometry
= dlib::find_min_global(objF,
xGeometryMin,
xGeometryMax,
dlib::max_function_calls(optimizationFunctionCalls_geometry),
std::chrono::hours(24 * 365 * 290),
settings.solverAccuracy);
//Material optimization of the reduced pattern
// std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl;
// std::cout << "opt size:" << global.minX[0] << std::endl;
dlib::function_evaluation result_dlib_material;
if (numberOfMaterialOptimizationParameters != 0) {
dlib::matrix<double, 0, 1> xMaterialMin(numberOfMaterialOptimizationParameters);
dlib::matrix<double, 0, 1> xMaterialMax(numberOfMaterialOptimizationParameters);
for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) {
xMaterialMin(i) = settings.xRanges[i].min;
xMaterialMax(i) = settings.xRanges[i].max;
}
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
function_updateReducedPattern_geometry(
result_dlib_geometry.x,
pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this
pReducedPatternSimulationMesh->reset();
};
result_dlib_material = dlib::find_min_global(objF,
xMaterialMin,
xMaterialMax,
dlib::max_function_calls(static_cast<int>(
2.0 * settings.numberOfFunctionCalls
/ 3.0)),
std::chrono::hours(24 * 365 * 290),
settings.solverAccuracy);
}
// constexpr bool useBOBYQA = false;
// if (useBOBYQA) {
// const size_t npt = 2 * global.numberOfOptimizationParameters;
// // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
// const double rhobeg = 0.1;
// // const double rhobeg = 10;
// const double rhoend = rhobeg * 1e-6;
// // const size_t maxFun = 10 * (x.size() ^ 2);
// const size_t bobyqa_maxFunctionCalls = 200;
// dlib::find_min_bobyqa(objF,
// result_dlib.x,
// npt,
// xMaterialMin,
// xMaterialMax,
// rhobeg,
// rhoend,
// bobyqa_maxFunctionCalls);
// }
optimalResult.x.set_size(global.numberOfOptimizationParameters);
std::copy(result_dlib_material.x.begin(),
result_dlib_material.x.end(),
optimalResult.x.begin());
// debug_x.begin());
std::copy(result_dlib_geometry.x.begin(),
result_dlib_geometry.x.end(),
optimalResult.x.begin() + numberOfMaterialOptimizationParameters);
std::cout << "opt x:";
for (const auto optx : optimalResult.x) {
std::cout << optx << " ";
}
std::cout << std::endl;
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
// const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel();
// global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre");
// global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial);
// polyscope::show();
// global.reducedPatternSimulationJobs[0]->pMesh->unregister();
optimalResult.y = objective(optimalResult.x);
} else {
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;
}
if (numberOfGeometryOptimizationParameters != 0
&& numberOfMaterialOptimizationParameters != 0) {
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
} else if (numberOfGeometryOptimizationParameters == 0) {
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_material(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
} else {
function_updateReducedPattern =
[&](const dlib::matrix<double, 0, 1> &x,
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh);
pReducedPatternSimulationMesh->reset();
};
}
optimalResult = dlib::find_min_global(objF,
xMin,
xMax,
dlib::max_function_calls(
settings.numberOfFunctionCalls),
std::chrono::hours(24 * 365 * 290),
settings.solverAccuracy);
}
getResults(optimalResult, settings, results);
}
void ReducedModelOptimizer::constructAxialSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job)
{
for (auto viPairIt = oppositeInterfaceViPairs.begin();
viPairIt != oppositeInterfaceViPairs.end();
viPairIt++) {
if (viPairIt != oppositeInterfaceViPairs.begin()) {
CoordType forceDirection(1, 0, 0);
job.nodalExternalForces[viPairIt->first]
= Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
* forceMagnitude;
job.constrainedVertices[viPairIt->second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
}
void ReducedModelOptimizer::constructShearSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job)
{
for (auto viPairIt = oppositeInterfaceViPairs.begin();
viPairIt != oppositeInterfaceViPairs.end();
viPairIt++) {
if (viPairIt != oppositeInterfaceViPairs.begin()) {
CoordType forceDirection(0, 1, 0);
const auto viPair = *viPairIt;
job.nodalExternalForces[viPair.first]
= Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0})
* forceMagnitude;
job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
}
void ReducedModelOptimizer::constructBendingSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job)
{
for (const auto &viPair : oppositeInterfaceViPairs) {
job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude;
job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
}
/*NOTE: From the results it seems as if the forced displacements are different in the linear and in the drm
* */
void ReducedModelOptimizer::constructDomeSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job)
{
for (auto viPairIt = oppositeInterfaceViPairs.begin();
viPairIt != oppositeInterfaceViPairs.end();
viPairIt++) {
const auto viPair = *viPairIt;
CoordType interfaceVector = (job.pMesh->vert[viPair.first].cP()
- 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.005
* 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.005 * 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};
}
}
}
void ReducedModelOptimizer::constructSaddleSimulationScenario(
const double &forceMagnitude,
const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
&oppositeInterfaceViPairs,
SimulationJob &job)
{
for (auto viPairIt = oppositeInterfaceViPairs.begin();
viPairIt != oppositeInterfaceViPairs.end();
viPairIt++) {
const auto &viPair = *viPairIt;
CoordType v = (job.pMesh->vert[viPair.first].cP() - job.pMesh->vert[viPair.second].cP())
^ CoordType(0, 0, -1).Normalize();
if (viPairIt == oppositeInterfaceViPairs.begin()) {
job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0})
* forceMagnitude;
job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0})
* forceMagnitude;
} else {
job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0})
* forceMagnitude / 2;
job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0})
* forceMagnitude / 2;
}
}
}
double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude)
{
SimulationJob job;
job.pMesh = global.pFullPatternSimulationMesh;
global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
// ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
// global.fullPatternInterfaceViPairs,
// job);
DRMSimulationModel simulator;
DRMSimulationModel::Settings settings;
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-8;
settings.viscousDampingFactor = 1e-2;
settings.useKineticDamping = true;
settings.shouldDraw = false;
settings.debugModeStep = 200;
// settings.averageResidualForcesCriterionThreshold = 1e-5;
settings.maxDRMIterations = 200000;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
settings);
const double &desiredRotationAngle = global.desiredMaxRotationAngle;
double error;
#ifdef POLYSCOPE_DEFINED
job.pMesh->setLabel("initial");
// job.pMesh->registerForDrawing();
// results.registerForDrawing();
// polyscope::show();
std::string saveJobToPath;
if (!results.converged) {
// std::cout << "Force used:" << forceMagnitude << std::endl;
error = std::numeric_limits<double>::max();
// DRMSimulationModel::Settings debugSimulationSettings;
// debugSimulationSettings.isDebugMode = true;
// debugSimulationSettings.debugModeStep = 1000;
// // debugSimulationSettings.maxDRMIterations = 100000;
// debugSimulationSettings.shouldDraw = true;
// debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
// debugSimulationSettings.shouldCreatePlots = true;
// // debugSimulationSettings.Dtini = 0.06;
// debugSimulationSettings.beVerbose = true;
// debugSimulationSettings.useAverage = true;
// // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
// debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
// auto debugResults = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
// debugSimulationSettings);
// std::terminate();
saveJobToPath = "../nonConvergingJobs";
} else {
error = std::abs(
results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError]
.angularDistance(Eigen::Quaterniond::Identity())
- desiredRotationAngle);
saveJobToPath = "../convergingJobs";
}
std::filesystem::path outputPath(std::filesystem::path(saveJobToPath)
.append(job.pMesh->getLabel())
.append("mag_" + global.currentScenarioName));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
settings.save(outputPath);
std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl;
#endif
return error;
}
double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude)
{
SimulationJob job;
job.pMesh = global.pFullPatternSimulationMesh;
global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
// ReducedModelOptimizer::axialSimulationScenario(forceMagnitude,
// global.fullPatternInterfaceViPairs,
// job);
DRMSimulationModel simulator;
DRMSimulationModel::Settings settings;
// settings.totalResidualForcesNormThreshold = 1e-3;
settings.totalExternalForcesNormPercentageTermination = 1e-2;
settings.totalTranslationalKineticEnergyThreshold = 1e-8;
settings.viscousDampingFactor = 5e-3;
settings.useKineticDamping = true;
// settings.averageResidualForcesCriterionThreshold = 1e-5;
// settings.useAverage = true;
// settings.totalTranslationalKineticEnergyThreshold = 1e-10;
// settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.shouldDraw = true;
// settings.isDebugMode = true;
// settings.drawingStep = 200000;
// settings.beVerbose = true;
// settings.debugModeStep = 200;
settings.maxDRMIterations = 200000;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
settings);
const double &desiredDisplacementValue = global.desiredMaxDisplacementValue;
double error;
if (!results.converged) {
error = std::numeric_limits<double>::max();
std::filesystem::path outputPath(std::filesystem::path("../nonConvergingJobs")
.append(job.pMesh->getLabel())
.append("mag_" + global.currentScenarioName));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
settings.save(outputPath);
// std::terminate();
} else {
error = std::abs(
results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm()
- desiredDisplacementValue);
}
#ifdef POLYSCOPE_DEFINED
std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl;
#endif
return error;
}
double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
const BaseSimulationScenario &scenario)
{
double forceMagnitude = 1;
double minimumError;
bool wasSuccessful = false;
global.constructScenarioFunction = constructBaseScenarioFunctions[scenario];
const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
std::function<double(const double &)> objectiveFunction;
double translationalOptimizationEpsilon{baseTriangleHeight * 0.001};
double objectiveEpsilon = translationalOptimizationEpsilon;
objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
global.desiredMaxDisplacementValue = 0.1 * baseTriangleHeight;
global.currentScenarioName = baseSimulationScenarioNames[scenario];
double forceMagnitudeEpsilon = 1e-4;
switch (scenario) {
case Axial:
global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight;
break;
case Shear:
global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight;
break;
case Bending:
break;
case Dome:
global.desiredMaxRotationAngle = vcg::math::ToRad(25.0);
objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective;
forceMagnitudeEpsilon *= 1e-2;
objectiveEpsilon = vcg::math::ToRad(3.0);
forceMagnitude = 0.6;
break;
case Saddle:
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
break;
}
minimumError = dlib::find_min_single_variable(objectiveFunction,
forceMagnitude,
1e-4,
1e4,
forceMagnitudeEpsilon,
1000);
wasSuccessful = minimumError < objectiveEpsilon;
#ifdef POLYSCOPE_DEFINED
std::cout << "Max " << ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
<< " magnitude:" << forceMagnitude << std::endl;
if (!wasSuccessful) {
SimulationJob job;
job.pMesh = global.pFullPatternSimulationMesh;
global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job);
std::cout << ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
+ " max scenario magnitude was not succefully determined."
<< std::endl;
std::filesystem::path outputPath(
std::filesystem::path("../nonConvergingJobs")
.append(m_pFullPatternSimulationMesh->getLabel())
.append("magFinal_"
+ ReducedPatternOptimization::baseSimulationScenarioNames[scenario]));
std::filesystem::create_directories(outputPath);
job.save(outputPath);
std::terminate();
}
#endif
return forceMagnitude;
}
std::vector<std::shared_ptr<SimulationJob>> ReducedModelOptimizer::createFullPatternSimulationJobs(
const std::shared_ptr<SimulationMesh> &pMesh,
const std::vector<std::pair<BaseSimulationScenario, double>> &scenarioMaxForceMagnitudePairs)
{
std::vector<std::shared_ptr<SimulationJob>> scenarios;
scenarios.resize(totalNumberOfSimulationScenarios);
SimulationJob job;
job.pMesh = pMesh;
for (std::pair<BaseSimulationScenario, double> scenarioMaxForceMagnitudePair :
scenarioMaxForceMagnitudePairs) {
const BaseSimulationScenario scenario = scenarioMaxForceMagnitudePair.first;
const double maxForceMagnitude = scenarioMaxForceMagnitudePair.second;
const int numberOfSimulationScenarios = simulationScenariosResolution[scenario];
const double minForceMagnitude = scenarioIsSymmetrical[scenario]
? maxForceMagnitude / numberOfSimulationScenarios
: -maxForceMagnitude;
const double forceMagnitudeStep = numberOfSimulationScenarios == 1
? maxForceMagnitude
: (maxForceMagnitude - minForceMagnitude)
/ (numberOfSimulationScenarios);
const int baseSimulationScenarioIndexOffset
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + scenario,
0);
for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationScenarios;
simulationScenarioIndex++) {
job.nodalExternalForces.clear();
job.constrainedVertices.clear();
job.nodalForcedDisplacements.clear();
job.label = baseSimulationScenarioNames[scenario] + "_"
+ std::to_string(simulationScenarioIndex);
const double forceMagnitude = (forceMagnitudeStep * simulationScenarioIndex
+ minForceMagnitude);
constructBaseScenarioFunctions[scenario](forceMagnitude,
m_fullPatternOppositeInterfaceViPairs,
job);
scenarios[baseSimulationScenarioIndexOffset + simulationScenarioIndex]
= std::make_shared<SimulationJob>(job);
}
}
#ifdef POLYSCOPE_DEFINED
std::cout << "Computed full pattern scenario magnitudes" << std::endl;
#endif
return scenarios;
}
void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
// Compute the sum of the displacement norms
std::vector<double> fullPatternTranslationalDisplacementNormSum(
totalNumberOfSimulationScenarios);
std::vector<double> fullPatternAngularDistance(totalNumberOfSimulationScenarios);
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
double translationalDisplacementNormSum = 0;
for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
const int fullPatternVi = interfaceViPair.second;
//If the full pattern vertex is translationally constrained dont take it into account
if (global.fullPatternSimulationJobs[simulationScenarioIndex]
->constrainedVertices.contains(fullPatternVi)) {
const std::unordered_set<int> constrainedDof
= global.fullPatternSimulationJobs[simulationScenarioIndex]
->constrainedVertices.at(fullPatternVi);
if (constrainedDof.contains(0) && constrainedDof.contains(1)
&& constrainedDof.contains(2)) {
continue;
}
}
const Vector6d &vertexDisplacement = global
.fullPatternResults[simulationScenarioIndex]
.displacements[fullPatternVi];
translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm();
}
double angularDistanceSum = 0;
for (auto interfaceViPair : global.reducedToFullInterfaceViMap) {
const int fullPatternVi = interfaceViPair.second;
//If the full pattern vertex is rotationally constrained dont take it into account
if (global.fullPatternSimulationJobs[simulationScenarioIndex]
->constrainedVertices.contains(fullPatternVi)) {
const std::unordered_set<int> constrainedDof
= global.fullPatternSimulationJobs[simulationScenarioIndex]
->constrainedVertices.at(fullPatternVi);
if (constrainedDof.contains(3) && constrainedDof.contains(5)
&& constrainedDof.contains(4)) {
continue;
}
}
angularDistanceSum += global.fullPatternResults[simulationScenarioIndex]
.rotationalDisplacementQuaternion[fullPatternVi]
.angularDistance(Eigen::Quaterniond::Identity());
}
fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex]
= translationalDisplacementNormSum;
fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum;
}
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
if (global.optimizationSettings.normalizationStrategy ==
Settings::NormalizationStrategy::Epsilon) {
const double epsilon_translationalDisplacement
= global.optimizationSettings.translationNormalizationParameter;
global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
epsilon_translationalDisplacement);
const double epsilon_rotationalDisplacement = global.optimizationSettings
.rotationNormalizationParameter;
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternAngularDistance[simulationScenarioIndex],
epsilon_rotationalDisplacement);
} else {
global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1;
}
}
}
void ReducedModelOptimizer::optimize(
const Settings &optimizationSettings,
ReducedPatternOptimization::Results &results,
const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
{
for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) {
//Increase the size of the vector holding the simulation scenario indices
global.simulationScenarioIndices.resize(
global.simulationScenarioIndices.size()
+ simulationScenariosResolution[baseSimulationScenarioIndex]);
//Add the simulation scenarios indices that correspond to this base simulation scenario
std::iota(global.simulationScenarioIndices.end()
- simulationScenariosResolution[baseSimulationScenarioIndex],
global.simulationScenarioIndices.end(),
std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin()
+ baseSimulationScenarioIndex,
0));
}
global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios);
global.fullPatternResults.resize(totalNumberOfSimulationScenarios);
global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
global.rotationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
global.minY = std::numeric_limits<double>::max();
global.numOfSimulationCrashes = 0;
global.numberOfFunctionCalls = 0;
global.optimizationSettings = optimizationSettings;
global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh;
std::vector<std::pair<BaseSimulationScenario, double>> fullPatternSimulationScenarioMaxMagnitudes
= getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices);
global.fullPatternSimulationJobs
= createFullPatternSimulationJobs(m_pFullPatternSimulationMesh,
fullPatternSimulationScenarioMaxMagnitudes);
// polyscope::removeAllStructures();
results.baseTriangle = global.baseTriangle;
DRMSimulationModel::Settings simulationSettings;
// simulationSettings.maxDRMIterations = 200000;
// simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
// simulationSettings.viscousDampingFactor = 5e-3;
// simulationSettings.useKineticDamping = true;
// 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) {
global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
ReducedPatternOptimization::Colors::fullInitial);
}
#endif
results.wasSuccessful = true;
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings);
// 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
if (!fullPatternResults.converged) {
results.wasSuccessful = false;
#ifdef POLYSCOPE_DEFINED
std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
<< " did not converge." << std::endl;
// 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);
#endif
std::terminate();
return;
// }
}
#ifdef POLYSCOPE_DEFINED
if (drawFullPatternSimulationResults) {
// SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation(
// pFullPatternSimulationJob);
fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed,
true);
// fullPatternResults_linear.labelPrefix += "_linear";
// fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed,
// true,
// true);
polyscope::show();
fullPatternResults.unregister();
// fullPatternResults_linear.unregister();
}
#endif
global.fullPatternResults[simulationScenarioIndex] = fullPatternResults;
SimulationJob reducedPatternSimulationJob;
reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh;
computeReducedModelSimulationJob(*pFullPatternSimulationJob,
m_fullToReducedInterfaceViMap,
reducedPatternSimulationJob);
global.reducedPatternSimulationJobs[simulationScenarioIndex]
= std::make_shared<SimulationJob>(reducedPatternSimulationJob);
// std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl;
}
#ifdef POLYSCOPE_DEFINED
if (drawFullPatternSimulationResults) {
global.fullPatternSimulationJobs[0]->pMesh->unregister();
}
#endif
// if (global.optimizationSettings.normalizationStrategy
// != Settings::NormalizationStrategy::NonNormalized) {
computeObjectiveValueNormalizationFactors();
// }
#ifdef POLYSCOPE_DEFINED
std::cout << "Running reduced model optimization" << std::endl;
#endif
runOptimization(optimizationSettings, results);
results.notes = optimizationNotes;
}
Reference in New Issue View Git Blame Copy Permalink