1865 lines
96 KiB
C++
1865 lines
96 KiB
C++
#include "reducedmodeloptimizer.hpp"
|
|
#include "linearsimulationmodel.hpp"
|
|
#include "simulationhistoryplotter.hpp"
|
|
#include "trianglepatterngeometry.hpp"
|
|
#include "trianglepattterntopology.hpp"
|
|
#include <armadillo>
|
|
#include <chrono>
|
|
#include <execution>
|
|
#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<size_t> objectiveValueHistory_iteration;
|
|
std::vector<double> objectiveValueHistory;
|
|
std::vector<double> plotColors;
|
|
std::array<double, 7> parametersInitialValue;
|
|
std::array<double, 7> optimizationInitialValue;
|
|
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::array<std::function<void(const double &newValue,
|
|
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh)>,
|
|
7>
|
|
functions_updateReducedPatternParameter;
|
|
std::vector<double> xMin;
|
|
std::vector<double> xMax;
|
|
std::vector<double> scenarioEqualizationWeight;
|
|
std::vector<double> scenarioUserWeights;
|
|
} 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());
|
|
//}
|
|
|
|
#ifdef USE_ENSMALLEN
|
|
struct EnsmallenOptimizationObjective
|
|
{
|
|
static double Evaluate(const arma::mat &x_arma)
|
|
{
|
|
std::vector<double> x(x_arma.begin(), x_arma.end());
|
|
for (int xi = 0; xi < x.size(); xi++) {
|
|
if (x[xi] > global.xMax[xi]) {
|
|
return std::numeric_limits<double>::max();
|
|
x[xi] = global.xMax[xi];
|
|
// std::cout << "Out of range" << std::endl;
|
|
} else if (x[xi] < global.xMin[xi]) {
|
|
return std::numeric_limits<double>::max();
|
|
x[xi] = global.xMin[xi];
|
|
// std::cout << "Out of range" << std::endl;
|
|
}
|
|
}
|
|
return ReducedModelOptimizer::objective(x);
|
|
}
|
|
};
|
|
#endif
|
|
|
|
#ifdef DLIB_DEFINED
|
|
double ReducedModelOptimizer::objective(const dlib::matrix<double, 0, 1> &x)
|
|
{
|
|
return objective(std::vector<double>(x.begin(), x.end()));
|
|
}
|
|
#endif
|
|
|
|
double ReducedModelOptimizer::objective(const double &xValue)
|
|
{
|
|
return objective({xValue});
|
|
}
|
|
|
|
double ReducedModelOptimizer::objective(const std::vector<double> &x)
|
|
{
|
|
// std::cout.precision(17);
|
|
// for (int i = 0; i < x.size(); i++) {
|
|
// std::cout << x[i] << " ";
|
|
// }
|
|
// std::cout << std::endl;
|
|
|
|
// std::cout << x(x.size() - 2) << " " << x(x.size() - 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;
|
|
simulator.setStructure(pReducedPatternSimulationMesh);
|
|
// simulator.initialize();
|
|
// FormFinder simulator;
|
|
|
|
std::for_each(
|
|
std::execution::par_unseq,
|
|
global.simulationScenarioIndices.begin(),
|
|
global.simulationScenarioIndices.end(),
|
|
[&](const int &simulationScenarioIndex) {
|
|
// 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 {
|
|
// double simulationScenarioError;
|
|
// constexpr bool usePotentialEnergy = false;
|
|
// if (usePotentialEnergy) {
|
|
// simulationScenarioError = std::abs(
|
|
// reducedModelResults.internalPotentialEnergy
|
|
// - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy);
|
|
// } else {
|
|
const double 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 += /*global.scenarioUserWeights[simulationScenarioIndex]
|
|
**/
|
|
/*global.scenarioEqualizationWeight[simulationScenarioIndex]
|
|
**/
|
|
simulationScenarioError * simulationScenarioError;
|
|
}
|
|
});
|
|
// std::cout << totalError << std::endl;
|
|
// global.objectiveValueHistory.push_back(totalError);
|
|
// global.plotColors.push_back(10);
|
|
++global.numberOfFunctionCalls;
|
|
if (totalError < global.minY) {
|
|
global.minY = totalError;
|
|
// global.objectiveValueHistory.push_back(totalError);
|
|
// global.objectiveValueHistory_iteration.push_back(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();
|
|
}
|
|
// Full pattern
|
|
pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
|
|
pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
|
|
pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus);
|
|
|
|
// Reduced pattern
|
|
pReducedPatternSimulationMesh = std::make_shared<SimulationMesh>(reducedModel);
|
|
pReducedPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
|
|
pReducedPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus);
|
|
}
|
|
|
|
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()
|
|
{
|
|
global.functions_updateReducedPatternParameter[R] =
|
|
[](const double &newR, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR,
|
|
newR / 2,
|
|
newR / 2);
|
|
const CoordType hexagonBaseTriangleVertexPosition
|
|
= global.baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0]
|
|
+ global.baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1]
|
|
+ global.baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2];
|
|
|
|
for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
|
|
rotationCounter++) {
|
|
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
|
|
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
|
|
vcg::math::ToRad(60.0 * rotationCounter))
|
|
* hexagonBaseTriangleVertexPosition;
|
|
}
|
|
};
|
|
|
|
global.functions_updateReducedPatternParameter[Theta] =
|
|
[](const double &newTheta, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
const CoordType baseTriangleHexagonVertexPosition
|
|
= pReducedPatternSimulationMesh->vert[0].cP();
|
|
|
|
const CoordType thetaRotatedHexagonBaseTriangleVertexPosition
|
|
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
|
|
vcg::math::ToRad(newTheta))
|
|
* baseTriangleHexagonVertexPosition;
|
|
|
|
for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize;
|
|
rotationCounter++) {
|
|
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P()
|
|
= vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal,
|
|
vcg::math::ToRad(60.0 * rotationCounter))
|
|
* thetaRotatedHexagonBaseTriangleVertexPosition;
|
|
}
|
|
};
|
|
|
|
global.functions_updateReducedPatternParameter[E] =
|
|
[](const double &newE, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
|
|
Element &e = pReducedPatternSimulationMesh->elements[ei];
|
|
e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE));
|
|
}
|
|
};
|
|
global.functions_updateReducedPatternParameter[A] =
|
|
[](const double &newA, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
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));
|
|
}
|
|
};
|
|
|
|
global.functions_updateReducedPatternParameter[I2] =
|
|
[](const double &newI2, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
|
|
Element &e = pReducedPatternSimulationMesh->elements[ei];
|
|
e.inertia.I2 = newI2;
|
|
e.updateRigidity();
|
|
}
|
|
};
|
|
global.functions_updateReducedPatternParameter[I3] =
|
|
[](const double &newI3, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
|
|
Element &e = pReducedPatternSimulationMesh->elements[ei];
|
|
e.inertia.I3 = newI3;
|
|
e.updateRigidity();
|
|
}
|
|
};
|
|
global.functions_updateReducedPatternParameter[J] =
|
|
[](const double &newJ, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
|
|
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
|
|
Element &e = pReducedPatternSimulationMesh->elements[ei];
|
|
e.inertia.J = newJ;
|
|
e.updateRigidity();
|
|
}
|
|
};
|
|
}
|
|
|
|
void ReducedModelOptimizer::initializeOptimizationParameters(
|
|
const std::shared_ptr<SimulationMesh> &mesh, const std::vector<xRange> &optimizationParameters)
|
|
{
|
|
global.numberOfOptimizationParameters = NumberOfOptimizationParameters;
|
|
|
|
for (int optimizationParameterIndex = 0;
|
|
optimizationParameterIndex < optimizationParameters.size();
|
|
optimizationParameterIndex++) {
|
|
for (int optimizationParameterIndex = E;
|
|
optimizationParameterIndex != NumberOfOptimizationParameters;
|
|
optimizationParameterIndex++) {
|
|
// const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex];
|
|
switch (optimizationParameterIndex) {
|
|
case E:
|
|
global.parametersInitialValue[optimizationParameterIndex]
|
|
= mesh->elements[0].material.youngsModulus;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 1;
|
|
break;
|
|
case A:
|
|
global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].A;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 1;
|
|
break;
|
|
case I2:
|
|
global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0]
|
|
.inertia.I2;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 1;
|
|
break;
|
|
case I3:
|
|
global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0]
|
|
.inertia.I3;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 1;
|
|
break;
|
|
case J:
|
|
global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0]
|
|
.inertia.J;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 1;
|
|
break;
|
|
case R:
|
|
global.parametersInitialValue[optimizationParameterIndex] = 0;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 0.5;
|
|
break;
|
|
case Theta:
|
|
global.parametersInitialValue[optimizationParameterIndex] = 0;
|
|
global.optimizationInitialValue[optimizationParameterIndex] = 0;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
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;
|
|
}
|
|
|
|
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 FunctionEvaluation &optimalObjective,
|
|
const Settings &settings,
|
|
ReducedPatternOptimization::Results &results)
|
|
{
|
|
//Number of crashes
|
|
// results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
|
|
//Value of optimal objective Y
|
|
results.objectiveValue.total = optimalObjective.y;
|
|
// results.objectiveValue.total = 0;
|
|
if (optimalObjective.y != global.minY) {
|
|
std::cout << "Different optimal y:" << optimalObjective.y << " " << global.minY
|
|
<< std::endl;
|
|
}
|
|
//Optimal X values
|
|
results.optimalXNameValuePairs.resize(settings.parameterRanges.size());
|
|
std::vector<double> optimalX(settings.parameterRanges.size());
|
|
for (int optimizationParameterIndex = E;
|
|
optimizationParameterIndex != NumberOfOptimizationParameters;
|
|
optimizationParameterIndex++) {
|
|
if (settings.parameterRanges[optimizationParameterIndex].max
|
|
- settings.parameterRanges[optimizationParameterIndex].min
|
|
< 1e-10) {
|
|
continue;
|
|
}
|
|
if (optimizationParameterIndex < settings.parameterRanges.size() - 2) {
|
|
results.optimalXNameValuePairs[optimizationParameterIndex]
|
|
= std::make_pair(settings.parameterRanges[optimizationParameterIndex].label,
|
|
optimalObjective.x[optimizationParameterIndex]
|
|
* global.parametersInitialValue[optimizationParameterIndex]);
|
|
} else {
|
|
//Hex size and angle are pure values (not multipliers of the initial values)
|
|
results.optimalXNameValuePairs[optimizationParameterIndex]
|
|
= std::make_pair(settings.parameterRanges[optimizationParameterIndex].label,
|
|
optimalObjective.x[optimizationParameterIndex]);
|
|
}
|
|
|
|
optimalX[optimizationParameterIndex] = optimalObjective.x[optimizationParameterIndex];
|
|
#ifdef POLYSCOPE_DEFINED
|
|
std::cout << results.optimalXNameValuePairs[optimizationParameterIndex].first << ":"
|
|
<< optimalX[optimizationParameterIndex] << " ";
|
|
#endif
|
|
}
|
|
#ifdef POLYSCOPE_DEFINED
|
|
std::cout << std::endl;
|
|
#endif
|
|
|
|
results.fullPatternYoungsModulus = youngsModulus;
|
|
// 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(std::vector<double>(optimalObjective.x.begin(),
|
|
optimalObjective.x.end()),
|
|
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]);
|
|
|
|
// results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i];
|
|
//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.objectiveValueHistory = global.objectiveValueHistory;
|
|
results.objectiveValueHistory_iteration = global.objectiveValueHistory_iteration;
|
|
// 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)
|
|
{
|
|
global.objectiveValueHistory.clear();
|
|
global.objectiveValueHistory_iteration.clear();
|
|
global.objectiveValueHistory.reserve(settings.numberOfFunctionCalls / 100);
|
|
global.objectiveValueHistory_iteration.reserve(settings.numberOfFunctionCalls / 100);
|
|
|
|
#if POLYSCOPE_DEFINED
|
|
// global.plotColors.reserve(settings.numberOfFunctionCalls);
|
|
#endif
|
|
#ifdef USE_ENSMALLEN
|
|
#else
|
|
double (*objF)(const dlib::matrix<double, 0, 1> &) = &objective;
|
|
#endif
|
|
enum OptimizationParameterComparisonScenarioIndex {
|
|
AllVar,
|
|
YMAndGeo,
|
|
noYM,
|
|
SplittedMatGeo,
|
|
SplittedYM_before,
|
|
SplittedYM_after,
|
|
SplittedGeoYMMat,
|
|
NumberOfScenarios
|
|
};
|
|
const std::vector<std::vector<std::vector<OptimizationParameterIndex>>> scenarioParameters =
|
|
[&]() {
|
|
std::vector<std::vector<std::vector<OptimizationParameterIndex>>> scenarioParameters(
|
|
NumberOfScenarios);
|
|
scenarioParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
|
|
scenarioParameters[YMAndGeo] = {{R, Theta, E}};
|
|
scenarioParameters[noYM] = {{A, I2, I3, J, R, Theta}};
|
|
scenarioParameters[SplittedMatGeo] = {{E, A, I2, I3, J}, {R, Theta}};
|
|
scenarioParameters[SplittedYM_before] = {{E}, {A, I2, I3, J, R, Theta}};
|
|
scenarioParameters[SplittedYM_after] = {{A, I2, I3, J, R, Theta}, {E}};
|
|
scenarioParameters[SplittedGeoYMMat] = {{R, Theta}, {E}, {A, I2, I3, J}};
|
|
return scenarioParameters;
|
|
}();
|
|
|
|
constexpr OptimizationParameterComparisonScenarioIndex scenario = noYM;
|
|
const std::vector<std::vector<OptimizationParameterIndex>> scenarioParameterGroups
|
|
= scenarioParameters[scenario];
|
|
|
|
//TODO:Ensure that the reduced pattern mesh has the initial parameter values before starting the optimization
|
|
FunctionEvaluation optimalResult;
|
|
for (const std::vector<OptimizationParameterIndex> ¶meterGroup : scenarioParameterGroups) {
|
|
//Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup
|
|
function_updateReducedPattern = [&](const std::vector<double> &x,
|
|
std::shared_ptr<SimulationMesh> &pMesh) {
|
|
for (int xIndex = 0; xIndex < x.size(); xIndex++) {
|
|
const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
|
|
// const double parameterInitialValue=optimizationSettings.parameterRanges[parameterIndex].initialValue;
|
|
const double parameterInitialValue = global.parametersInitialValue[parameterIndex];
|
|
const double parameterNewValue = [&]() {
|
|
if (parameterIndex == R
|
|
|| parameterIndex
|
|
== Theta) { //NOTE:Here I explore the geometry parameters linearly
|
|
return x[xIndex] /*+ parameterInitialValue*/;
|
|
}
|
|
//and the material parameters exponentially(?).TODO: Check what happens if I make all linear
|
|
return x[xIndex] * parameterInitialValue;
|
|
}();
|
|
global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue,
|
|
pMesh);
|
|
}
|
|
pMesh->reset(); //NOTE: I could put this code into each updateParameter function for avoiding unessecary calculations
|
|
};
|
|
|
|
std::vector<double> xMin;
|
|
std::vector<double> xMax;
|
|
xMin.resize(parameterGroup.size());
|
|
xMax.resize(parameterGroup.size());
|
|
for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
|
|
const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
|
|
|
|
xMin[xIndex] = settings.parameterRanges[parameterIndex].min;
|
|
xMax[xIndex] = settings.parameterRanges[parameterIndex].max;
|
|
}
|
|
#ifdef USE_ENSMALLEN
|
|
arma::mat x(parameterGroup.size(), 1);
|
|
for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
|
|
const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
|
|
|
|
x(xIndex, 0) = global.optimizationInitialValue[parameterIndex];
|
|
}
|
|
|
|
global.xMin = xMin;
|
|
global.xMax = xMax;
|
|
|
|
// Create simulated annealing optimizer with default options.
|
|
// The ens::SA<> type can be replaced with any suitable ensmallen optimizer
|
|
// that is able to handle arbitrary functions.
|
|
EnsmallenOptimizationObjective optimizationFunction;
|
|
//Set min max values
|
|
// ens::SA optimizer;
|
|
// ens::CNE optimizer;
|
|
// ens::DE optimizer;
|
|
// ens::SPSA optimizer;
|
|
// arma::mat xMin_arma(global.xMin);
|
|
// arma::mat xMax_arma(global.xMax);
|
|
// ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000);
|
|
ens::LBestPSO optimizer;
|
|
const double minima = optimizer.Optimize(optimizationFunction, x);
|
|
for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
|
|
if (x[xIndex] > global.xMax[xIndex]) {
|
|
x[xIndex] = global.xMax[xIndex];
|
|
} else if (x[xIndex] < global.xMin[xIndex]) {
|
|
x[xIndex] = global.xMin[xIndex];
|
|
}
|
|
}
|
|
|
|
#ifdef POLYSCOPE_DEFINED
|
|
std::cout << "optimal x:"
|
|
<< "\n"
|
|
<< x << std::endl;
|
|
std::cout << "Minima:" << minima << std::endl;
|
|
#endif
|
|
optimalResult.x.clear();
|
|
optimalResult.x.resize(parameterGroup.size());
|
|
std::copy(x.begin(), x.end(), optimalResult.x.begin());
|
|
optimalResult.y = minima;
|
|
|
|
#else
|
|
//Set min max values
|
|
dlib::matrix<double, 0, 1> xMin_dlib(parameterGroup.size());
|
|
dlib::matrix<double, 0, 1> xMax_dlib(parameterGroup.size());
|
|
for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
|
|
const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
|
|
|
|
xMin_dlib(xIndex) = settings.parameterRanges[parameterIndex].min;
|
|
xMax_dlib(xIndex) = settings.parameterRanges[parameterIndex].max;
|
|
}
|
|
const dlib::function_evaluation optimalResult_dlib
|
|
= dlib::find_min_global(objF,
|
|
xMin_dlib,
|
|
xMax_dlib,
|
|
dlib::max_function_calls(settings.numberOfFunctionCalls),
|
|
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 = 10000;
|
|
// dlib::matrix<double, 1, 0> x;
|
|
// dlib::function_evaluation optimalResult_dlib;
|
|
// optimalResult_dlib.x.set_size(parameterGroup.size());
|
|
// for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
|
|
// const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
|
|
// optimalResult_dlib.x(xIndex) = global.optimizationInitialValue[parameterIndex];
|
|
|
|
// std::cout << "xIndex:" << xIndex << std::endl;
|
|
// std::cout << "xInit:" << optimalResult_dlib.x(xIndex) << std::endl;
|
|
// }
|
|
|
|
// optimalResult_dlib.y = dlib::find_min_bobyqa(objF,
|
|
// optimalResult_dlib.x,
|
|
// npt,
|
|
// xMin,
|
|
// xMax,
|
|
// rhobeg,
|
|
// rhoend,
|
|
// bobyqa_maxFunctionCalls);
|
|
// }
|
|
optimalResult.x.clear();
|
|
optimalResult.x.resize(parameterGroup.size());
|
|
std::copy(optimalResult_dlib.x.begin(), optimalResult_dlib.x.end(), optimalResult.x.begin());
|
|
optimalResult.y = optimalResult_dlib.y;
|
|
|
|
#endif
|
|
function_updateReducedPattern(
|
|
std::vector<double>(optimalResult.x.begin(), optimalResult.x.end()),
|
|
global.reducedPatternSimulationJobs[0]
|
|
->pMesh); //TODO: Check if its ok to update only the changed parameters
|
|
// std::cout << "GLOBAL MIN:" << global.minY << std::endl;
|
|
// std::cout << "opt res y:" << optimalResult.y << std::endl;
|
|
}
|
|
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;
|
|
}
|
|
|
|
void search(const std::function<double(const double &)> objectiveFunction,
|
|
const double &targetY,
|
|
double &optimalX,
|
|
double &error,
|
|
const double &epsilon,
|
|
const double &maxIterations)
|
|
{
|
|
error = std::numeric_limits<double>::max();
|
|
int iterationIndex = 0;
|
|
double low = optimalX / 1e-3, high = optimalX * 1e3;
|
|
while (error > epsilon && iterationIndex < maxIterations) {
|
|
const double y = objectiveFunction(optimalX);
|
|
error = std::abs(targetY - y);
|
|
if (y > targetY) {
|
|
high = optimalX;
|
|
} else {
|
|
low = optimalX;
|
|
}
|
|
optimalX = low + (high - low) / 2;
|
|
|
|
iterationIndex++;
|
|
}
|
|
|
|
if (iterationIndex == maxIterations) {
|
|
std::cerr << "Max iterations reached." << std::endl;
|
|
}
|
|
}
|
|
|
|
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;
|
|
}
|
|
|
|
#ifdef USE_ENSMALLEN
|
|
struct ForceMagnitudeOptimization
|
|
{
|
|
std::function<double(const double &)> objectiveFunction;
|
|
ForceMagnitudeOptimization(std::function<double(const double &)> &f) : objectiveFunction(f) {}
|
|
double Evaluate(const arma::mat &x) { return objectiveFunction(x(0, 0)); }
|
|
};
|
|
#endif
|
|
|
|
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;
|
|
}
|
|
|
|
constexpr int maxIterations = 1000;
|
|
minimumError = dlib::find_min_single_variable(objectiveFunction,
|
|
forceMagnitude,
|
|
1e-4,
|
|
1e4,
|
|
forceMagnitudeEpsilon,
|
|
maxIterations);
|
|
#ifdef DLIB_DEFINED
|
|
#else
|
|
// ens::SA<> optimizer;
|
|
// arma::vec lowerBound("0.00001");
|
|
// arma::vec upperBound("10000");
|
|
// ens::LBestPSO optimizer(64, lowerBound, upperBound, maxIterations, 350, forceMagnitudeEpsilon);
|
|
// ForceMagnitudeOptimization f(objectiveFunction); // Create function to be optimized.
|
|
// arma::mat forceMagnitude_mat({forceMagnitude});
|
|
// minimumError = optimizer.Optimize(f, forceMagnitude_mat);
|
|
// std::cout << ReducedPatternOptimization::baseSimulationScenarioNames[scenario] << ": "
|
|
// << optimalObjective << std::endl;
|
|
// forceMagnitude = forceMagnitude_mat(0, 0);
|
|
|
|
// search(objectiveFunction,
|
|
// global.desiredMaxDisplacementValue,
|
|
// forceMagnitude,
|
|
// minimumError,
|
|
// objectiveEpsilon,
|
|
// maxIterations);
|
|
#endif
|
|
|
|
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;
|
|
|
|
global.scenarioEqualizationWeight.resize(
|
|
totalNumberOfSimulationScenarios); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization
|
|
global.scenarioUserWeights.resize(totalNumberOfSimulationScenarios);
|
|
const double userWeightsSum = std::accumulate(baseScenarioWeights.begin(),
|
|
baseScenarioWeights.end(),
|
|
0);
|
|
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);
|
|
const double baseScenarioEqualizationWeight = static_cast<double>(
|
|
numberOfSimulationScenarios)
|
|
/ totalNumberOfSimulationScenarios;
|
|
const double baseScenarioUserWeight = static_cast<double>(baseScenarioWeights[scenario])
|
|
/ userWeightsSum;
|
|
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);
|
|
global.scenarioEqualizationWeight[simulationScenarioIndex]
|
|
= baseScenarioEqualizationWeight;
|
|
global.scenarioUserWeights[simulationScenarioIndex] = baseScenarioUserWeight;
|
|
}
|
|
}
|
|
|
|
#ifdef POLYSCOPE_DEFINED
|
|
std::cout << "Computed full pattern scenario magnitudes" << std::endl;
|
|
#endif
|
|
return scenarios;
|
|
}
|
|
|
|
void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors()
|
|
{
|
|
// m_pFullPatternSimulationMesh->registerForDrawing();
|
|
// m_pFullPatternSimulationMesh->save(std::filesystem::current_path().append("initial.ply"));
|
|
// 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];
|
|
// std::cout << "vi:" << fullPatternVi << std::endl;
|
|
// std::cout << "displacement:" << vertexDisplacement.getTranslation().norm() << std::endl;
|
|
|
|
translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm();
|
|
}
|
|
// global.fullPatternResults[simulationScenarioIndex].saveDeformedModel(
|
|
// std::filesystem::current_path());
|
|
// global.fullPatternResults[simulationScenarioIndex].registerForDrawing();
|
|
// polyscope::show();
|
|
// global.fullPatternResults[simulationScenarioIndex].unregister();
|
|
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];
|
|
|
|
// std::filesystem::path patternMaxForceMagnitudesFilePath(
|
|
// std::filesystem::path(forceMagnitudesDirectoryPath)
|
|
// .append(pFullPatternSimulationJob->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;
|
|
// }
|
|
// }
|
|
#ifdef POLYSCOPE_DEFINED
|
|
LinearSimulationModel linearSimulator;
|
|
SimulationResults fullPatternResults = linearSimulator.executeSimulation(
|
|
pFullPatternSimulationJob);
|
|
|
|
#else
|
|
SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
|
|
simulationSettings);
|
|
#endif
|
|
// if (!fullPatternResults.converged) {
|
|
// DRMSimulationModel::Settings simulationSettings_secondRound;
|
|
// simulationSettings_secondRound.viscousDampingFactor = 2e-3;
|
|
// simulationSettings_secondRound.useKineticDamping = true;
|
|
// simulationSettings.maxDRMIterations = 200000;
|
|
// fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob,
|
|
// simulationSettings_secondRound);
|
|
#ifdef POLYSCOPE_DEFINED
|
|
// std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel()
|
|
// << " used viscous damping." << std::endl;
|
|
#endif
|
|
|
|
if (!fullPatternResults.converged) {
|
|
results.wasSuccessful = false;
|
|
std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel()
|
|
<< " did not converge." << std::endl;
|
|
#ifdef POLYSCOPE_DEFINED
|
|
// DRMSimulationModel::Settings debugSimulationSettings;
|
|
// debugSimulationSettings.debugModeStep = 50;
|
|
// // // debugSimulationSettings.maxDRMIterations = 100000;
|
|
// debugSimulationSettings.shouldDraw = true;
|
|
// // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep;
|
|
// debugSimulationSettings.shouldCreatePlots = true;
|
|
// // // debugSimulationSettings.Dtini = 0.06;
|
|
// debugSimulationSettings.beVerbose = true;
|
|
// debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5;
|
|
// debugSimulationSettings.maxDRMIterations = 100000;
|
|
// debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
|
|
// debugSimulationSettings.viscousDampingFactor = 1e-2;
|
|
// // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3;
|
|
// // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true;
|
|
// auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob,
|
|
// debugSimulationSettings);
|
|
// debugResults.setLabelPrefix("debugResults");
|
|
// debugResults.registerForDrawing();
|
|
// polyscope::show();
|
|
// debugResults.unregister();
|
|
std::filesystem::path outputPath(
|
|
std::filesystem::path("../nonConvergingJobs")
|
|
.append(m_pFullPatternSimulationMesh->getLabel())
|
|
.append("final_" + pFullPatternSimulationJob->getLabel()));
|
|
std::filesystem::create_directories(outputPath);
|
|
pFullPatternSimulationJob->save(outputPath);
|
|
simulationSettings.save(outputPath);
|
|
#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::Epsilon) {
|
|
computeObjectiveValueNormalizationFactors();
|
|
}
|
|
#ifdef POLYSCOPE_DEFINED
|
|
std::cout << "Running reduced model optimization" << std::endl;
|
|
#endif
|
|
runOptimization(optimizationSettings, results);
|
|
results.notes = optimizationNotes;
|
|
}
|