Optimizes the reduced model for a set of simulation scenarios instead of a single scenario
This commit is contained in:
parent
9fe6b83ee8
commit
81ac6c20b8
|
@ -0,0 +1,4 @@
|
||||||
|
/build
|
||||||
|
/TestSet
|
||||||
|
*.user
|
||||||
|
|
|
@ -75,7 +75,6 @@ target_include_directories(${PROJECT_NAME}
|
||||||
PRIVATE ${BOBYQA_SOURCE_DIR}/include/
|
PRIVATE ${BOBYQA_SOURCE_DIR}/include/
|
||||||
)
|
)
|
||||||
|
|
||||||
#Use C++17
|
|
||||||
if(MSVC)
|
if(MSVC)
|
||||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++20")
|
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /std:c++20")
|
||||||
else(MSVC)
|
else(MSVC)
|
||||||
|
|
11
src/main.cpp
11
src/main.cpp
|
@ -93,7 +93,6 @@ int main(int argc, char *argv[]) {
|
||||||
// "/home/iason/Models/TestSet_validPatterns/59.ply");
|
// "/home/iason/Models/TestSet_validPatterns/59.ply");
|
||||||
entry.path();
|
entry.path();
|
||||||
const auto filepathString = filepath.string();
|
const auto filepathString = filepath.string();
|
||||||
std::cout << "Full pattern:" << filepathString << std::endl;
|
|
||||||
// Use only the base triangle version
|
// Use only the base triangle version
|
||||||
const std::string tiledSuffix = "_tiled.ply";
|
const std::string tiledSuffix = "_tiled.ply";
|
||||||
if (filepathString.compare(filepathString.size() - tiledSuffix.size(),
|
if (filepathString.compare(filepathString.size() - tiledSuffix.size(),
|
||||||
|
@ -102,11 +101,11 @@ int main(int argc, char *argv[]) {
|
||||||
}
|
}
|
||||||
FlatPattern pattern(filepathString);
|
FlatPattern pattern(filepathString);
|
||||||
pattern.setLabel(filepath.stem().string());
|
pattern.setLabel(filepath.stem().string());
|
||||||
std::cout << "Testing Pattern:" << filepathString << std::endl;
|
|
||||||
// for (int reducedPatternIndex = 0;
|
// for (int reducedPatternIndex = 0;
|
||||||
// reducedPatternIndex < reducedModels.size();
|
// reducedPatternIndex < reducedModels.size();
|
||||||
// reducedPatternIndex++) {
|
// reducedPatternIndex++) {
|
||||||
// FlatPattern *pReducedModel = reducedModels[reducedPatternIndex];
|
// FlatPattern *pReducedModel =
|
||||||
|
// reducedModels[reducedPatternIndex];
|
||||||
std::unordered_set<size_t> optimizationExcludedEi;
|
std::unordered_set<size_t> optimizationExcludedEi;
|
||||||
// if (pReducedModel !=
|
// if (pReducedModel !=
|
||||||
// reducedModels[0]) { // assumes that the singleBar reduced model
|
// reducedModels[0]) { // assumes that the singleBar reduced model
|
||||||
|
@ -116,9 +115,9 @@ int main(int argc, char *argv[]) {
|
||||||
// optimizationExcludedEi.insert(0);
|
// optimizationExcludedEi.insert(0);
|
||||||
// }
|
// }
|
||||||
std::cout << "Full pattern:" << filepathString << std::endl;
|
std::cout << "Full pattern:" << filepathString << std::endl;
|
||||||
std::cout << "Reduced pattern:" << 1 << std::endl;
|
optimizer.initializePatterns(*reducedModels[2], *reducedModels[0],
|
||||||
optimizer.initialize(pattern, *reducedModels[1], optimizationExcludedEi);
|
optimizationExcludedEi);
|
||||||
Eigen::VectorXd optimalParameters = optimizer.optimize();
|
optimizer.optimize();
|
||||||
// }
|
// }
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -1,251 +0,0 @@
|
||||||
#include "beamformfinder.hpp"
|
|
||||||
#include "edgemesh.hpp"
|
|
||||||
#include "flatpattern.hpp"
|
|
||||||
#include "polyscope/curve_network.h"
|
|
||||||
#include "polyscope/point_cloud.h"
|
|
||||||
#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");
|
|
||||||
|
|
||||||
// Create reduced models
|
|
||||||
// FormFinder::runUnitTests();
|
|
||||||
const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
|
|
||||||
std::vector<vcg::Point2i> singleBarReducedModelEdges{vcg::Point2i(0, 3)};
|
|
||||||
FlatPattern singleBarReducedModel(numberOfNodesPerSlot,
|
|
||||||
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> 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};
|
|
||||||
|
|
||||||
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
|
|
||||||
|
|
||||||
// FlatPattern pattern("/home/iason/Models/TestSet_validPatterns/59.ply");
|
|
||||||
// optimizer.initialize(pattern, *reducedModels[1], {});
|
|
||||||
// optimizer.setInitialGuess
|
|
||||||
// ({0.050000000000000003, 4.0001308096448325, 4.2377893536538567,
|
|
||||||
// 5.6122373437520334});
|
|
||||||
// Eigen::VectorXd optimalParameters = optimizer.optimize(0);
|
|
||||||
|
|
||||||
std::string fullPatternsTestSetDirectory =
|
|
||||||
// "/home/iason/Models/TestSet_validPatterns";
|
|
||||||
"/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
|
|
||||||
"1v_0v_2e_1e_1c_6fan/3/Valid";
|
|
||||||
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/222.ply");
|
|
||||||
// std::filesystem::path(
|
|
||||||
// "/home/iason/Documents/PhD/Research/Pattern_enumerator/Results/"
|
|
||||||
// "1v_0v_2e_1e_1c_6fan/2/Valid/33.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/59.ply");
|
|
||||||
entry.path();
|
|
||||||
const auto filepathString = filepath.string();
|
|
||||||
std::cout << "Full pattern:" << filepathString << std::endl;
|
|
||||||
// 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 (int reducedPatternIndex = 0;
|
|
||||||
// reducedPatternIndex < reducedModels.size();
|
|
||||||
// reducedPatternIndex++) {
|
|
||||||
// FlatPattern *pReducedModel = reducedModels[reducedPatternIndex];
|
|
||||||
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);
|
|
||||||
// }
|
|
||||||
std::cout << "Full pattern:" << filepathString << std::endl;
|
|
||||||
std::cout << "Reduced pattern:" << 1 << std::endl;
|
|
||||||
optimizer.initialize(pattern, *reducedModels[1], optimizationExcludedEi);
|
|
||||||
Eigen::VectorXd optimalParameters = optimizer.optimize();
|
|
||||||
// }
|
|
||||||
}
|
|
||||||
|
|
||||||
// // Full model simulation
|
|
||||||
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
||||||
// fixedVertices;
|
|
||||||
// // fixedVertices[0] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
|
||||||
// fixedVertices[3] = std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
|
|
||||||
// // fixedVertices[7] = std::unordered_set<DoFType>{0, 1, 2};
|
|
||||||
// std::unordered_map<VertexIndex, Vector6d> nodalForces{
|
|
||||||
// {15, {0, 0, 2000, 0, 0, 0}}};
|
|
||||||
// SimulationJob fullModelSimulationJob{
|
|
||||||
// std::make_shared<ElementalMesh>(pattern), fixedVertices, nodalForces,
|
|
||||||
// {}};
|
|
||||||
// Simulator formFinder;
|
|
||||||
// SimulationResults fullModelResults =
|
|
||||||
// formFinder.executeSimulation(fullModelSimulationJob);
|
|
||||||
// fullModelResults.simulationLabel = "Full Model";
|
|
||||||
// fullModelResults.draw(fullModelSimulationJob);
|
|
||||||
// fullModelSimulationJob.draw();
|
|
||||||
// double stiffnessFactor = 1;
|
|
||||||
|
|
||||||
// while (true) {
|
|
||||||
// Reduced model simulation
|
|
||||||
// SimulationJob reducedModelSimulationJob =
|
|
||||||
// optimizer.getReducedSimulationJob(fullModelSimulationJob);
|
|
||||||
// const std::vector<double> stiffnessVector{
|
|
||||||
// fullModelSimulationJob.mesh->elements[0].properties.A,
|
|
||||||
// stiffnessFactor *
|
|
||||||
// fullModelSimulationJob.mesh->elements[0].properties.J, stiffnessFactor
|
|
||||||
// * fullModelSimulationJob.mesh->elements[0].properties.I2,
|
|
||||||
// stiffnessFactor *
|
|
||||||
// fullModelSimulationJob.mesh->elements[0].properties.I3};
|
|
||||||
// for (EdgeIndex ei = 0; ei < reducedModelSimulationJob.mesh->EN(); ei++) {
|
|
||||||
// BeamFormFinder::Element &e =
|
|
||||||
// reducedModelSimulationJob.mesh->elements[ei]; e.properties.A =
|
|
||||||
// 0.00035185827018667374;
|
|
||||||
// // stifnessVector[0];
|
|
||||||
// e.properties.J = // stiffnessVector[1];
|
|
||||||
// 7.709325104874406e-08;
|
|
||||||
// e.properties.I2 = -0.0000015661453308127776;
|
|
||||||
// // stiffnessVector[2];
|
|
||||||
// e.properties.I3 = 3.7099813776947167e-07; // stiffnessVector[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;
|
|
||||||
// }
|
|
||||||
|
|
||||||
// SimulationResults reducedModelResults =
|
|
||||||
// formFinder.executeSimulation(reducedModelSimulationJob, true);
|
|
||||||
// reducedModelResults.simulationLabel =
|
|
||||||
// "Reduced Model_" + std::to_string(stiffnessFactor);
|
|
||||||
// reducedModelResults.draw(reducedModelSimulationJob);
|
|
||||||
// SimulationResultsReporter resultsReporter;
|
|
||||||
// resultsReporter.reportResults(
|
|
||||||
// {reducedModelResults},
|
|
||||||
// std::filesystem::current_path().append("Results"),
|
|
||||||
// "_" + std::to_string(stiffnessFactor));
|
|
||||||
// if (reducedModelResults.history.numberOfSteps ==
|
|
||||||
// BeamFormFinder::Simulator::maxDRMIterations) {
|
|
||||||
// break;
|
|
||||||
// }
|
|
||||||
// stiffnessFactor *= 1.5;
|
|
||||||
// }
|
|
||||||
|
|
||||||
// // Beam
|
|
||||||
// // registerWorldAxes();
|
|
||||||
// VCGEdgeMesh mesh;
|
|
||||||
// // mesh.loadFromPly("/home/iason/Models/simple_beam_model_10elem_1m.ply");
|
|
||||||
// mesh.loadFromPly("/home/iason/Coding/Projects/Approximating shapes with
|
|
||||||
// flat "
|
|
||||||
// "patterns/RodModelOptimizationForPatterns/build/Debug/"
|
|
||||||
// "Fanned_Single bar reduced model.ply");
|
|
||||||
// // vcg::Matrix44d R;
|
|
||||||
// // R.SetRotateDeg(45, {0, 0, 1});
|
|
||||||
// // vcg::tri::UpdatePosition<VCGEdgeMesh>::Matrix(mesh, R);
|
|
||||||
// // mesh.updateEigenEdgeAndVertices();
|
|
||||||
// FormFinder formFinder;
|
|
||||||
// formFinder.runUnitTests();
|
|
||||||
// std::unordered_map<VertexIndex, std::unordered_set<DoFType>>
|
|
||||||
// fixedVertices; fixedVertices[1] = std::unordered_set<DoFType>{2};
|
|
||||||
// fixedVertices[2] = std::unordered_set<DoFType>{2};
|
|
||||||
// fixedVertices[3] = std::unordered_set<DoFType>{1, 2};
|
|
||||||
// fixedVertices[4] = std::unordered_set<DoFType>{2};
|
|
||||||
// fixedVertices[5] = std::unordered_set<DoFType>{2};
|
|
||||||
// fixedVertices[6] = std::unordered_set<DoFType>{0, 1, 2};
|
|
||||||
// // Forced displacements
|
|
||||||
// std::unordered_map<size_t, Eigen::Vector3d> nodalForcedDisplacements;
|
|
||||||
// // nodalForcedDisplacements[10] = Eigen::Vector3d(0, 0.1, 0.1);
|
|
||||||
// std::unordered_map<size_t, VectorType> nodalForcedNormal;
|
|
||||||
|
|
||||||
// CoordType v14 = (mesh.vert[4].cP() - mesh.vert[0].cP()).Normalize();
|
|
||||||
// CoordType v25 = (mesh.vert[5].cP() - mesh.vert[0].cP()).Normalize();
|
|
||||||
// CoordType v36 = (mesh.vert[6].cP() - mesh.vert[0].cP()).Normalize();
|
|
||||||
// const double forceMagnitude = 0.4;
|
|
||||||
// std::unordered_map<VertexIndex, Vector6d> nodalForces{
|
|
||||||
// // {4, Vector6d({0, 0, 0, v14[0], v14[1], 0}) * forceMagnitude},
|
|
||||||
// // {1, Vector6d({0, 0, 0, -v14[0], -v14[1], 0}) *
|
|
||||||
// forceMagnitude},
|
|
||||||
// // {5, Vector6d({0, 0, 0, v25[0], v25[1], 0}) * forceMagnitude},
|
|
||||||
// // {2, Vector6d({0, 0, 0, -v25[0], -v25[1], 0}) *
|
|
||||||
// forceMagnitude},
|
|
||||||
// // {6, Vector6d({0, 0, 0, v36[0], v36[1], 0}) * forceMagnitude},
|
|
||||||
// {3, Vector6d({0, 0, 0, -v36[0], -v36[1], 0}) * forceMagnitude}};
|
|
||||||
// // nodalForcedNormal[10] = v;
|
|
||||||
// // nodalForcedNormal[0] = -v;
|
|
||||||
// // fixedVertices[10] = {0, 1};
|
|
||||||
|
|
||||||
// SimulationJob beamSimulationJob{std::make_shared<SimulationMesh>(mesh),
|
|
||||||
// fixedVertices,
|
|
||||||
// nodalForces,
|
|
||||||
// {},
|
|
||||||
// {}};
|
|
||||||
// beamSimulationJob.mesh->setBeamMaterial(0.3, 1);
|
|
||||||
// beamSimulationJob.mesh->setBeamCrossSection(
|
|
||||||
// RectangularBeamDimensions{0.002, 0.002});
|
|
||||||
// beamSimulationJob.registerForDrawing();
|
|
||||||
// registerWorldAxes();
|
|
||||||
// polyscope::show();
|
|
||||||
// SimulationResults beamSimulationResults =
|
|
||||||
// formFinder.executeSimulation(beamSimulationJob, true, true);
|
|
||||||
// beamSimulationResults.label = "Beam";
|
|
||||||
// beamSimulationResults.registerForDrawing(beamSimulationJob);
|
|
||||||
// polyscope::show();
|
|
||||||
|
|
||||||
return 0;
|
|
||||||
}
|
|
|
@ -8,8 +8,8 @@
|
||||||
const bool gShouldDraw = true;
|
const bool gShouldDraw = true;
|
||||||
|
|
||||||
FormFinder simulator;
|
FormFinder simulator;
|
||||||
Eigen::MatrixX3d g_optimalReducedModelDisplacements;
|
std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
|
||||||
SimulationJob g_reducedPatternSimulationJob;
|
std::vector<SimulationJob> g_reducedPatternSimulationJob;
|
||||||
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
|
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
|
||||||
g_reducedToFullInterfaceViMap;
|
g_reducedToFullInterfaceViMap;
|
||||||
matplot::line_handle gPlotHandle;
|
matplot::line_handle gPlotHandle;
|
||||||
|
@ -18,6 +18,7 @@ Eigen::Vector4d g_initialX;
|
||||||
std::unordered_set<size_t> g_reducedPatternExludedEdges;
|
std::unordered_set<size_t> g_reducedPatternExludedEdges;
|
||||||
// double g_initialParameters;
|
// double g_initialParameters;
|
||||||
Eigen::VectorXd g_initialParameters;
|
Eigen::VectorXd g_initialParameters;
|
||||||
|
ReducedModelOptimizer::SimulationScenario g_chosenSimulationScenarioName;
|
||||||
|
|
||||||
// struct OptimizationCallback {
|
// struct OptimizationCallback {
|
||||||
// double operator()(const size_t &iterations, const Eigen::VectorXd &x,
|
// double operator()(const size_t &iterations, const Eigen::VectorXd &x,
|
||||||
|
@ -96,16 +97,35 @@ Eigen::VectorXd g_initialParameters;
|
||||||
// }
|
// }
|
||||||
//};
|
//};
|
||||||
|
|
||||||
double ReducedModelOptimizer::objective(long n, const double *x) {
|
double ReducedModelOptimizer::computeError(
|
||||||
std::cout.precision(17);
|
const SimulationResults &reducedPatternResults,
|
||||||
|
const Eigen::MatrixX3d &optimalReducedPatternDisplacements) {
|
||||||
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
|
double error = 0;
|
||||||
std::cout << "x[" + std::to_string(parameterIndex) + "]="
|
for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
|
||||||
<< x[parameterIndex] << std::endl;
|
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(
|
||||||
|
reducedPatternResults.displacements[reducedModelVi][0],
|
||||||
|
reducedPatternResults.displacements[reducedModelVi][1],
|
||||||
|
reducedPatternResults.displacements[reducedModelVi][2]);
|
||||||
|
Eigen::Vector3d errorVector =
|
||||||
|
Eigen::Vector3d(
|
||||||
|
optimalReducedPatternDisplacements.row(reducedModelVi)) -
|
||||||
|
vertexDisplacement;
|
||||||
|
// error += errorVector.squaredNorm();
|
||||||
|
error += errorVector.norm();
|
||||||
}
|
}
|
||||||
|
return error;
|
||||||
|
}
|
||||||
|
|
||||||
for (EdgeIndex ei = 0; ei < g_reducedPatternSimulationJob.mesh->EN(); ei++) {
|
void updateMesh(long n, const double *x) {
|
||||||
Element &e = g_reducedPatternSimulationJob.mesh->elements[ei];
|
std::shared_ptr<SimulationMesh> pReducedPatternSimulationMesh =
|
||||||
|
g_reducedPatternSimulationJob[0].mesh;
|
||||||
|
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
|
||||||
|
Element &e = pReducedPatternSimulationMesh->elements[ei];
|
||||||
// if (g_reducedPatternExludedEdges.contains(ei)) {
|
// if (g_reducedPatternExludedEdges.contains(ei)) {
|
||||||
// continue;
|
// continue;
|
||||||
// }
|
// }
|
||||||
|
@ -124,30 +144,56 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
|
||||||
e.secondBendingConstFactor =
|
e.secondBendingConstFactor =
|
||||||
2 * e.properties.E * e.properties.I3 / e.initialLength;
|
2 * e.properties.E * e.properties.I3 / e.initialLength;
|
||||||
}
|
}
|
||||||
// run simulation
|
}
|
||||||
SimulationResults reducedModelResults =
|
|
||||||
simulator.executeSimulation(g_reducedPatternSimulationJob, false, false);
|
double ReducedModelOptimizer::objectiveAllScenarios(long n, const double *x) {
|
||||||
|
std::cout.precision(17);
|
||||||
|
|
||||||
|
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
|
||||||
|
std::cout << "x[" + std::to_string(parameterIndex) + "]="
|
||||||
|
<< x[parameterIndex] << std::endl;
|
||||||
|
}
|
||||||
|
updateMesh(n, x);
|
||||||
|
|
||||||
|
// run simulations
|
||||||
|
double error = 0;
|
||||||
|
for (int simulationScenarioIndex = SimulationScenario::Axial;
|
||||||
|
simulationScenarioIndex !=
|
||||||
|
SimulationScenario::NumberOfSimulationScenarios;
|
||||||
|
simulationScenarioIndex++) {
|
||||||
|
SimulationResults reducedModelResults = simulator.executeSimulation(
|
||||||
|
g_reducedPatternSimulationJob[simulationScenarioIndex], false, false);
|
||||||
|
error += computeError(
|
||||||
|
reducedModelResults,
|
||||||
|
g_optimalReducedModelDisplacements[simulationScenarioIndex]);
|
||||||
|
}
|
||||||
|
|
||||||
// compute error and return it
|
// compute error and return it
|
||||||
double error = 0;
|
gObjectiveValueHistory.push_back(error);
|
||||||
for (const auto reducedFullViPair : g_reducedToFullInterfaceViMap) {
|
auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
|
||||||
VertexIndex reducedModelVi = reducedFullViPair.first;
|
gObjectiveValueHistory.size());
|
||||||
// const auto pos =
|
gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory);
|
||||||
// g_reducedPatternSimulationJob.mesh->vert[reducedModelVi].cP();
|
|
||||||
// std::cout << "Interface vi " << reducedModelVi << " is at position "
|
return error;
|
||||||
// << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
|
}
|
||||||
Eigen::Vector3d vertexDisplacement(
|
|
||||||
reducedModelResults.displacements[reducedModelVi][0],
|
double ReducedModelOptimizer::objectiveSingleScenario(long n, const double *x) {
|
||||||
reducedModelResults.displacements[reducedModelVi][1],
|
std::cout.precision(17);
|
||||||
reducedModelResults.displacements[reducedModelVi][2]);
|
|
||||||
Eigen::Vector3d errorVector =
|
for (size_t parameterIndex = 0; parameterIndex < n; parameterIndex++) {
|
||||||
Eigen::Vector3d(
|
std::cout << "x[" + std::to_string(parameterIndex) + "]="
|
||||||
g_optimalReducedModelDisplacements.row(reducedModelVi)) -
|
<< x[parameterIndex] << std::endl;
|
||||||
vertexDisplacement;
|
|
||||||
// error += errorVector.squaredNorm();
|
|
||||||
error += errorVector.norm();
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
updateMesh(n, x);
|
||||||
|
|
||||||
|
SimulationResults reducedModelResults = simulator.executeSimulation(
|
||||||
|
g_reducedPatternSimulationJob[g_chosenSimulationScenarioName], false,
|
||||||
|
false);
|
||||||
|
const double error = computeError(
|
||||||
|
reducedModelResults,
|
||||||
|
g_optimalReducedModelDisplacements[g_chosenSimulationScenarioName]);
|
||||||
|
|
||||||
gObjectiveValueHistory.push_back(error);
|
gObjectiveValueHistory.push_back(error);
|
||||||
auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
|
auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(),
|
||||||
gObjectiveValueHistory.size());
|
gObjectiveValueHistory.size());
|
||||||
|
@ -302,7 +348,7 @@ ReducedModelOptimizer::ReducedModelOptimizer(
|
||||||
FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
|
FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode);
|
||||||
}
|
}
|
||||||
|
|
||||||
void ReducedModelOptimizer::initialize(
|
void ReducedModelOptimizer::initializePatterns(
|
||||||
FlatPattern &fullPattern, FlatPattern &reducedPattern,
|
FlatPattern &fullPattern, FlatPattern &reducedPattern,
|
||||||
const std::unordered_set<size_t> &reducedModelExcludedEdges) {
|
const std::unordered_set<size_t> &reducedModelExcludedEdges) {
|
||||||
assert(fullPattern.VN() == reducedPattern.VN() &&
|
assert(fullPattern.VN() == reducedPattern.VN() &&
|
||||||
|
@ -316,6 +362,8 @@ void ReducedModelOptimizer::initialize(
|
||||||
computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
|
computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
|
||||||
createSimulationMeshes(copyFullPattern, copyReducedPattern);
|
createSimulationMeshes(copyFullPattern, copyReducedPattern);
|
||||||
initializeStiffnesses();
|
initializeStiffnesses();
|
||||||
|
copyFullPattern.registerForDrawing();
|
||||||
|
polyscope::show();
|
||||||
}
|
}
|
||||||
|
|
||||||
void ReducedModelOptimizer::initializeStiffnesses() {
|
void ReducedModelOptimizer::initializeStiffnesses() {
|
||||||
|
@ -393,7 +441,6 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
|
||||||
const SimulationResults &fullModelResults,
|
const SimulationResults &fullModelResults,
|
||||||
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {
|
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) {
|
||||||
|
|
||||||
optimalDisplacementsOfReducedModel.resize(0, 3);
|
|
||||||
optimalDisplacementsOfReducedModel.resize(
|
optimalDisplacementsOfReducedModel.resize(
|
||||||
m_pReducedPatternSimulationMesh->VN(), 3);
|
m_pReducedPatternSimulationMesh->VN(), 3);
|
||||||
optimalDisplacementsOfReducedModel.setZero(
|
optimalDisplacementsOfReducedModel.setZero(
|
||||||
|
@ -411,28 +458,11 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
|
Eigen::VectorXd ReducedModelOptimizer::runOptimization(
|
||||||
const SimulationJob &fullPatternSimulationJob) {
|
double (*pObjectiveFunction)(long, const double *)) {
|
||||||
|
|
||||||
// fullPatternSimulationJob.registerForDrawing();
|
|
||||||
// polyscope::show();
|
|
||||||
gObjectiveValueHistory.clear();
|
gObjectiveValueHistory.clear();
|
||||||
fullPatternSimulationJob.mesh->savePly(
|
|
||||||
"Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply");
|
|
||||||
SimulationResults fullModelResults =
|
|
||||||
simulator.executeSimulation(fullPatternSimulationJob, false);
|
|
||||||
fullModelResults.label = "fullModel";
|
|
||||||
fullModelResults.registerForDrawing(fullPatternSimulationJob);
|
|
||||||
registerWorldAxes();
|
|
||||||
polyscope::show();
|
|
||||||
computeDesiredReducedModelDisplacements(fullModelResults,
|
|
||||||
g_optimalReducedModelDisplacements);
|
|
||||||
computeReducedModelSimulationJob(fullPatternSimulationJob,
|
|
||||||
g_reducedPatternSimulationJob);
|
|
||||||
// gReducedPatternSimulationJob.registerForDrawing();
|
|
||||||
// polyscope::show();
|
|
||||||
|
|
||||||
double (*pObjectiveFunction)(long, const double *) = &objective;
|
|
||||||
// const size_t n = m_pReducedPatternSimulationMesh->EN();
|
// const size_t n = m_pReducedPatternSimulationMesh->EN();
|
||||||
const size_t n = 4;
|
const size_t n = 4;
|
||||||
const size_t npt = 2 * n;
|
const size_t npt = 2 * n;
|
||||||
|
@ -442,7 +472,7 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
|
||||||
"conditions is not recommended.");
|
"conditions is not recommended.");
|
||||||
// Set initial guess of solution
|
// Set initial guess of solution
|
||||||
|
|
||||||
std::vector<double> x(n, 4);
|
std::vector<double> x(n, 2);
|
||||||
if (!initialGuess.empty()) {
|
if (!initialGuess.empty()) {
|
||||||
x = initialGuess;
|
x = initialGuess;
|
||||||
} // {0.10000000000000 001, 2, 1.9999999971613847, 6.9560343643347764};
|
} // {0.10000000000000 001, 2, 1.9999999971613847, 6.9560343643347764};
|
||||||
|
@ -450,8 +480,8 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
|
||||||
// {0.0001, 2, 2.000000005047502, 1.3055270196964464};
|
// {0.0001, 2, 2.000000005047502, 1.3055270196964464};
|
||||||
// {initialGuess(0), initialGuess(1), initialGuess(2),
|
// {initialGuess(0), initialGuess(1), initialGuess(2),
|
||||||
// initialGuess(3)};
|
// initialGuess(3)};
|
||||||
const double xMin = 0.12;
|
const double xMin = 1e-2;
|
||||||
const double xMax = 10;
|
const double xMax = 5000;
|
||||||
// assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
|
// assert(x.end() == find_if(x.begin(), x.end(), [&](const double &d) {
|
||||||
// return d >= xMax || d <= xMin;
|
// return d >= xMax || d <= xMin;
|
||||||
// }));
|
// }));
|
||||||
|
@ -470,29 +500,6 @@ Eigen::VectorXd ReducedModelOptimizer::optimizeForSimulationJob(
|
||||||
bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
|
bobyqa(pObjectiveFunction, n, npt, x.data(), xLow.data(), xUpper.data(),
|
||||||
rhobeg, rhoend, maxFun, w.data());
|
rhobeg, rhoend, maxFun, w.data());
|
||||||
|
|
||||||
SimulationResults reducedModelOptimizedResults =
|
|
||||||
simulator.executeSimulation(g_reducedPatternSimulationJob);
|
|
||||||
|
|
||||||
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.label = "reducedModel";
|
|
||||||
reducedModelOptimizedResults.registerForDrawing(
|
|
||||||
g_reducedPatternSimulationJob);
|
|
||||||
polyscope::show();
|
|
||||||
|
|
||||||
Eigen::VectorXd eigenX(x.size(), 1);
|
Eigen::VectorXd eigenX(x.size(), 1);
|
||||||
for (size_t xi = 0; xi < x.size(); xi++) {
|
for (size_t xi = 0; xi < x.size(); xi++) {
|
||||||
eigenX(xi) = x[xi];
|
eigenX(xi) = x[xi];
|
||||||
|
@ -659,29 +666,77 @@ std::vector<SimulationJob> ReducedModelOptimizer::createScenarios(
|
||||||
return scenarios;
|
return scenarios;
|
||||||
}
|
}
|
||||||
|
|
||||||
Eigen::VectorXd ReducedModelOptimizer::optimize(const int &simulationScenario) {
|
void ReducedModelOptimizer::optimize(const int &simulationScenario) {
|
||||||
|
|
||||||
std::vector<SimulationJob> simulationJobs =
|
std::vector<SimulationJob> simulationJobs =
|
||||||
createScenarios(m_pFullModelSimulationMesh);
|
createScenarios(m_pFullModelSimulationMesh);
|
||||||
std::vector<Eigen::VectorXd> results;
|
g_optimalReducedModelDisplacements.resize(6);
|
||||||
if (simulationScenario == -1) { // run all scenarios
|
g_reducedPatternSimulationJob.resize(6);
|
||||||
for (int i = 0; i < simulationJobs.size(); i++) {
|
polyscope::removeAllStructures();
|
||||||
std::cout << "Running simulation job:" << i << std::endl;
|
|
||||||
const SimulationJob &job = simulationJobs[i];
|
for (int simulationScenarioIndex = SimulationScenario::Axial;
|
||||||
polyscope::removeAllStructures();
|
simulationScenarioIndex !=
|
||||||
auto result = optimizeForSimulationJob(job);
|
SimulationScenario::NumberOfSimulationScenarios;
|
||||||
results.push_back(result);
|
simulationScenarioIndex++) {
|
||||||
}
|
const SimulationJob &fullPatternSimulationJob =
|
||||||
|
simulationJobs[simulationScenarioIndex];
|
||||||
|
// fullPatternSimulationJob.mesh->savePly(
|
||||||
|
// "Fanned_" + m_pFullModelSimulationMesh->getLabel() + ".ply");
|
||||||
|
SimulationResults fullModelResults =
|
||||||
|
simulator.executeSimulation(fullPatternSimulationJob, false);
|
||||||
|
// fullModelResults.label =
|
||||||
|
// "fullModel_" + SimulationScenarioStrings[simulationScenarioIndex];
|
||||||
|
// fullModelResults.registerForDrawing(fullPatternSimulationJob);
|
||||||
|
computeDesiredReducedModelDisplacements(
|
||||||
|
fullModelResults,
|
||||||
|
g_optimalReducedModelDisplacements[simulationScenarioIndex]);
|
||||||
|
computeReducedModelSimulationJob(
|
||||||
|
fullPatternSimulationJob,
|
||||||
|
g_reducedPatternSimulationJob[simulationScenarioIndex]);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (simulationScenario == -1) {
|
||||||
|
double (*pObjectiveFunction)(long, const double *) = &objectiveAllScenarios;
|
||||||
|
runOptimization(pObjectiveFunction);
|
||||||
} else { // run chosen
|
} else { // run chosen
|
||||||
std::cout << "Running simulation job:" << simulationScenario << std::endl;
|
g_chosenSimulationScenarioName = SimulationScenario(simulationScenario);
|
||||||
const SimulationJob &job = simulationJobs[simulationScenario];
|
double (*pObjectiveFunction)(long, const double *) =
|
||||||
polyscope::removeAllStructures();
|
&objectiveSingleScenario;
|
||||||
auto result = optimizeForSimulationJob(job);
|
runOptimization(pObjectiveFunction);
|
||||||
results.push_back(result);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
if (results.empty()) {
|
m_pFullModelSimulationMesh->registerForDrawing();
|
||||||
return Eigen::VectorXd();
|
double error = 0;
|
||||||
|
for (int simulationScenarioIndex = SimulationScenario::Axial;
|
||||||
|
simulationScenarioIndex !=
|
||||||
|
SimulationScenario::NumberOfSimulationScenarios;
|
||||||
|
simulationScenarioIndex++) {
|
||||||
|
const SimulationJob &fullPatternSimulationJob =
|
||||||
|
simulationJobs[simulationScenarioIndex];
|
||||||
|
SimulationResults fullModelResults =
|
||||||
|
simulator.executeSimulation(fullPatternSimulationJob, false);
|
||||||
|
fullModelResults.label =
|
||||||
|
"fullModel_" + SimulationScenarioStrings[simulationScenarioIndex];
|
||||||
|
fullModelResults.registerForDrawing(fullPatternSimulationJob);
|
||||||
|
const SimulationJob &reducedPatternSimulationJob =
|
||||||
|
g_reducedPatternSimulationJob[simulationScenarioIndex];
|
||||||
|
SimulationResults reducedModelResults =
|
||||||
|
simulator.executeSimulation(reducedPatternSimulationJob, false, false);
|
||||||
|
error += computeError(
|
||||||
|
reducedModelResults,
|
||||||
|
g_optimalReducedModelDisplacements[simulationScenarioIndex]);
|
||||||
|
reducedModelResults.label =
|
||||||
|
"reducedModel_" + SimulationScenarioStrings[simulationScenarioIndex];
|
||||||
|
reducedModelResults.registerForDrawing(reducedPatternSimulationJob);
|
||||||
|
// registerWorldAxes();
|
||||||
|
std::cout << "A full:"
|
||||||
|
<< m_pFullModelSimulationMesh->elements[0].properties.A
|
||||||
|
<< std::endl;
|
||||||
|
std::cout << "A reduced:"
|
||||||
|
<< m_pReducedPatternSimulationMesh->elements[0].properties.A
|
||||||
|
<< std::endl;
|
||||||
|
polyscope::show();
|
||||||
|
fullModelResults.unregister();
|
||||||
|
reducedModelResults.unregister();
|
||||||
}
|
}
|
||||||
|
|
||||||
return results[0];
|
|
||||||
}
|
}
|
||||||
|
|
|
@ -22,7 +22,17 @@ class ReducedModelOptimizer {
|
||||||
std::vector<double> initialGuess;
|
std::vector<double> initialGuess;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
Eigen::VectorXd optimize(const int &simulationScenario = -1);
|
enum SimulationScenario {
|
||||||
|
Axial,
|
||||||
|
Shear,
|
||||||
|
Bending,
|
||||||
|
Double,
|
||||||
|
Saddle,
|
||||||
|
NumberOfSimulationScenarios
|
||||||
|
};
|
||||||
|
inline static const std::string SimulationScenarioStrings[] = {
|
||||||
|
"Axial", "Shear", "Bending", "Double", "Saddle"};
|
||||||
|
void optimize(const int &simulationScenario = -1);
|
||||||
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
|
double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const;
|
||||||
|
|
||||||
ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
|
ReducedModelOptimizer(const std::vector<size_t> &numberOfNodesPerSlot);
|
||||||
|
@ -33,8 +43,9 @@ public:
|
||||||
SimulationJob
|
SimulationJob
|
||||||
getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
|
getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
|
||||||
|
|
||||||
void initialize(FlatPattern &fullPattern, FlatPattern &reducedPatterm,
|
void initializePatterns(
|
||||||
const std::unordered_set<size_t> &reducedModelExcludedEges);
|
FlatPattern &fullPattern, FlatPattern &reducedPatterm,
|
||||||
|
const std::unordered_set<size_t> &reducedModelExcludedEges);
|
||||||
|
|
||||||
void setInitialGuess(std::vector<double> v);
|
void setInitialGuess(std::vector<double> v);
|
||||||
|
|
||||||
|
@ -42,8 +53,8 @@ private:
|
||||||
void computeDesiredReducedModelDisplacements(
|
void computeDesiredReducedModelDisplacements(
|
||||||
const SimulationResults &fullModelResults,
|
const SimulationResults &fullModelResults,
|
||||||
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
|
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
|
||||||
Eigen::VectorXd
|
Eigen::VectorXd runOptimization(double (*pObjectiveFunction)(long,
|
||||||
optimizeForSimulationJob(const SimulationJob &fullModelSimulationJob);
|
const double *));
|
||||||
std::vector<SimulationJob>
|
std::vector<SimulationJob>
|
||||||
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
|
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
|
||||||
void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,
|
void computeMaps(FlatPattern &fullModel, FlatPattern &reducedPattern,
|
||||||
|
@ -52,7 +63,11 @@ private:
|
||||||
FlatPattern &reducedModel);
|
FlatPattern &reducedModel);
|
||||||
void initializeStiffnesses();
|
void initializeStiffnesses();
|
||||||
|
|
||||||
static double objective(long n, const double *x);
|
static double
|
||||||
|
computeError(const SimulationResults &reducedPatternResults,
|
||||||
|
const Eigen::MatrixX3d &optimalReducedPatternDisplacements);
|
||||||
|
static double objectiveSingleScenario(long n, const double *x);
|
||||||
|
static double objectiveAllScenarios(long n, const double *x);
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif // REDUCEDMODELOPTIMIZER_HPP
|
#endif // REDUCEDMODELOPTIMIZER_HPP
|
||||||
|
|
Loading…
Reference in New Issue