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

1540 lines
79 KiB
C++
Raw Blame History

#include "reducedmodeloptimizer.hpp"
#include "linearsimulationmodel.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepatterngeometry.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
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> objectiveValueHistory;
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;
} 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);
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);
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(long n, const double *x) {
// std::cout.precision(17);
// for (int i = 0; i < n; 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;
updateMesh(n, x);
// global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing();
// global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing();
// polyscope::show();
// global.reducedPatternSimulationJobs[0]->pMesh->unregister();
// std::cout << e.axialConstFactor << " " << 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];
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 {
const 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
// std::cout << "sim error:" << simulationScenarioError << std::endl;
// global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(
// ReducedModelOptimization::Colors::reducedInitial);
// reducedModelResults.registerForDrawing(
// ReducedModelOptimization::Colors::reducedDeformed);
// polyscope::show();
// reducedModelResults.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;
}
}
// std::cout << error << std::endl;
if (totalError < global.minY) {
global.minY = totalError;
global.minX.assign(x, x + n);
}
#ifdef POLYSCOPE_DEFINED
++global.numberOfFunctionCalls;
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
// global.objectiveValueHistory.push_back(totalError);
// auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(),
// global.objectiveValueHistory.size());
// std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
// if (global.g_firstRoundIterationIndex != 0) {
// for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
// [](double &c) { c = 0.7; });
// }
// global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
// SimulationResultsReporter::createPlot("Number of Steps", "Objective value",
// global.objectiveValueHistory);
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();
}
// Full pattern
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
pFullPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pFullPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
// Reduced pattern
pReducedPatternSimulationMesh =
std::make_shared<SimulationMesh>(reducedModel);
pReducedPatternSimulationMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pReducedPatternSimulationMesh->setBeamMaterial(0.3, 1 * 1e9);
}
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);
}
void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
PatternGeometry &reducedPattern,
const int &optimizationParameters)
{
assert(fullPattern.VN() == reducedPattern.VN() &&
fullPattern.EN() >= reducedPattern.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);
initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters);
}
void updateMesh(long n, const double *x) {
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
->pMesh;
const double E=global.initialParameters(0)*x[0];
const double A=global.initialParameters(1) * x[1];
const double beamWidth=std::sqrt(A);
const double beamHeight=beamWidth;
const double J=global.initialParameters(2) * x[2];
const double I2=global.initialParameters(3) * x[3];
const double I3=global.initialParameters(4) * x[4];
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(
RectangularBeamDimensions(beamWidth,
beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio,
E));
e.J = J;
e.I2 = I2;
e.I3 = I3;
}
assert(pReducedPatternSimulationMesh->EN() == 12);
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;
}
pReducedPatternSimulationMesh->reset();
//#ifdef POLYSCOPE_DEFINED
// pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
// pReducedPatternSimulationMesh->registerForDrawing();
// std::cout << "Angle:" + std::to_string(x[n - 1]) + " size:" + std::to_string(x[n - 2])
// << std::endl;
// std::cout << "Verts:" << pReducedPatternSimulationMesh->VN() << std::endl;
// polyscope::show();
//#endif
}
void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh,const int& optimizationParamters) {
global.numberOfOptimizationParameters = optimizationParamters;
global.initialParameters.resize(global.numberOfOptimizationParameters);
global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
global.initialParameters(1) = mesh->elements[0].A;
global.initialParameters(2) = mesh->elements[0].J;
global.initialParameters(3) = mesh->elements[0].I2;
global.initialParameters(4) = mesh->elements[0].I3;
}
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 < 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);
#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
updateMesh(optimalX.size(), optimalX.data());
// 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());
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;
#endif
}
const bool printDebugInfo = false;
if (printDebugInfo) {
std::cout << "Finished optimizing." << endl;
std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl;
assert(global.minY == optimizationResult_dlib.y);
if (global.minY != optimizationResult_dlib.y) {
std::cerr << "ERROR in objective value" << std::endl;
}
}
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);
}
}
void ReducedModelOptimizer::runOptimization(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;
}
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.solverAccuracy);
getResults(result_dlib, 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.01
* std::abs(
forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude?
// * std::abs(forceMagnitude / maxForceMagnitude_dome);
job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0],
interfaceVector[1],
0)
* 0.001 * std::abs(forceMagnitude);
// * std::abs(forceMagnitude / maxForceMagnitude_dome);
// CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// ^ CoordType(0, 0, -1).Normalize();
// nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude
// * 0.0001;
// nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude
// * 0.0001;
} else {
job.nodalExternalForces[viPair.first]
= Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude
/ 3;
job.nodalExternalForces[viPair.second]
= Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]})
* forceMagnitude / 3;
job.constrainedVertices[viPair.first] = std::unordered_set<DoFType>{2};
job.constrainedVertices[viPair.second] = std::unordered_set<DoFType>{2};
}
}
}
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.totalResidualForcesNormThreshold = 1e-3;
// settings.totalTranslationalKineticEnergyThreshold = 1e-6;
// settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.shouldDraw = true;
// settings.isDebugMode = true;
// settings.drawingStep = 500;
// settings.beVerbose = true;
// settings.debugModeStep = 200000;
// settings.shouldCreatePlots = true;
settings.maxDRMIterations = 100000;
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job),
settings);
const double &desiredRotationAngle = global.desiredMaxRotationAngle;
if (!results.converged) {
return std::numeric_limits<double>::max();
}
const double error = std::abs(
// results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError]
.angularDistance(Eigen::Quaterniond::Identity())
- desiredRotationAngle);
#ifdef POLYSCOPE_DEFINED
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.totalTranslationalKineticEnergyThreshold = 1e-10;
// settings.shouldUseTranslationalKineticEnergyThreshold = true;
// settings.shouldDraw = true;
// settings.isDebugMode = true;
// settings.drawingStep = 200000;
// settings.beVerbose = true;
// settings.debugModeStep = 200000;
settings.maxDRMIterations = 50000;
#ifdef POLYSCOPE_DEFINED
// settings.shouldDraw = true;
// settings.isDebugMode = true;
// settings.drawingStep = 100000;
// settings.debugModeStep = 10000;
// settings.shouldCreatePlots = true;
#endif
SimulationResults results = simulator.executeSimulation(std::make_shared<SimulationJob>(job));
const double &desiredDisplacementValue = global.desiredMaxDisplacementValue;
if (!results.converged) {
return std::numeric_limits<double>::max();
}
const double error = std::abs(
// results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm()
results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm()
- desiredDisplacementValue);
#ifdef POLYSCOPE_DEFINED
std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl;
#endif
return error;
}
double ReducedModelOptimizer::getFullPatternMaxSimulationForce(
const BaseSimulationScenario &scenario, bool &wasSuccessful)
{
double forceMagnitude = 1;
const double forceMagnitudeEpsilon = 1e-2;
double minimumError;
double translationalOptimizationEpsilon;
dlib::function_evaluation result;
switch (scenario) {
case Axial:
global.desiredMaxDisplacementValue = 0.03
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.constructScenarioFunction = &ReducedModelOptimizer::constructAxialSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
break;
case Shear:
global.desiredMaxDisplacementValue = 0.04
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.constructScenarioFunction = &ReducedModelOptimizer::constructShearSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
break;
case Bending:
global.desiredMaxDisplacementValue = 0.1
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.constructScenarioFunction = &ReducedModelOptimizer::constructBendingSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
break;
case Dome:
global.desiredMaxRotationAngle = vcg::math::ToRad(35.0);
global.constructScenarioFunction = &ReducedModelOptimizer::constructDomeSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceRotationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
wasSuccessful = minimumError < vcg::math::ToRad(3.0);
// result = dlib::find_min_global(&fullPatternMaxSimulationForceRotationalObjective,
// 1e-2,
// 1,
// dlib::max_function_calls(50));
// wasSuccessful = result.y < vcg::math::ToRad(3.0);
break;
case Saddle:
global.desiredMaxDisplacementValue = 0.1
* vcg::Distance(global.baseTriangle.cP(0),
(global.baseTriangle.cP(1)
+ global.baseTriangle.cP(2))
/ 2);
global.constructScenarioFunction = &ReducedModelOptimizer::constructSaddleSimulationScenario;
global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first;
minimumError
= dlib::find_min_single_variable(&fullPatternMaxSimulationForceTranslationalObjective,
forceMagnitude,
1e-2,
1e2,
forceMagnitudeEpsilon);
translationalOptimizationEpsilon = 1e-2 * global.desiredMaxDisplacementValue;
wasSuccessful = minimumError < translationalOptimizationEpsilon;
break;
}
#ifdef POLYSCOPE_DEFINED
std::cout << "Max " << ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
<< " magnitude:" << forceMagnitude << std::endl;
if (!wasSuccessful)
std::cout << "Was not successfull" << std::endl;
#endif
if (!wasSuccessful) {
const std::string debugMessage
= ReducedPatternOptimization::baseSimulationScenarioNames[scenario]
+ " max scenario magnitude was not succefully determined.";
optimizationNotes.append(debugMessage);
}
return forceMagnitude;
}
//TODO: Make more compact
std::vector<std::shared_ptr<SimulationJob>>
ReducedModelOptimizer::createFullPatternSimulationScenarios(
const std::shared_ptr<SimulationMesh> &pMesh)
{
std::vector<std::shared_ptr<SimulationJob>> scenarios;
scenarios.resize(totalNumberOfSimulationScenarios);
SimulationJob job;
job.pMesh = pMesh;
//// Axial
bool wasSuccessful;
const double maxForceMagnitude_axial
= getFullPatternMaxSimulationForce(BaseSimulationScenario::Axial, wasSuccessful);
const double minForceMagnitude_axial = -maxForceMagnitude_axial;
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++) {
job.nodalExternalForces.clear();
job.constrainedVertices.clear();
job.nodalForcedDisplacements.clear();
job.label = baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_"
+ std::to_string(axialSimulationScenarioIndex);
const double forceMagnitude = (forceMagnitudeStep_axial * axialSimulationScenarioIndex
+ minForceMagnitude_axial);
constructAxialSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
scenarios[baseSimulationScenarioIndexOffset_axial + axialSimulationScenarioIndex]
= std::make_shared<SimulationJob>(job);
}
//// Shear
const double maxForceMagnitude_shear
= getFullPatternMaxSimulationForce(BaseSimulationScenario::Shear, wasSuccessful);
const double minForceMagnitude_shear = -maxForceMagnitude_shear;
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++) {
job.constrainedVertices.clear();
job.nodalExternalForces.clear();
job.label = baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_"
+ std::to_string(shearSimulationScenarioIndex);
const double forceMagnitude = (forceMagnitudeStep_shear * shearSimulationScenarioIndex
+ minForceMagnitude_shear);
constructShearSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
scenarios[baseSimulationScenarioIndexOffset_shear + shearSimulationScenarioIndex]
= std::make_shared<SimulationJob>(job);
}
//// Bending
const double maxForceMagnitude_bending
= getFullPatternMaxSimulationForce(BaseSimulationScenario::Bending, wasSuccessful);
const double minForceMagnitude_bending = 0;
const int numberOfSimulationScenarios_bending
= simulationScenariosResolution[BaseSimulationScenario::Bending];
const double forceMagnitudeStep_bending = (maxForceMagnitude_bending - minForceMagnitude_bending)
/ (numberOfSimulationScenarios_bending);
const int baseSimulationScenarioIndexOffset_bending
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Bending,
0);
for (int bendingSimulationScenarioIndex = 0;
bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending;
bendingSimulationScenarioIndex++) {
job.nodalExternalForces.clear();
job.constrainedVertices.clear();
job.label = baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_"
+ std::to_string(bendingSimulationScenarioIndex);
const double forceMagnitude = (forceMagnitudeStep_bending
* (bendingSimulationScenarioIndex + 1)
+ minForceMagnitude_bending);
constructBendingSimulationScenario(forceMagnitude,
m_fullPatternOppositeInterfaceViPairs,
job);
scenarios[baseSimulationScenarioIndexOffset_bending + bendingSimulationScenarioIndex]
= std::make_shared<SimulationJob>(job);
}
//// Dome
const double maxForceMagnitude_dome
= getFullPatternMaxSimulationForce(BaseSimulationScenario::Dome, wasSuccessful);
const double minForceMagnitude_dome = 0;
const int numberOfSimulationScenarios_dome
= simulationScenariosResolution[BaseSimulationScenario::Dome];
const double forceMagnitudeStep_dome = (maxForceMagnitude_dome - minForceMagnitude_dome)
/ (numberOfSimulationScenarios_dome);
const int baseSimulationScenarioIndexOffset_dome
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Dome,
0);
for (int domeSimulationScenarioIndex = 0;
domeSimulationScenarioIndex < numberOfSimulationScenarios_dome;
domeSimulationScenarioIndex++) {
job.constrainedVertices.clear();
job.nodalExternalForces.clear();
job.nodalForcedDisplacements.clear();
job.label = baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_"
+ std::to_string(domeSimulationScenarioIndex);
const double forceMagnitude = (forceMagnitudeStep_dome * (domeSimulationScenarioIndex + 1)
+ minForceMagnitude_dome);
constructDomeSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job);
scenarios[baseSimulationScenarioIndexOffset_dome + domeSimulationScenarioIndex]
= std::make_shared<SimulationJob>(job);
}
//// Saddle
const double maxForceMagnitude_saddle
= getFullPatternMaxSimulationForce(BaseSimulationScenario::Saddle, wasSuccessful);
const double minForceMagnitude_saddle = 0;
const int numberOfSimulationScenarios_saddle
= simulationScenariosResolution[BaseSimulationScenario::Saddle];
const double forceMagnitudeStep_saddle = (maxForceMagnitude_saddle - minForceMagnitude_saddle)
/ numberOfSimulationScenarios_saddle;
const int baseSimulationScenarioIndexOffset_saddle
= std::accumulate(simulationScenariosResolution.begin(),
simulationScenariosResolution.begin() + BaseSimulationScenario::Saddle,
0);
for (int saddleSimulationScenarioIndex = 0;
saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle;
saddleSimulationScenarioIndex++) {
job.constrainedVertices.clear();
job.nodalExternalForces.clear();
job.nodalForcedDisplacements.clear();
job.label = baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_"
+ std::to_string(saddleSimulationScenarioIndex);
const double forceMagnitude = (forceMagnitudeStep_saddle
* (saddleSimulationScenarioIndex + 1)
+ minForceMagnitude_saddle);
constructSaddleSimulationScenario(forceMagnitude,
m_fullPatternOppositeInterfaceViPairs,
job);
scenarios[baseSimulationScenarioIndexOffset_saddle + saddleSimulationScenarioIndex]
= 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
.normalizationParameter;
global.translationalDisplacementNormalizationValues[simulationScenarioIndex]
= std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex],
epsilon_translationalDisplacement);
// const double epsilon_rotationalDisplacement = vcg::math::ToRad(10.0);
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));
}
if (desiredBaseSimulationScenarioIndices.empty()) {
global.simulationScenarioIndices.resize(totalNumberOfSimulationScenarios);
std::iota(global.simulationScenarioIndices.begin(),
global.simulationScenarioIndices.end(),
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;
global.fullPatternSimulationJobs = createFullPatternSimulationScenarios(
m_pFullPatternSimulationMesh);
// polyscope::removeAllStructures();
results.baseTriangle = global.baseTriangle;
DRMSimulationModel::Settings simulationSettings;
simulationSettings.maxDRMIterations = 200000;
// simulationSettings.shouldDraw = true;
// simulationSettings.isDebugMode = true;
// simulationSettings.debugModeStep = 100000;
#ifdef POLYSCOPE_DEFINED
const bool drawFullPatternSimulationResults = false;
if (drawFullPatternSimulationResults) {
global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
ReducedPatternOptimization::Colors::fullInitial);
}
#endif
// LinearSimulationModel linearSimulator;
for (int simulationScenarioIndex : global.simulationScenarioIndices) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob
= global.fullPatternSimulationJobs[simulationScenarioIndex];
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
simulationSettings);
if (!fullPatternResults.converged) {
results.wasSuccessful = false;
#ifdef POLYSCOPE_DEFINED
std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
<< " did not converge." << std::endl;
DRMSimulationModel::Settings debugSimulationSettings;
debugSimulationSettings.isDebugMode = true;
debugSimulationSettings.debugModeStep = 1000;
debugSimulationSettings.shouldCreatePlots = true;
simulator.executeSimulation(pFullPatternSimulationJob, simulationSettings);
#endif
return;
}
results.wasSuccessful = true;
#ifdef POLYSCOPE_DEFINED
if (drawFullPatternSimulationResults) {
// 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);
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