manolas
/
ReducedModelOptimization
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Created simulation scenarios.

This commit is contained in:
Iason 2020-12-09 17:58:48 +02:00
parent b5ef56dcb7
commit 4d4aa346b5
3 changed files with 420 additions and 356 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

181
src/main.cpp
View File

@ -6,18 +6,42 @@
#include "polyscope/polyscope.h"
#include "reducedmodeloptimizer.hpp"
#include "simulationhistoryplotter.hpp"
#include "trianglepattterntopology.hpp"
#include <chrono>
#include <filesystem>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vcg/complex/algorithms/update/position.h>
// void scale(const std::vector<size_t>& numberOfNodesPerSlot,
// FlatPattern &pattern) {
// const double desiredSize = 0.0025; // center to boundary
// const size_t interfaceSlotIndex = 4; // bottom edge
// std::unordered_map<size_t, size_t> nodeToSlot;
// std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
// FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot,
// nodeToSlot); FlatPatternTopology::constructSlotToNodeMap(nodeToSlot,
// slotToNode); 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());
// const double currentSize =
// (pattern.vert[baseTriangleInterfaceVi].cP() - pattern.vert[0].cP())
// .Norm();
// const double scaleFactor = desiredSize / currentSize;
// vcg::tri::UpdatePosition<FlatPattern>::Scale(full, scaleFactor);
//}
int main(int argc, char *argv[]) {
// FlatPattern pattern("/home/iason/Models/valid_6777.ply");
// FlatPattern pattern("/home/iason/Models/simple_beam_paper_example.ply");
// pattern.savePly("fannedValid.ply");
registerWorldAxes();
// Create reduced models
const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
std::vector<vcg::Point2i> singleBarReducedModelEdges{vcg::Point2i(0, 3)};
@ -25,50 +49,56 @@ int main(int argc, char *argv[]) {
singleBarReducedModelEdges);
singleBarReducedModel.setLabel("Single bar reduced model");
std::vector<vcg::Point2i> CCWreducedModelEdges{vcg::Point2i(1, 5),
vcg::Point2i(3, 5)};
FlatPattern CWReducedModel(numberOfNodesPerSlot, CCWreducedModelEdges);
CWReducedModel.setLabel("CW reduced model");
// std::vector<vcg::Point2i> CCWreducedModelEdges{vcg::Point2i(1, 5),
// vcg::Point2i(3, 5)};
// FlatPattern CWReducedModel(numberOfNodesPerSlot, CCWreducedModelEdges);
// CWReducedModel.setLabel("CW reduced model");
std::vector<vcg::Point2i> CWreducedModelEdges{vcg::Point2i(1, 5),
vcg::Point2i(3, 1)};
FlatPattern CCWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges);
CCWReducedModel.setLabel("CCW reduced model");
// std::vector<vcg::Point2i> CWreducedModelEdges{vcg::Point2i(1, 5),
// vcg::Point2i(3, 1)};
// FlatPattern CCWReducedModel(numberOfNodesPerSlot, CWreducedModelEdges);
// CCWReducedModel.setLabel("CCW reduced model");
std::vector<FlatPattern *> reducedModels{&singleBarReducedModel,
&CWReducedModel, &CCWReducedModel};
// std::vector<FlatPattern *> reducedModels{&singleBarReducedModel,
// &CWReducedModel,
// &CCWReducedModel};
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
std::string fullPatternsTestSetDirectory =
"/home/iason/Models/TestSet_validPatterns";
for (const auto &entry :
filesystem::directory_iterator(fullPatternsTestSetDirectory)) {
const auto filepath =
// entry.path();
std::filesystem::path("/home/iason/Models/valid_6777.ply");
const auto filepathString = filepath.string();
// Use only the base triangle version
const std::string tiledSuffix = "_tiled.ply";
if (filepathString.compare(filepathString.size() - tiledSuffix.size(),
tiledSuffix.size(), tiledSuffix) == 0) {
continue;
}
FlatPattern pattern(filepathString);
pattern.setLabel(filepath.stem().string());
std::cout << "Testing Pattern:" << filepathString << std::endl;
for (FlatPattern *pReducedModel : reducedModels) {
pReducedModel = reducedModels[0];
std::unordered_set<size_t> optimizationExcludedEi;
if (pReducedModel !=
reducedModels[0]) { // assumes that the singleBar reduced model is the
// first in the reducedModels vector
optimizationExcludedEi.insert(0);
}
optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi);
optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi);
Eigen::VectorXd optimalParameters = optimizer.optimize();
}
}
// ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
// std::string fullPatternsTestSetDirectory =
// "/home/iason/Models/TestSet_validPatterns";
// for (const auto &entry :
// filesystem::directory_iterator(fullPatternsTestSetDirectory)) {
// const auto filepath =
// std::filesystem::path("/home/iason/Models/valid_6777.ply");
// // std::filesystem::path(
// // "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
// // "1v_0v_2e_1e_1c_6fan/3/Valid/431.ply");
// // std::filesystem::path(
// // "/home/iason/Models/TestSet_validPatterns/865.ply");
// // entry.path();
// const auto filepathString = filepath.string();
// // Use only the base triangle version
// const std::string tiledSuffix = "_tiled.ply";
// if (filepathString.compare(filepathString.size() - tiledSuffix.size(),
// tiledSuffix.size(), tiledSuffix) == 0) {
// continue;
// }
// FlatPattern pattern(filepathString);
// pattern.setLabel(filepath.stem().string());
// std::cout << "Testing Pattern:" << filepathString << std::endl;
// for (FlatPattern *pReducedModel : reducedModels) {
// // pReducedModel = reducedModels[1];
// std::unordered_set<size_t> optimizationExcludedEi;
// if (pReducedModel !=
// reducedModels[0]) { // assumes that the singleBar reduced model is
// // the
// // first in the reducedModels vector
// optimizationExcludedEi.insert(0);
// }
// optimizer.initialize(pattern, *pReducedModel, optimizationExcludedEi);
// Eigen::VectorXd optimalParameters = optimizer.optimize();
// }
// }
// // Full model simulation
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
@ -135,43 +165,34 @@ int main(int argc, char *argv[]) {
// stiffnessFactor *= 1.5;
// }
// Beam
// VCGEdgeMesh beamMesh;
// beamMesh.loadFromPly("/home/iason/Models/simple_beam_model_2elem_1m.ply");
// const VertexIndex reducedModelOpposite_vi =
// std::ceil(reducedModelSimulationJob.mesh->VN() / 2.0);
// auto v0 =
// reducedModelSimulationJob.mesh->vert[reducedModelOpposite_vi].cP() -
// reducedModelSimulationJob.mesh->vert[1].cP();
// auto v1 = beamMesh.vert[2].cP() - beamMesh.vert[0].cP();
// vcg::Matrix44d R;
// const double rotationTheta =
// std::asin((v0 ^ v1).Norm() / (v0.Norm() * v1.Norm()));
// R.SetRotateRad(-rotationTheta, v0 ^ v1);
// vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(beamMesh, R);
// vcg::tri::UpdatePosition<VCGEdgeMesh>::Scale(beamMesh, v0.Norm() /
// v1.Norm()); vcg::tri::UpdatePosition<VCGEdgeMesh>::Translate(
// beamMesh, reducedModelSimulationJob.mesh->vert[1].cP());
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
// beamFixedVertices;
// beamFixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// std::unordered_map<VertexIndex, Eigen::Vector3d> beamNodalForces{
// {2, Eigen::Vector3d(0, 0, 1000)}};
// SimulationJob beamSimulationJob{std::make_shared<ElementalMesh>(beamMesh),
// beamFixedVertices,
// beamNodalForces,
// {}};
// // for (EdgeIndex ei = 0; ei < beamSimulationJob.mesh->EN(); ei++) {
// // BeamFormFinder::Element &e = beamSimulationJob.mesh->elements[ei];
// // e.properties.A = stifnessVector[0];
// // e.properties.J = stifnessVector[1];
// // e.properties.I2 = stifnessVector[2];
// // e.properties.I3 = stifnessVector[3];
// // }
// // beamSimulationJob.draw();
// SimulationResults beamSimulationResults =
// formFinder.executeSimulation(beamSimulationJob);
// beamSimulationResults.simulationLabel = "Beam";
// beamSimulationResults.draw(beamSimulationJob);
// Beam
VCGEdgeMesh mesh;
mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
FormFinder formFinder;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
std::unordered_map<VertexIndex, Vector6d> nodalForces{
{10, Vector6d({0, 0, 0, 0.00001, 0, 0})}};
// Forced displacements
std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
std::unordered_map<size_t, VectorType> nodalForcedNormal;
// CoordType v = (mesh.vert[10].cP() - mesh.vert[0].cP()).Normalize();
// nodalForcedNormal[10] = v;
fixedVertices[10] = {0, 1, 2};
SimulationJob beamSimulationJob{std::make_shared<SimulationMesh>(mesh),
fixedVertices, nodalForces,
nodalForcedDisplacements, nodalForcedNormal};
beamSimulationJob.mesh->setBeamMaterial(0.3, 1);
beamSimulationJob.mesh->setBeamCrossSection(
RectangularBeamDimensions{0.002, 0.002});
// registerWorldAxes();
SimulationResults beamSimulationResults =
formFinder.executeSimulation(beamSimulationJob, true, true);
beamSimulationResults.simulationLabel = "Beam";
beamSimulationResults.registerForDrawing(beamSimulationJob);
polyscope::show();
return 0;
}

592
src/reducedmodeloptimizer.cpp
View File

@ -2,126 +2,137 @@
#include "bobyqa.h"
#include "flatpattern.hpp"
#include "gradientDescent.h"
#include "matplot/matplot.h"
#include "simulationhistoryplotter.hpp"
#include "trianglepattterntopology.hpp"
const bool gShouldDraw = true;
size_t g_numberOfOptimizationRounds{0};
FormFinder simulator;
Eigen::MatrixX3d g_optimalReducedModelDisplacements;
SimulationJob gReducedPatternSimulationJob;
SimulationJob g_reducedPatternSimulationJob;
std::unordered_map<ReducedModelVertexIndex, FullModelVertexIndex>
g_reducedToFullInterfaceViMap;
matplot::line_handle gPlotHandle;
std::vector<double> gObjectiveValueHistory;
Eigen::Vector4d g_initialX;
std::unordered_set<size_t> g_reducedPatternExludedEdges;
Eigen::MatrixX4d g_initialStiffnessFactors;
double g_initialParameters;
struct OptimizationCallback {
double operator()(const size_t &iterations, const Eigen::VectorXd &x,
const double &fval, Eigen::VectorXd &gradient) const {
// run simulation
// SimulationResults reducedModelResults =
// simulator.executeSimulation(reducedModelSimulationJob);
// reducedModelResults.draw(reducedModelSimulationJob);
gObjectiveValueHistory.push_back(fval);
auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
gObjectiveValueHistory.size());
gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
// const std::string plotImageFilename = "objectivePlot.png";
// matplot::save(plotImageFilename);
// if (numberOfOptimizationRounds % 30 == 0) {
// std::filesystem::copy_file(
// std::filesystem::path(plotImageFilename),
// std::filesystem::path("objectivePlot_copy.png"));
// struct OptimizationCallback {
// double operator()(const size_t &iterations, const Eigen::VectorXd &x,
// const double &fval, Eigen::VectorXd &gradient) const {
// // run simulation
// // SimulationResults reducedModelResults =
// // simulator.executeSimulation(reducedModelSimulationJob);
// // reducedModelResults.draw(reducedModelSimulationJob);
// gObjectiveValueHistory.push_back(fval);
// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
// gObjectiveValueHistory.size());
// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
// // const std::string plotImageFilename = "objectivePlot.png";
// // matplot::save(plotImageFilename);
// // if (numberOfOptimizationRounds % 30 == 0) {
// // std::filesystem::copy_file(
// // std::filesystem::path(plotImageFilename),
// // std::filesystem::path("objectivePlot_copy.png"));
// // }
// // std::stringstream ss;
// // ss << x;
// // reducedModelResults.simulationLabel = ss.str();
// // SimulationResultsReporter resultsReporter;
// // resultsReporter.reportResults(
// // {reducedModelResults},
// // std::filesystem::current_path().append("Results"));
// return true;
// }
//};
// struct Objective {
// double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const {
// assert(x.rows() == 4);
// // drawSimulationJob(simulationJob);
// // Set mesh from x
// std::shared_ptr<SimulationMesh> reducedModel =
// g_reducedPatternSimulationJob.mesh;
// for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
// if (g_reducedPatternExludedEdges.contains(ei)) {
// continue;
// }
// Element &e = reducedModel->elements[ei];
// e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * x(0);
// e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * x(1);
// e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * x(2);
// e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * x(3);
// }
// // run simulation
// SimulationResults reducedModelResults =
// simulator.executeSimulation(g_reducedPatternSimulationJob);
// // std::stringstream ss;
// // ss << x;
// // reducedModelResults.simulationLabel = ss.str();
// // SimulationResultsReporter resultsReporter;
// // resultsReporter.reportResults(
// // {reducedModelResults},
// // std::filesystem::current_path().append("Results"));
// // compute error and return it
// double error = 0;
// for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
// VertexIndex reducedModelVi = reducedFullViPair.first;
// Eigen::Vector3d vertexDisplacement(
// reducedModelResults.displacements[reducedModelVi][0],
// reducedModelResults.displacements[reducedModelVi][1],
// reducedModelResults.displacements[reducedModelVi][2]);
// Eigen::Vector3d errorVector =
// Eigen::Vector3d(
// g_optimalReducedModelDisplacements.row(reducedModelVi)) -
// vertexDisplacement;
// error += errorVector.norm();
// }
// return error;
// }
//};
double ReducedModelOptimizer::objective(long n, const double *x) {
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
std::cout << "x[" + std::to_string(parameterIndex) + "]="
<< x[parameterIndex] << std::endl;
}
for (EdgeIndex ei = 0; ei < g_reducedPatternSimulationJob.mesh->EN(); ei++) {
Element &e = g_reducedPatternSimulationJob.mesh->elements[ei];
// if (g_reducedPatternExludedEdges.contains(ei)) {
// continue;
// }
// std::stringstream ss;
// ss << x;
// reducedModelResults.simulationLabel = ss.str();
// SimulationResultsReporter resultsReporter;
// resultsReporter.reportResults(
// {reducedModelResults},
// std::filesystem::current_path().append("Results"));
return true;
}
};
struct Objective {
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const {
assert(x.rows() == 4);
// drawSimulationJob(simulationJob);
// Set mesh from x
std::shared_ptr<SimulationMesh> reducedModel =
gReducedPatternSimulationJob.mesh;
for (EdgeIndex ei = 0; ei < reducedModel->EN(); ei++) {
if (g_reducedPatternExludedEdges.contains(ei)) {
continue;
}
Element &e = reducedModel->elements[ei];
e.axialConstFactor *= x(0);
e.torsionConstFactor *= x(1);
e.firstBendingConstFactor *= x(2);
e.secondBendingConstFactor *= x(3);
}
// run simulation
SimulationResults reducedModelResults =
simulator.executeSimulation(gReducedPatternSimulationJob);
// std::stringstream ss;
// ss << x;
// reducedModelResults.simulationLabel = ss.str();
// SimulationResultsReporter resultsReporter;
// resultsReporter.reportResults(
// {reducedModelResults},
// std::filesystem::current_path().append("Results"));
// compute error and return it
double error = 0;
for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
VertexIndex reducedModelVi = reducedFullViPair.first;
Eigen::Vector3d vertexDisplacement(
reducedModelResults.displacements[reducedModelVi][0],
reducedModelResults.displacements[reducedModelVi][1],
reducedModelResults.displacements[reducedModelVi][2]);
Eigen::Vector3d errorVector =
Eigen::Vector3d(
g_optimalReducedModelDisplacements.row(reducedModelVi)) -
vertexDisplacement;
error += errorVector.norm();
}
return error;
}
};
double objective(long n, const double *x) {
Eigen::VectorXd eigenX(n, 1);
for (size_t xi = 0; xi < n; xi++) {
eigenX(xi) = x[xi];
}
std::shared_ptr<SimulationMesh> reducedPattern =
gReducedPatternSimulationJob.mesh;
for (EdgeIndex ei = 0; ei < reducedPattern->EN(); ei++) {
Element &e = reducedPattern->elements[ei];
if (g_reducedPatternExludedEdges.contains(ei)) {
continue;
}
e.axialConstFactor = g_initialStiffnessFactors(ei, 0) * eigenX(0);
e.torsionConstFactor = g_initialStiffnessFactors(ei, 1) * eigenX(1);
e.firstBendingConstFactor = g_initialStiffnessFactors(ei, 2) * eigenX(2);
e.secondBendingConstFactor = g_initialStiffnessFactors(ei, 3) * eigenX(3);
e.properties.E = g_initialParameters * x[0];
// e.properties.G = g_initialParameters(1) * x[1];
// e.properties.A = g_initialParameters(0, 0) * x[0];
// e.properties.J = g_initialParameters(0, 1) * x[1];
// e.properties.I2 = g_initialParameters(0, 2) * x[2];
// e.properties.I3 = g_initialParameters(0, 3) * x[3];
e.properties.G = e.properties.E / (2 * (1 + 0.3));
e.axialConstFactor = e.properties.E * e.properties.A / e.initialLength;
e.torsionConstFactor = e.properties.G * e.properties.J / e.initialLength;
e.firstBendingConstFactor =
2 * e.properties.E * e.properties.I2 / e.initialLength;
e.secondBendingConstFactor =
2 * e.properties.E * e.properties.I3 / e.initialLength;
}
// run simulation
SimulationResults reducedModelResults =
simulator.executeSimulation(gReducedPatternSimulationJob, false, false);
simulator.executeSimulation(g_reducedPatternSimulationJob, false, false);
// compute error and return it
double error = 0;
for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
VertexIndex reducedModelVi = reducedFullViPair.first;
// const auto pos =
// g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
// std::cout << "Interface vi " << reducedModelVi << " is at position "
// << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
Eigen::Vector3d vertexDisplacement(
reducedModelResults.displacements[reducedModelVi][0],
reducedModelResults.displacements[reducedModelVi][1],
@ -130,6 +141,7 @@ double objective(long n, const double *x) {
Eigen::Vector3d(
g_optimalReducedModelDisplacements.row(reducedModelVi)) -
vertexDisplacement;
// error += errorVector.squaredNorm();
error += errorVector.norm();
}
@ -267,7 +279,12 @@ void ReducedModelOptimizer::computeMaps(
void ReducedModelOptimizer::createSimulationMeshes(FlatPattern &fullModel,
FlatPattern &reducedModel) {
pReducedModelElementalMesh = std::make_shared<SimulationMesh>(reducedModel);
pReducedModelElementalMesh->setBeamCrossSection(
CrossSectionType{0.002, 0.002});
pReducedModelElementalMesh->setBeamMaterial(0.3, 1);
pFullModelElementalMesh = std::make_shared<SimulationMesh>(fullModel);
pFullModelElementalMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002});
pFullModelElementalMesh->setBeamMaterial(0.3, 1);
}
ReducedModelOptimizer::ReducedModelOptimizer(
@ -296,22 +313,29 @@ void ReducedModelOptimizer::initialize(
}
void ReducedModelOptimizer::initializeStiffnesses() {
g_initialStiffnessFactors.resize(pReducedModelElementalMesh->EN(), 4);
// g_initialParameters.resize(1, 2);
// Save save the beam stiffnesses
for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
Element &e = pReducedModelElementalMesh->elements[ei];
// if (g_reducedPatternExludedEdges.contains(ei)) {
// const double stiffnessFactor = 1;
// e.axialConstFactor *= stiffnessFactor;
// e.torsionConstFactor *= stiffnessFactor;
// e.firstBendingConstFactor *= stiffnessFactor;
// e.secondBendingConstFactor *= stiffnessFactor;
// }
g_initialStiffnessFactors(ei, 0) = e.axialConstFactor;
g_initialStiffnessFactors(ei, 1) = e.torsionConstFactor;
g_initialStiffnessFactors(ei, 2) = e.firstBendingConstFactor;
g_initialStiffnessFactors(ei, 3) = e.secondBendingConstFactor;
}
// for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
// Element &e = pReducedModelElementalMesh->elements[ei];
// if (g_reducedPatternExludedEdges.contains(ei)) {
// const double stiffnessFactor = 5;
// e.axialConstFactor *= stiffnessFactor;
// e.torsionConstFactor *= stiffnessFactor;
// e.firstBendingConstFactor *= stiffnessFactor;
// e.secondBendingConstFactor *= stiffnessFactor;
// }
g_initialParameters = pReducedModelElementalMesh->elements[0].properties.E;
// g_initialParameters(1) =
// pReducedModelElementalMesh->elements[0].properties.G;
// g_initialParameters(0, 0) =
// pReducedModelElementalMesh->elements[0].properties.A;
// g_initialParameters(0, 1) =
// pReducedModelElementalMesh->elements[0].properties.J;
// g_initialParameters(0, 2) =
// pReducedModelElementalMesh->elements[0].properties.I2;
// g_initialParameters(0, 3) =
// pReducedModelElementalMesh->elements[0].properties.I3;
// }
}
void ReducedModelOptimizer::computeReducedModelSimulationJob(
@ -331,10 +355,19 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob(
fullModelNodalForce.first)] = fullModelNodalForce.second;
}
std::unordered_map<VertexIndex, VectorType> reducedModelNodalForcedNormals;
for (auto fullModelNodalForcedRotation :
simulationJobOfFullModel.nodalForcedNormals) {
reducedModelNodalForcedNormals[m_fullToReducedInterfaceViMap.at(
fullModelNodalForcedRotation.first)] =
fullModelNodalForcedRotation.second;
}
simulationJobOfReducedModel = SimulationJob{pReducedModelElementalMesh,
reducedModelFixedVertices,
reducedModelNodalForces,
{}};
{},
reducedModelNodalForcedNormals};
}
SimulationJob ReducedModelOptimizer::getReducedSimulationJob(
@ -349,6 +382,7 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {
optimalDisplacementsOfReducedModel.resize(0, 3);
optimalDisplacementsOfReducedModel.resize(pReducedModelElementalMesh->VN(),
3);
optimalDisplacementsOfReducedModel.setZero(
@ -367,115 +401,88 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
}
Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
const SimulationJob &fullModelSimulationJob) {
const SimulationJob &fullPatternSimulationJob) {
gObjectiveValueHistory.clear();
// fullPatternSimulationJob.registerForDrawing();
// polyscope::show();
fullPatternSimulationJob.mesh->savePly(
"Fanned_" + pFullModelElementalMesh->getLabel() + ".ply");
SimulationResults fullModelResults =
simulator.executeSimulation(fullModelSimulationJob, false, false);
simulator.executeSimulation(fullPatternSimulationJob, false, true);
fullModelResults.simulationLabel = "fullModel";
computeDesiredReducedModelDisplacements(fullModelResults,
g_optimalReducedModelDisplacements);
computeReducedModelSimulationJob(fullModelSimulationJob,
gReducedPatternSimulationJob);
computeReducedModelSimulationJob(fullPatternSimulationJob,
g_reducedPatternSimulationJob);
// fullModelSimulationJob.registerForDrawing();
// polyscope::show();
// gReducedPatternSimulationJob.registerForDrawing();
// polyscope::show();
fullModelResults.registerForDrawing(fullModelSimulationJob);
fullModelResults.registerForDrawing(fullPatternSimulationJob);
polyscope::show();
double (*pObjectiveFunction)(long, const double *) = &objective;
const size_t n = 1;
// g_initialParameters.rows();
const size_t npt = ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2;
assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
assert(npt <= 2 * n + 1 && "The choice of the number of interpolation "
"conditions is not recommended.");
// Set initial guess of solution
Eigen::VectorXd initialGuess(4);
const double stifnessFactor = 1;
initialGuess(0) = stifnessFactor;
initialGuess(1) = stifnessFactor;
initialGuess(2) = stifnessFactor;
initialGuess(3) = stifnessFactor;
const double initialGuess = 1;
std::vector<double> x(n, initialGuess);
// {1, 5.9277};
// {initialGuess(0), initialGuess(1), initialGuess(2),
// initialGuess(3)};
const double xMin = 1;
const double xMax = 100;
// assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
// return d >= xMax || d <= xMin;
// }));
std::vector<double> xLow(x.size(), xMin);
std::vector<double> xUpper(x.size(), xMax);
const double maxX = *std::max_element(
x.begin(), x.end(),
[](const double &a, const double &b) { return abs(a) < abs(b); });
const double rhobeg = std::min(0.95, 0.2 * maxX);
// const double rhobeg = 10;
const double rhoend = rhobeg * 1e-6;
const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
std::vector<double> w(wSize);
// const size_t maxFun = 10 * (x.size() ^ 2);
const size_t maxFun = 120;
bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
rhobeg, rhoend, maxFun, w.data());
const bool useGradientDescent = false;
if (useGradientDescent) {
SimulationResults reducedModelOptimizedResults =
simulator.executeSimulation(g_reducedPatternSimulationJob);
// gdc::GradientDescent<double, Objective,
// gdc::DecreaseBacktracking<double>,
// OptimizationCallback>
gdc::GradientDescent<double, Objective, gdc::BarzilaiBorwein<double>,
OptimizationCallback>
// gdc::GradientDescent<double, Objective,
// gdc::DecreaseBacktracking<double>,
// OptimizationCallback>
// gdc::GradientDescent<double, Objective,
// gdc::DecreaseBacktracking<double>,
// OptimizationCallback>
optimizer;
// Turn verbosity on, so the optimizer prints status updates after each
// iteration.
optimizer.setVerbosity(1);
// Set initial guess.
matplot::xlabel("Optimization iterations");
matplot::ylabel("Objective value");
// matplot::figure(false);
matplot::grid(matplot::on);
// Start the optimization
auto result = optimizer.minimize(initialGuess);
std::cout << "Done! Converged: " << (result.converged ? "true" : "false")
<< " Iterations: " << result.iterations << std::endl;
// do something with final function value
std::cout << "Final fval: " << result.fval << std::endl;
// do something with final x-value
std::cout << "Final xval: " << result.xval.transpose() << std::endl;
SimulationResults reducedModelOptimizedResults =
simulator.executeSimulation(gReducedPatternSimulationJob);
reducedModelOptimizedResults.simulationLabel = "reducedModel";
reducedModelOptimizedResults.registerForDrawing(
gReducedPatternSimulationJob);
return result.xval;
} else { // use bobyqa
double (*pObjectiveFunction)(long, const double *) = &objective;
const size_t n = 4;
const size_t npt = 8;
assert(npt <= (n + 1) * (n + 2) / 2 && npt >= n + 2);
assert(npt < 2 * n + 1 && "The choice of the number of interpolation "
"conditions is not recommended.");
std::vector<double> x
// {1.03424, 0.998456, 0.619916, -0.202997};
{initialGuess(0), initialGuess(1), initialGuess(2), initialGuess(3)};
std::vector<double> xLow(x.size(), -100);
std::vector<double> xUpper(x.size(), 100);
const double maxX = *std::max_element(
x.begin(), x.end(),
[](const double &a, const double &b) { return abs(a) < abs(b); });
const double rhobeg = std::min(0.95, 0.2 * maxX);
const double rhoend = rhobeg * 1e-6;
const size_t wSize = (npt + 5) * (npt + n) + 3 * n * (n + 5) / 2;
std::vector<double> w(wSize);
bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
rhobeg, rhoend, 100, w.data());
std::cout << "Final objective value:" << objective(n, x.data())
<< std::endl;
Eigen::VectorXd eigenX(x.size(), 1);
for (size_t xi = 0; xi < x.size(); xi++) {
eigenX(xi) = x[xi];
}
SimulationResults reducedModelOptimizedResults =
simulator.executeSimulation(gReducedPatternSimulationJob);
reducedModelOptimizedResults.simulationLabel = "reducedModel";
reducedModelOptimizedResults.registerForDrawing(
gReducedPatternSimulationJob);
polyscope::show();
return eigenX;
double error = 0;
for (const auto reducedToFullViPair : g_reducedToFullInterfaceViMap) {
const size_t reducedInterfaceVi = reducedToFullViPair.first;
error +=
(Eigen::Vector3d(
g_optimalReducedModelDisplacements.row(reducedInterfaceVi)) -
Eigen::Vector3d(
reducedModelOptimizedResults.displacements[reducedInterfaceVi][0],
reducedModelOptimizedResults.displacements[reducedInterfaceVi][1],
reducedModelOptimizedResults.displacements[reducedInterfaceVi][2]))
.norm();
}
std::cout << "Final objective value:" << error << std::endl;
reducedModelOptimizedResults.simulationLabel = "reducedModel";
reducedModelOptimizedResults.registerForDrawing(
g_reducedPatternSimulationJob);
polyscope::show();
Eigen::VectorXd eigenX(x.size(), 1);
for (size_t xi = 0; xi < x.size(); xi++) {
eigenX(xi) = x[xi];
}
return eigenX;
}
std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
@ -483,7 +490,10 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
std::vector<SimulationJob> scenarios;
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalForces;
const double forceMagnitude = 250;
std::unordered_map<VertexIndex, VectorType> nodalForcedNormals;
const double forceMagnitude = 1;
// Assuming the patterns lays on the x-y plane
const CoordType patternPlaneNormal(0, 0, 1);
// // Axial
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
@ -493,82 +503,100 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
// nodalForces[viPair.first] = Vector6d({forceDirection[0],
// forceDirection[1],
// forceDirection[2], 0, 0, 0}) *
// forceMagnitude;
// forceMagnitude * 10;
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// // In-plane Bending
// // Assuming the patterns lay on the x-y plane
// const CoordType patternPlane(0, 0, 1);
// fixedVertices.clear();
// nodalForces.clear();
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType forceDirection = patternPlane ^ v;
// CoordType forceDirection = (v ^ patternPlaneNormal).Normalize();
// nodalForces[viPair.first] = Vector6d({forceDirection[0],
// forceDirection[1],
// forceDirection[2], 0, 0, 0}) *
// forceMagnitude;
// 0.40 * forceMagnitude;
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// // Torsion
// {
// fixedVertices.clear();
// nodalForces.clear();
// for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
// viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++)
// {
// const auto &viPair = *viPairIt;
// if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() -
// pMesh->vert[viPair.second].cP())
// .Normalize();
// nodalForces[viPair.first] =
// Vector6d({0, 0, 0, v[1], v[0], 0}) * forceMagnitude;
// } else {
// fixedVertices[viPair.first] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// }
// Out-of-plane bending . Pull towards Z
fixedVertices.clear();
nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0});
fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
}
scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// // Dou??
// fixedVertices.clear();
// nodalForces.clear();
// for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
// viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
// const auto &viPair = *viPairIt;
// if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType normalVDerivativeDir = (v ^ patternPlaneNormal).Normalize();
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
// nodalForcedNormals[viPair.first] = normalVDerivativeDir;
// } else {
// fixedVertices[viPair.first] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back(
// {pMesh, fixedVertices, nodalForces, {}, nodalForcedNormals});
// // Out - of - plane bending.Pull towards Z
// fixedVertices.clear();
// nodalForces.clear();
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) *
// 2; fixedVertices[viPair.second] =
// std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// Double
// fixedVertices.clear();
// nodalForces.clear();
// nodalForcedNormals.clear();
// int counter = 0;
// for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType momentDirection = patternPlane ^ v;
// nodalForces[viPair.first] =
// Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) *
// forceMagnitude;
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
// nodalForcedNormals[viPair.first] = -v;
// if (counter++ == 1) {
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0};
// } else {
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
// }
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// scenarios.push_back(
// {pMesh, fixedVertices, nodalForces, {}, nodalForcedNormals});
// // Saddle
// Double
fixedVertices.clear();
nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
.Normalize();
CoordType momentDirection = patternPlaneNormal ^ v;
nodalForces[viPair.first] =
Vector6d({0, 0, 0, momentDirection[1], momentDirection[0], 0}) *
forceMagnitude;
fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
}
scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
//// Saddle using moments
// fixedVertices.clear();
// nodalForces.clear();
// for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
@ -577,38 +605,50 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
// CoordType v =
// (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
// .Normalize();
// CoordType momentDirection = patternPlane ^ v;
// // CoordType momentDirection = (patternPlane ^ v).Normalize();
// if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
// nodalForces[viPair.first] =
// Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * 3
// * forceMagnitude;
// Vector6d({0, 0, 0, v[0], v[1], 0}) * 5 * forceMagnitude;
// nodalForces[viPair.second] =
// Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) * 3
// *
// (-forceMagnitude);
// Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 5 * forceMagnitude;
// } else {
// fixedVertices[viPair.first] = std::unordered_set<DoFType>{2};
// // nodalForces[viPair.first] =
// // Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 2 * forceMagnitude;
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{2};
// nodalForces[viPair.first] =
// Vector6d({0, 0, 0, momentDirection[0], momentDirection[1], 0}) *
// (-forceMagnitude);
// fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
// Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 2 * forceMagnitude;
// nodalForces[viPair.second] =
// Vector6d({0, 0, 0, v[0], v[1], 0}) * 2 * forceMagnitude;
// }
// }
// scenarios.push_back({pMesh, fixedVertices, nodalForces, {}});
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
// saddleFixedVertices;
// // saddle_fixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2};
// saddleFixedVertices[7] = std::unordered_set<DoFType>{0, 1, 2};
// saddleFixedVertices[11] = std::unordered_set<DoFType>{0, 1, 2};
// // saddle_fixedVertices[15] = std::unordered_set<DoFType>{0, 1, 2};
// // saddle_fixedVertices[19] = std::unordered_set<DoFType>{0, 1, 2};
// // saddle_fixedVertices[23] = std::unordered_set<DoFType>{0, 1, 2};
// std::unordered_map<VertexIndex, Vector6d> saddleNodalForces{
// {15, {0, 0, 0, 0, -4 * 90, 0}}, {3, {0, 0, 0, 0, 4 * 90, 0}},
// {7, {0, 0, 0, 4 * 70, 0, 0}}, {11, {0, 0, 0, 4 * 70, 0, 0}},
// {19, {0, 0, 0, -4 * 70, 0, 0}}, {23, {0, 0, 0, -4 * 70, 0, 0}}};
// scenarios.push_back({pMesh, saddleFixedVertices, saddleNodalForces, {}});
//// Saddle using forced normals
fixedVertices.clear();
nodalForces.clear();
nodalForcedNormals.clear();
for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin();
viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) {
const auto &viPair = *viPairIt;
CoordType v =
(pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP())
.Normalize();
// CoordType momentDirection = (patternPlane ^ v).Normalize();
if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) {
nodalForcedNormals[viPair.first] = v;
nodalForcedNormals[viPair.second] = -v;
} else {
nodalForcedNormals[viPair.first] = -v;
nodalForcedNormals[viPair.second] = v;
}
fixedVertices[viPair.first] = std::unordered_set<DoFType>{0, 1, 2};
fixedVertices[viPair.second] = std::unordered_set<DoFType>{0, 1, 2};
}
scenarios.push_back(
{pMesh, fixedVertices, nodalForces, {}, nodalForcedNormals});
return scenarios;
}

3
src/reducedmodeloptimizer.hpp
View File

@ -4,6 +4,7 @@
#include "beamformfinder.hpp"
#include "edgemesh.hpp"
#include "elementalmesh.hpp"
#include "matplot/matplot.h"
#include <Eigen/Dense>
using FullModelVertexIndex = VertexIndex;
@ -47,6 +48,8 @@ private:
void createSimulationMeshes(FlatPattern &fullModel,
FlatPattern &reducedModel);
void initializeStiffnesses();
static double objective(long n, const double *x);
};
#endif // REDUCEDMODELOPTIMIZER_HPP