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

Refactoring

This commit is contained in:
iasonmanolas 2021-03-23 19:20:46 +02:00
parent 4c0c5307b9
commit 4daed81fef
6 changed files with 137 additions and 652 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
CMakeLists.txt
View File

@ -68,15 +68,6 @@ download_project(PROJ vcglib_devel
) )
file(GLOB EXT_SOURCES ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp) file(GLOB EXT_SOURCES ${vcglib_devel_SOURCE_DIR}/wrap/ply/plylib.cpp)
##threed-beam-fea
download_project(PROJ threed-beam-fea
GIT_REPOSITORY https://github.com/IasonManolas/threed-beam-fea.git
GIT_TAG master
PREFIX ${EXTERNAL_DEPS_DIR}
${UPDATE_DISCONNECTED_IF_AVAILABLE}
)
add_subdirectory(${threed-beam-fea_SOURCE_DIR} ${threed-beam-fea_BINARY_DIR})
##Eigen 3 NOTE: Eigen is required on the system the code is ran ##Eigen 3 NOTE: Eigen is required on the system the code is ran
find_package(Eigen3 3.3 REQUIRED) find_package(Eigen3 3.3 REQUIRED)
@ -94,7 +85,6 @@ target_compile_features(${PROJECT_NAME} PUBLIC cxx_std_20)
target_include_directories(${PROJECT_NAME} target_include_directories(${PROJECT_NAME}
PUBLIC ${vcglib_devel_SOURCE_DIR} PUBLIC ${vcglib_devel_SOURCE_DIR}
PUBLIC ${threed-beam-fea_SOURCE_DIR}/include
PUBLIC ${MYSOURCES_SOURCE_DIR} PUBLIC ${MYSOURCES_SOURCE_DIR}
PUBLIC ${MYSOURCES_SOURCE_DIR}/boost_graph PUBLIC ${MYSOURCES_SOURCE_DIR}/boost_graph
) )
@ -102,7 +92,7 @@ target_include_directories(${PROJECT_NAME}
target_link_directories(${PROJECT_NAME} PRIVATE ${MYSOURCES_SOURCE_DIR}/boost_graph/libs/) target_link_directories(${PROJECT_NAME} PRIVATE ${MYSOURCES_SOURCE_DIR}/boost_graph/libs/)
if(${USE_POLYSCOPE}) if(${USE_POLYSCOPE})
target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen dlib::dlib ThreedBeamFEA MySources) target_link_libraries(${PROJECT_NAME} polyscope Eigen3::Eigen dlib::dlib MySources)
else() else()
target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen dlib::dlib ThreedBeamFEA MySources) target_link_libraries(${PROJECT_NAME} -static Eigen3::Eigen dlib::dlib MySources)
endif() endif()

0
src/linearsimulationmodel.cpp
View File

215
src/linearsimulationmodel.hpp
View File

@ -1,215 +0,0 @@
#ifndef LINEARSIMULATIONMODEL_HPP
#define LINEARSIMULATIONMODEL_HPP
//#include "beam.hpp"
#include "simulationresult.hpp"
#include "threed_beam_fea.h"
#include <filesystem>
#include <vector>
// struct BeamSimulationProperties {
// float crossArea;
// float I2;
// float I3;
// float polarInertia;
// float G;
// // Properties used by fea
// float EA;
// float EIz;
// float EIy;
// float GJ;
// BeamSimulationProperties(const BeamDimensions &dimensions,
// const BeamMaterial &material);
//};
// struct NodalForce {
// int index;
// int dof;
// double magnitude;
//};
// struct SimulationJob {
// Eigen::MatrixX3d nodes;
// Eigen::MatrixX2i elements;
// Eigen::MatrixX3d elementalNormals;
// Eigen::VectorXi fixedNodes;
// std::vector<NodalForce> nodalForces;
// std::vector<BeamDimensions> beamDimensions;
// std::vector<BeamMaterial> beamMaterial;
//};
// struct SimulationResults {
// std::vector<Eigen::VectorXd> edgeForces; ///< Force values per force
// component
// ///< #force components x #edges
// Eigen::MatrixXd
// nodalDisplacements; ///< The displacement of each node #nodes x 3
// SimulationResults(const fea::Summary &feaSummary);
// SimulationResults() {}
//};
class LinearSimulationModel {
public:
LinearSimulationModel(){
}
static std::vector<fea::Elem>
getFeaElements(const std::shared_ptr<SimulationJob> &job) {
const int numberOfEdges = job->pMesh->EN();
std::vector<fea::Elem> elements(numberOfEdges);
for (int edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
const SimulationMesh::CoordType &evn0 =
job->pMesh->edge[edgeIndex].cV(0)->cN();
const SimulationMesh::CoordType &evn1 =
job->pMesh->edge[edgeIndex].cV(1)->cN();
const std::vector<double> nAverageVector{(evn0[0] + evn1[0]) / 2,
(evn0[1] + evn1[1]) / 2,
(evn0[2] + evn1[2]) / 2};
const Element &element = job->pMesh->elements[edgeIndex];
const double E = element.material.youngsModulus;
fea::Props feaProperties(E * element.A, E * element.I3, E * element.I2,
element.G * element.J, nAverageVector);
const int vi0 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(0));
const int vi1 = job->pMesh->getIndex(job->pMesh->edge[edgeIndex].cV(1));
elements[edgeIndex] = fea::Elem(vi0, vi1, feaProperties);
}
return elements;
}
static std::vector<fea::Node>
getFeaNodes(const std::shared_ptr<SimulationJob> &job) {
const int numberOfNodes = job->pMesh->VN();
std::vector<fea::Node> feaNodes(numberOfNodes);
for (int vi = 0; vi < numberOfNodes; vi++) {
const CoordType &p = job->pMesh->vert[vi].cP();
feaNodes[vi] = fea::Node(p[0], p[1], p[2]);
}
return feaNodes;
}
static std::vector<fea::BC>
getFeaFixedNodes(const std::shared_ptr<SimulationJob> &job) {
std::vector<fea::BC> boundaryConditions;
boundaryConditions.reserve(job->constrainedVertices.size() * 6);
for (auto fixedVertex : job->constrainedVertices) {
const int vertexIndex = fixedVertex.first;
for (int dofIndex : fixedVertex.second) {
boundaryConditions.emplace_back(
fea::BC(vertexIndex, static_cast<fea::DOF>(dofIndex), 0));
}
}
return boundaryConditions;
}
static std::vector<fea::Force>
getFeaNodalForces(const std::shared_ptr<SimulationJob> &job) {
std::vector<fea::Force> nodalForces;
nodalForces.reserve(job->nodalExternalForces.size() * 6);
for (auto nodalForce : job->nodalExternalForces) {
for (int dofIndex = 0; dofIndex < 6; dofIndex++) {
if (nodalForce.second[dofIndex] == 0) {
continue;
}
nodalForces.emplace_back(
fea::Force(nodalForce.first, dofIndex, nodalForce.second[dofIndex]));
}
}
return nodalForces;
}
static SimulationResults getResults(const fea::Summary &feaSummary) {
SimulationResults results;
results.executionTime = feaSummary.total_time_in_ms * 1000;
// displacements
results.displacements.resize(feaSummary.num_nodes);
for (int vi = 0; vi < feaSummary.num_nodes; vi++) {
results.displacements[vi] = Vector6d(feaSummary.nodal_displacements[vi]);
}
// // Convert forces
// // Convert to vector of eigen matrices of the form force component-> per
// // Edge
// const int numDof = 6;
// const size_t numberOfEdges = feaSummary.element_forces.size();
// edgeForces =
// std::vector<Eigen::VectorXd>(numDof, Eigen::VectorXd(2 *
// numberOfEdges));
// for (gsl::index edgeIndex = 0; edgeIndex < numberOfEdges; edgeIndex++) {
// for (gsl::index forceComponentIndex = 0; forceComponentIndex < numDof;
// forceComponentIndex++) {
// (edgeForces[forceComponentIndex])(2 * edgeIndex) =
// feaSummary.element_forces[edgeIndex][forceComponentIndex];
// (edgeForces[forceComponentIndex])(2 * edgeIndex + 1) =
// feaSummary.element_forces[edgeIndex][numDof +
// forceComponentIndex];
// }
// }
return results;
}
SimulationResults
executeSimulation(const std::shared_ptr<SimulationJob> &simulationJob) {
assert(simulationJob->pMesh->VN() != 0);
fea::Job job(getFeaNodes(simulationJob), getFeaElements(simulationJob));
// printInfo(job);
// create the default options
fea::Options opts;
opts.save_elemental_forces = false;
opts.save_nodal_displacements = false;
opts.save_nodal_forces = false;
opts.save_report = false;
opts.save_tie_forces = false;
// if (!elementalForcesOutputFilepath.empty()) {
// opts.save_elemental_forces = true;
// opts.elemental_forces_filename = elementalForcesOutputFilepath;
// }
// if (!nodalDisplacementsOutputFilepath.empty()) {
// opts.save_nodal_displacements = true;
// opts.nodal_displacements_filename = nodalDisplacementsOutputFilepath;
// }
// have the program output status updates
opts.verbose = false;
// form an empty vector of ties
std::vector<fea::Tie> ties;
// also create an empty list of equations
std::vector<fea::Equation> equations;
auto fixedVertices = getFeaFixedNodes(simulationJob);
auto nodalForces = getFeaNodalForces(simulationJob);
fea::Summary feaResults =
fea::solve(job, fixedVertices, nodalForces, ties, equations, opts);
SimulationResults results = getResults(feaResults);
results.job = simulationJob;
return results;
}
// SimulationResults getResults() const;
// void setResultsNodalDisplacementCSVFilepath(const std::string
// &outputPath); void setResultsElementalForcesCSVFilepath(const std::string
// &outputPath);
private:
// std::string nodalDisplacementsOutputFilepath{"nodal_displacement.csv"};
// std::string elementalForcesOutputFilepath{"elemental_forces.csv"};
// SimulationResults results;
static void printInfo(const fea::Job &job) {
std::cout << "Details regarding the fea::Job:" << std::endl;
std::cout << "Nodes:" << std::endl;
for (fea::Node n : job.nodes) {
std::cout << n << std::endl;
}
std::cout << "Elements:" << std::endl;
for (Eigen::Vector2i e : job.elems) {
std::cout << e << std::endl;
}
}
};
#endif // LINEARSIMULATIONMODEL_HPP

24
src/main.cpp
View File

@ -1,4 +1,4 @@
#include "beamformfinder.hpp" #include "drmsimulationmodel.hpp"
#include "csvfile.hpp" #include "csvfile.hpp"
#include "edgemesh.hpp" #include "edgemesh.hpp"
#include "reducedmodeloptimizer.hpp" #include "reducedmodeloptimizer.hpp"
@ -41,21 +41,21 @@ int main(int argc, char *argv[]) {
reducedPattern.scale(0.03,interfaceNodeIndex); reducedPattern.scale(0.03,interfaceNodeIndex);
// Set the optization settings // Set the optization settings
ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9}; ReducedModelOptimization::xRange beamE{"E", 0.001, 10000};
ReducedModelOptimizer::xRange beamA{"A", 0.25, 2.25}; ReducedModelOptimization::xRange beamA{"A", 0.001, 10000};
ReducedModelOptimizer::xRange beamI2{"I2", std::sqrt(beamA.min)*std::pow(std::sqrt(beamA.min),3)/12, std::sqrt(beamA.max)*std::pow(std::sqrt(beamA.max),3)/12}; ReducedModelOptimization::xRange beamI2{"I2", 0.001,10000};
ReducedModelOptimizer::xRange beamI3{"I3", std::sqrt(beamA.min)*std::pow(std::sqrt(beamA.min),3)/12, std::sqrt(beamA.max)*std::pow(std::sqrt(beamA.max),3)/12}; ReducedModelOptimization::xRange beamI3{"I3", 0.001,10000};
ReducedModelOptimizer::xRange beamJ{"J", beamI2.min+beamI3.min, beamI2.max+beamI3.max}; ReducedModelOptimization::xRange beamJ{"J", 0.001,10000};
ReducedModelOptimizer::xRange innerHexagonSize{"HexSize", 0.1, 0.8}; ReducedModelOptimization::xRange innerHexagonSize{"HexSize", 0.1, 0.8};
ReducedModelOptimizer::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5}; ReducedModelOptimization::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5};
ReducedModelOptimizer::Settings settings_optimization; ReducedModelOptimization::Settings settings_optimization;
settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3, settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3,
innerHexagonSize, innerHexagonAngle}; innerHexagonSize, innerHexagonAngle};
const bool input_numberOfFunctionCallsDefined = argc >= 4; const bool input_numberOfFunctionCallsDefined = argc >= 4;
settings_optimization.numberOfFunctionCalls = settings_optimization.numberOfFunctionCalls =
input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100;
settings_optimization.normalizationStrategy = settings_optimization.normalizationStrategy =
ReducedModelOptimizer::Settings::NormalizationStrategy::Epsilon; ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon;
settings_optimization.normalizationParameter = 0.0003; settings_optimization.normalizationParameter = 0.0003;
settings_optimization.solutionAccuracy = 0.01; settings_optimization.solutionAccuracy = 0.01;
@ -66,8 +66,8 @@ int main(int argc, char *argv[]) {
const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]); assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]);
ReducedModelOptimizer optimizer(numberOfNodesPerSlot); ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
optimizer.initializePatterns(fullPattern, reducedPattern, {}); optimizer.initializePatterns(fullPattern, reducedPattern, {},settings_optimization.xRanges.size());
ReducedModelOptimizer::Results optimizationResults = ReducedModelOptimization::Results optimizationResults =
optimizer.optimize(settings_optimization); optimizer.optimize(settings_optimization);
// Export results // Export results

205
src/reducedmodeloptimizer.cpp
View File

@ -7,6 +7,8 @@
#include <dlib/global_optimization.h> #include <dlib/global_optimization.h>
#include <dlib/optimization.h> #include <dlib/optimization.h>
using namespace ReducedModelOptimization;
struct GlobalOptimizationVariables { struct GlobalOptimizationVariables {
std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements; std::vector<Eigen::MatrixX3d> g_optimalReducedModelDisplacements;
std::vector<std::vector<Vector6d>> fullPatternDisplacements; std::vector<std::vector<Vector6d>> fullPatternDisplacements;
@ -16,9 +18,9 @@ struct GlobalOptimizationVariables {
std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
reducedToFullInterfaceViMap; reducedToFullInterfaceViMap;
matplot::line_handle gPlotHandle; matplot::line_handle gPlotHandle;
std::vector<double> gObjectiveValueHistory; std::vector<double> objectiveValueHistory;
Eigen::VectorXd initialParameters; Eigen::VectorXd initialParameters;
std::vector<ReducedModelOptimizer::SimulationScenario> std::vector<SimulationScenario>
simulationScenarioIndices; simulationScenarioIndices;
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)}; std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
double innerHexagonInitialRotationAngle{30}; double innerHexagonInitialRotationAngle{30};
@ -28,7 +30,7 @@ struct GlobalOptimizationVariables {
int numOfSimulationCrashes{false}; int numOfSimulationCrashes{false};
int numberOfFunctionCalls{0}; int numberOfFunctionCalls{0};
int numberOfOptimizationParameters{5}; int numberOfOptimizationParameters{5};
ReducedModelOptimizer::Settings optimizationSettings; ReducedModelOptimization::Settings optimizationSettings;
} global; } global;
std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs; std::vector<SimulationJob> reducedPatternMaximumDisplacementSimulationJobs;
@ -81,47 +83,6 @@ double ReducedModelOptimizer::computeRawError(
return error; return error;
} }
void updateMesh(long n, const double *x) {
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
->pMesh;
const double E=global.initialParameters(0)*x[0];
const double A=global.initialParameters(1) * x[1];
const double beamWidth=std::sqrt(A);
const double beamHeight=beamWidth;
const double J=global.initialParameters(2) * x[2];
const double I2=global.initialParameters(3) * x[3];
const double I3=global.initialParameters(4) * x[4];
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(
RectangularBeamDimensions(beamWidth,
beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio,
E));
e.A = A;
e.J = J;
e.I2 = I2;
e.I3 = I3;
}
assert(pReducedPatternSimulationMesh->EN() == 12);
auto R = vcg::RotationMatrix(
ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(x[6] - global.innerHexagonInitialRotationAngle));
for (int rotationCounter = 0;
rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() =
R * global.g_innerHexagonVectors[rotationCounter] * x[5];
}
pReducedPatternSimulationMesh->reset();
#ifdef POLYSCOPE_DEFINED
pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
#endif
}
double ReducedModelOptimizer::objective(double b, double r, double E) { double ReducedModelOptimizer::objective(double b, double r, double E) {
std::vector<double> x{b, r, E}; std::vector<double> x{b, r, E};
return ReducedModelOptimizer::objective(x.size(), x.data()); return ReducedModelOptimizer::objective(x.size(), x.data());
@ -134,8 +95,19 @@ double ReducedModelOptimizer::objective(double E,double A,double J,double I2,dou
return ReducedModelOptimizer::objective(x.size(), x.data()); return ReducedModelOptimizer::objective(x.size(), x.data());
} }
double ReducedModelOptimizer::objective(double E,double A,
double innerHexagonSize,
double innerHexagonRotationAngle) {
std::vector<double> x{E,A, innerHexagonSize, innerHexagonRotationAngle};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
double ReducedModelOptimizer::objective(long n, const double *x) { double ReducedModelOptimizer::objective(long n, const double *x) {
// std::cout.precision(17); // std::cout.precision(17);
// for(int i=0;i<n;i++){
// std::cout<<x[i]<<" ";
// }
// std::cout<<std::endl;
// const Element &e = // const Element &e =
// global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout <<
@ -219,19 +191,18 @@ double ReducedModelOptimizer::objective(long n, const double *x) {
<< std::endl; << std::endl;
} }
#endif #endif
// compute error and return it // compute error and return it
// global.gObjectiveValueHistory.push_back(error); // global.objectiveValueHistory.push_back(totalError);
// auto xPlot = matplot::linspace(0, gObjectiveValueHistory.size(), // auto xPlot = matplot::linspace(0, global.objectiveValueHistory.size(),
// gObjectiveValueHistory.size()); // global.objectiveValueHistory.size());
// std::vector<double> colors(gObjectiveValueHistory.size(), 2); // std::vector<double> colors(global.gObjectiveValueHistory.size(), 2);
// if (g_firstRoundIterationIndex != 0) { // if (global.g_firstRoundIterationIndex != 0) {
// for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(), // for_each(colors.begin() + g_firstRoundIterationIndex, colors.end(),
// [](double &c) { c = 0.7; }); // [](double &c) { c = 0.7; });
// } // }
// gPlotHandle = matplot::scatter(xPlot, gObjectiveValueHistory, 6, colors); // global.gPlotHandle = matplot::scatter(xPlot, global.objectiveValueHistory);
// SimulationResultsReporter::createPlot("Number of Steps", "Objective // SimulationResultsReporter::createPlot("Number of Steps", "Objective value",
// value", // global.objectiveValueHistory);
// gObjectiveValueHistory);
return totalError; return totalError;
} }
@ -416,7 +387,7 @@ ReducedModelOptimizer::ReducedModelOptimizer(
void ReducedModelOptimizer::initializePatterns( void ReducedModelOptimizer::initializePatterns(
PatternGeometry &fullPattern, PatternGeometry &reducedPattern, PatternGeometry &fullPattern, PatternGeometry &reducedPattern,
const std::unordered_set<size_t> &reducedModelExcludedEdges) { const std::unordered_set<size_t> &reducedModelExcludedEdges,const int& optimizationParameters) {
// fullPattern.setLabel("full_pattern_" + fullPattern.getLabel()); // fullPattern.setLabel("full_pattern_" + fullPattern.getLabel());
// reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel()); // reducedPattern.setLabel("reduced_pattern_" + reducedPattern.getLabel());
assert(fullPattern.VN() == reducedPattern.VN() && assert(fullPattern.VN() == reducedPattern.VN() &&
@ -455,31 +426,64 @@ void ReducedModelOptimizer::initializePatterns(
This is not very generic */ This is not very generic */
computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges); computeMaps(copyFullPattern, copyReducedPattern, reducedModelExcludedEdges);
createSimulationMeshes(copyFullPattern, copyReducedPattern); createSimulationMeshes(copyFullPattern, copyReducedPattern);
initializeOptimizationParameters(m_pReducedPatternSimulationMesh); initializeOptimizationParameters(m_pReducedPatternSimulationMesh,optimizationParameters);
} }
void updateMesh(long n, const double *x) {
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
->pMesh;
const double E=global.initialParameters(0)*x[0];
const double A=global.initialParameters(1) * x[1];
const double beamWidth=std::sqrt(A);
const double beamHeight=beamWidth;
const double J=global.initialParameters(2) * x[2];
const double I2=global.initialParameters(3) * x[3];
const double I3=global.initialParameters(4) * x[4];
for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
Element &e = pReducedPatternSimulationMesh->elements[ei];
e.setDimensions(
RectangularBeamDimensions(beamWidth,
beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio,
E));
e.J = J;
e.I2 = I2;
e.I3 = I3;
}
assert(pReducedPatternSimulationMesh->EN() == 12);
assert(n>=2);
auto R = vcg::RotationMatrix(
ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(x[n-1] - global.innerHexagonInitialRotationAngle));
for (int rotationCounter = 0;
rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) {
pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() =
R * global.g_innerHexagonVectors[rotationCounter] * x[n-2];
}
pReducedPatternSimulationMesh->reset();
#ifdef POLYSCOPE_DEFINED
pReducedPatternSimulationMesh->updateEigenEdgeAndVertices();
#endif
}
void ReducedModelOptimizer::initializeOptimizationParameters( void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh) { const std::shared_ptr<SimulationMesh> &mesh,const int& optimizationParamters) {
global.numberOfOptimizationParameters = 7; global.numberOfOptimizationParameters = optimizationParamters;
global.initialParameters.resize(global.numberOfOptimizationParameters); global.initialParameters.resize(global.numberOfOptimizationParameters);
// 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 = 5;
// e.axialConstFactor *= stiffnessFactor;
// e.torsionConstFactor *= stiffnessFactor;
// e.firstBendingConstFactor *= stiffnessFactor;
// e.secondBendingConstFactor *= stiffnessFactor;
// }
global.initialParameters(0) = mesh->elements[0].material.youngsModulus; global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
global.initialParameters(1) = mesh->elements[0].A; global.initialParameters(1) = mesh->elements[0].A;
global.initialParameters(2) = mesh->elements[0].J; global.initialParameters(2) = mesh->elements[0].J;
global.initialParameters(3) = mesh->elements[0].I2; global.initialParameters(3) = mesh->elements[0].I2;
global.initialParameters(4) = mesh->elements[0].I3; global.initialParameters(4) = mesh->elements[0].I3;
global.initialParameters(5) = global.innerHexagonInitialPos; global.initialParameters(optimizationParamters-2) = global.innerHexagonInitialPos;
global.innerHexagonInitialRotationAngle = 30; global.innerHexagonInitialRotationAngle = 30;
global.initialParameters(6) = global.innerHexagonInitialRotationAngle; global.initialParameters(optimizationParamters-1) = global.innerHexagonInitialRotationAngle;
} }
void ReducedModelOptimizer::computeReducedModelSimulationJob( void ReducedModelOptimizer::computeReducedModelSimulationJob(
@ -522,7 +526,7 @@ void ReducedModelOptimizer::visualizeResults(
&fullPatternSimulationJobs, &fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>> const std::vector<std::shared_ptr<SimulationJob>>
&reducedPatternSimulationJobs, &reducedPatternSimulationJobs,
const std::vector<SimulationScenario> &simulationScenarios, const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex> const std::unordered_map<ReducedPatternVertexIndex, FullPatternVertexIndex>
&reducedToFullInterfaceViMap) { &reducedToFullInterfaceViMap) {
FormFinder simulator; FormFinder simulator;
@ -549,7 +553,7 @@ void ReducedModelOptimizer::visualizeResults(
simulator.executeSimulation(pReducedPatternSimulationJob); simulator.executeSimulation(pReducedPatternSimulationJob);
double normalizationFactor = 1; double normalizationFactor = 1;
if (global.optimizationSettings.normalizationStrategy != if (global.optimizationSettings.normalizationStrategy !=
Settings::NormalizationStrategy::NonNormalized) { ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) {
normalizationFactor = normalizationFactor =
global.objectiveNormalizationValues[simulationScenarioIndex]; global.objectiveNormalizationValues[simulationScenarioIndex];
} }
@ -559,7 +563,7 @@ void ReducedModelOptimizer::visualizeResults(
reducedModelResults.displacements, fullModelResults.displacements, reducedModelResults.displacements, fullModelResults.displacements,
reducedToFullInterfaceViMap, normalizationFactor); reducedToFullInterfaceViMap, normalizationFactor);
std::cout << "Error of simulation scenario " std::cout << "Error of simulation scenario "
<< simulationScenarioStrings[simulationScenarioIndex] << " is " << ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex] << " is "
<< error << std::endl; << error << std::endl;
totalError += error; totalError += error;
reducedModelResults.registerForDrawing(); reducedModelResults.registerForDrawing();
@ -570,7 +574,7 @@ void ReducedModelOptimizer::visualizeResults(
"patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/"
"Images/" + "Images/" +
pFullPatternSimulationMesh->getLabel() + "_" + pFullPatternSimulationMesh->getLabel() + "_" +
simulationScenarioStrings[simulationScenarioIndex]; ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex];
polyscope::show(); polyscope::show();
polyscope::screenshot(screenshotFilename, false); polyscope::screenshot(screenshotFilename, false);
fullModelResults.unregister(); fullModelResults.unregister();
@ -636,10 +640,10 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements(
} }
} }
ReducedModelOptimizer::Results ReducedModelOptimization::Results
ReducedModelOptimizer::runOptimization(const Settings &settings) { ReducedModelOptimizer::runOptimization(const Settings &settings) {
global.gObjectiveValueHistory.clear(); global.objectiveValueHistory.clear();
dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters); dlib::matrix<double, 0, 1> xMin(global.numberOfOptimizationParameters);
dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters); dlib::matrix<double, 0, 1> xMax(global.numberOfOptimizationParameters);
for (int i = 0; i < global.numberOfOptimizationParameters; i++) { for (int i = 0; i < global.numberOfOptimizationParameters; i++) {
@ -657,36 +661,43 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
auto end = std::chrono::system_clock::now(); auto end = std::chrono::system_clock::now();
auto elapsed = auto elapsed =
std::chrono::duration_cast<std::chrono::milliseconds>(end - start); std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
Results results; ReducedModelOptimization::Results results;
results.numberOfSimulationCrashes = global.numOfSimulationCrashes; results.numberOfSimulationCrashes = global.numOfSimulationCrashes;
results.x = global.minX; results.optimalXNameValuePairs.reserve(settings.xRanges.size());
for(int xVariableIndex=0;xVariableIndex<settings.xRanges.size();xVariableIndex++){
results.optimalXNameValuePairs[settings.xRanges[xVariableIndex].label]=
global.minX[xVariableIndex];
}
results.objectiveValue = global.minY; results.objectiveValue = global.minY;
#ifdef POLYSCOPE_DEFINED
std::cout<<"Total optimal objective value:"<<global.minY<<std::endl;
if (global.minY != result.y) { if (global.minY != result.y) {
std::cerr << "Global min objective is not equal to result objective" std::cerr << "Global min objective is not equal to result objective"
<< std::endl; << std::endl;
} }
#endif
// Compute obj value per simulation scenario and the raw objective value // Compute obj value per simulation scenario and the raw objective value
results.rawObjectiveValue=0; results.rawObjectiveValue=0;
updateMesh(results.x.size(), results.x.data()); std::vector<double> optimalX(results.optimalXNameValuePairs.size());
for(int xVariableIndex=0;xVariableIndex<global.numberOfOptimizationParameters;xVariableIndex++){
optimalX[xVariableIndex]=
global.minX[xVariableIndex];
}
updateMesh(optimalX.size(), optimalX.data());
results.objectiveValuePerSimulationScenario.resize( results.objectiveValuePerSimulationScenario.resize(
NumberOfSimulationScenarios); ReducedModelOptimization::NumberOfSimulationScenarios);
FormFinder::Settings simulationSettings; LinearSimulationModel simulator;
FormFinder simulator; for (int simulationScenarioIndex:global.simulationScenarioIndices) {
for (int simulationScenarioIndex = 0;
simulationScenarioIndex < NumberOfSimulationScenarios;
simulationScenarioIndex++) {
SimulationResults reducedModelResults = simulator.executeSimulation( SimulationResults reducedModelResults = simulator.executeSimulation(
global.reducedPatternSimulationJobs[simulationScenarioIndex], global.reducedPatternSimulationJobs[simulationScenarioIndex]);
simulationSettings);
const double error = computeError( const double error = computeError(
reducedModelResults.displacements, reducedModelResults.displacements,
global.fullPatternDisplacements[simulationScenarioIndex], global.fullPatternDisplacements[simulationScenarioIndex],
global.reducedToFullInterfaceViMap, global.reducedToFullInterfaceViMap,
global.objectiveNormalizationValues[simulationScenarioIndex]); global.objectiveNormalizationValues[simulationScenarioIndex]);
results.rawObjectiveValue+=computeRawError(reducedModelResults.displacements,global.fullPatternDisplacements[simulationScenarioIndex],global.reducedToFullInterfaceViMap); results.rawObjectiveValue+=computeRawError(reducedModelResults.displacements,global.fullPatternDisplacements[simulationScenarioIndex],global.reducedToFullInterfaceViMap);
results.objectiveValuePerSimulationScenario[simulationScenarioIndex] = results.objectiveValuePerSimulationScenario[simulationScenarioIndex] =
error; error;
} }
@ -709,7 +720,7 @@ ReducedModelOptimizer::createScenarios(
scenarios.resize(SimulationScenario::NumberOfSimulationScenarios); scenarios.resize(SimulationScenario::NumberOfSimulationScenarios);
std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices; std::unordered_map<VertexIndex, std::unordered_set<DoFType>> fixedVertices;
std::unordered_map<VertexIndex, Vector6d> nodalForces; std::unordered_map<VertexIndex, Vector6d> nodalForces;
const double forceMagnitude = 1; const double forceMagnitude = 10;
//// Axial //// Axial
SimulationScenario scenarioName = SimulationScenario::Axial; SimulationScenario scenarioName = SimulationScenario::Axial;
@ -722,7 +733,7 @@ ReducedModelOptimizer::createScenarios(
nodalForces[viPair.first] = nodalForces[viPair.first] =
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0,
0, 0}) * 0, 0}) *
forceMagnitude * 2; forceMagnitude * 10;
fixedVertices[viPair.second] = fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5}; std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
} }
@ -756,7 +767,7 @@ ReducedModelOptimizer::createScenarios(
nodalForces[viPair.first] = nodalForces[viPair.first] =
Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0,
0, 0}) * 0, 0}) *
forceMagnitude * 1; forceMagnitude * 4;
fixedVertices[viPair.second] = fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5}; std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
} }
@ -806,7 +817,7 @@ ReducedModelOptimizer::createScenarios(
fixedVertices.clear(); fixedVertices.clear();
nodalForces.clear(); nodalForces.clear();
for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) {
nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 1; nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 0.5;
fixedVertices[viPair.second] = fixedVertices[viPair.second] =
std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5}; std::unordered_set<DoFType>{0, 1, 2, 3, 4, 5};
} }
@ -814,7 +825,7 @@ ReducedModelOptimizer::createScenarios(
SimulationJob(pMesh, simulationScenarioStrings[scenarioName], SimulationJob(pMesh, simulationScenarioStrings[scenarioName],
fixedVertices, nodalForces, {})); fixedVertices, nodalForces, {}));
//// Double using moments //// Dome using moments
scenarioName = SimulationScenario::Dome; scenarioName = SimulationScenario::Dome;
fixedVertices.clear(); fixedVertices.clear();
nodalForces.clear(); nodalForces.clear();
@ -903,7 +914,7 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
} }
} }
ReducedModelOptimizer::Results ReducedModelOptimizer::optimize( Results ReducedModelOptimizer::optimize(
const Settings &optimizationSettings, const Settings &optimizationSettings,
const std::vector<SimulationScenario> &simulationScenarios) { const std::vector<SimulationScenario> &simulationScenarios) {

331
src/reducedmodeloptimizer.hpp
View File

@ -1,13 +1,14 @@
#ifndef REDUCEDMODELOPTIMIZER_HPP #ifndef REDUCEDMODELOPTIMIZER_HPP
#define REDUCEDMODELOPTIMIZER_HPP #define REDUCEDMODELOPTIMIZER_HPP
#include "beamformfinder.hpp" #include "drmsimulationmodel.hpp"
#include "csvfile.hpp" #include "csvfile.hpp"
#include "edgemesh.hpp" #include "edgemesh.hpp"
#include "elementalmesh.hpp" #include "simulationmesh.hpp"
#include "linearsimulationmodel.hpp" #include "linearsimulationmodel.hpp"
#include "matplot/matplot.h" #include "matplot/matplot.h"
#include <Eigen/Dense> #include <Eigen/Dense>
#include "reducedmodeloptimizer_structs.hpp"
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
#include "polyscope/color_management.h" #include "polyscope/color_management.h"
@ -16,6 +17,7 @@ using FullPatternVertexIndex = VertexIndex;
using ReducedPatternVertexIndex = VertexIndex; using ReducedPatternVertexIndex = VertexIndex;
class ReducedModelOptimizer { class ReducedModelOptimizer {
std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh; std::shared_ptr<SimulationMesh> m_pReducedPatternSimulationMesh;
std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh; std::shared_ptr<SimulationMesh> m_pFullPatternSimulationMesh;
std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex> std::unordered_map<FullPatternVertexIndex, ReducedPatternVertexIndex>
@ -25,318 +27,14 @@ class ReducedModelOptimizer {
std::unordered_map<size_t, size_t> nodeToSlot; std::unordered_map<size_t, size_t> nodeToSlot;
std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode; std::unordered_map<size_t, std::unordered_set<size_t>> slotToNode;
#ifdef POLYSCOPE_DEFINED #ifdef POLYSCOPE_DEFINED
struct StaticColors {
glm::vec3 fullInitial;
glm::vec3 fullDeformed;
glm::vec3 reducedInitial;
glm::vec3 reducedDeformed;
StaticColors() {
fullInitial = {0.416666, 0.109804, 0.890196};
fullDeformed = {0.583333, 0.890196, 0.109804};
reducedInitial = {0.890196, 0.109804, 0.193138};
reducedDeformed = {0.109804, 0.890196, 0.806863};
}
};
inline static StaticColors colors;
#endif // POLYSCOPE_DEFINED #endif // POLYSCOPE_DEFINED
public: public:
inline static int fanSize{6}; inline static int fanSize{6};
inline static VectorType patternPlaneNormal{0, 0, 1}; inline static VectorType patternPlaneNormal{0, 0, 1};
enum SimulationScenario { ReducedModelOptimization::Results optimize(const ReducedModelOptimization::Settings &xRanges,
Axial, const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios =
Shear, std::vector<ReducedModelOptimization::SimulationScenario>());
Bending,
Dome,
Saddle,
NumberOfSimulationScenarios
};
struct xRange{
std::string label;
double min;
double max;
std::string toString() const {
return label + "=[" + std::to_string(min) + "," + std::to_string(max) +
"]";
}
};
struct Results;
struct Settings {
enum NormalizationStrategy {
NonNormalized,
Epsilon,
MaxDisplacement,
EqualDisplacements
};
inline static vector<std::string> normalizationStrategyStrings{
"NonNormalized", "Epsilon", "MaxDsiplacement", "EqualDisplacements"};
std::vector<xRange> xRanges;
int numberOfFunctionCalls{100};
double solutionAccuracy{1e-2};
NormalizationStrategy normalizationStrategy{Epsilon};
double normalizationParameter{0.003};
std::string toString() const {
std::string settingsString;
if (!xRanges.empty()) {
std::string xRangesString;
for (const xRange &range : xRanges) {
xRangesString += range.toString() + " ";
}
settingsString += xRangesString;
}
settingsString +=
"FuncCalls=" + std::to_string(numberOfFunctionCalls) +
" Accuracy=" + std::to_string(solutionAccuracy) +
" Norm=" + normalizationStrategyStrings[normalizationStrategy];
return settingsString;
}
void writeHeaderTo(csvFile &os) const {
if (!xRanges.empty()) {
for (const xRange &range : xRanges) {
os << range.label + " max";
os << range.label + " min";
}
}
os << "Function Calls";
os << "Solution Accuracy";
os << "Normalization strategy";
// os << std::endl;
}
void writeSettingsTo(csvFile &os) const {
if (!xRanges.empty()) {
for (const xRange &range : xRanges) {
os << range.max;
os << range.min;
}
}
os << numberOfFunctionCalls;
os << solutionAccuracy;
os << normalizationStrategyStrings[normalizationStrategy] + "_" +
std::to_string(normalizationParameter);
}
};
struct Results {
double time{-1};
int numberOfSimulationCrashes{0};
std::vector<double> x;
double objectiveValue;
double rawObjectiveValue;
std::vector<double> objectiveValuePerSimulationScenario;
std::vector<std::shared_ptr<SimulationJob>> fullPatternSimulationJobs;
std::vector<std::shared_ptr<SimulationJob>> reducedPatternSimulationJobs;
void save(const string &saveToPath) const {
assert(std::filesystem::is_directory(saveToPath));
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
std::filesystem::path simulationJobFolderPath(
std::filesystem::path(saveToPath)
.append(pFullPatternSimulationJob->label));
std::filesystem::create_directory(simulationJobFolderPath);
const auto fullPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Full");
std::filesystem::create_directory(fullPatternDirectoryPath);
pFullPatternSimulationJob->save(fullPatternDirectoryPath.string());
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
const auto reducedPatternDirectoryPath =
std::filesystem::path(simulationJobFolderPath).append("Reduced");
if (!std::filesystem::exists(reducedPatternDirectoryPath)) {
std::filesystem::create_directory(reducedPatternDirectoryPath);
}
pReducedPatternSimulationJob->save(
reducedPatternDirectoryPath.string());
}
}
void load(const string &loadFromPath) {
assert(std::filesystem::is_directory(loadFromPath));
for (const auto &directoryEntry :
filesystem::directory_iterator(loadFromPath)) {
const auto simulationScenarioPath = directoryEntry.path();
if (!std::filesystem::is_directory(simulationScenarioPath)) {
continue;
}
// Load reduced pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath)
.append("Full"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
fullPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
}
// Load full pattern files
for (const auto &fileEntry : filesystem::directory_iterator(
std::filesystem::path(simulationScenarioPath)
.append("Reduced"))) {
const auto filepath = fileEntry.path();
if (filepath.extension() == ".json") {
SimulationJob job;
job.load(filepath.string());
reducedPatternSimulationJobs.push_back(
std::make_shared<SimulationJob>(job));
}
}
}
}
#if POLYSCOPE_DEFINED
void draw() const {
initPolyscope();
FormFinder simulator;
LinearSimulationModel linearSimulator;
assert(fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->registerForDrawing(
colors.fullInitial);
reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(
colors.reducedInitial, false);
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
pFullPatternSimulationJob->registerForDrawing(
fullPatternSimulationJobs[0]->pMesh->getLabel());
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.registerForDrawing(colors.fullDeformed);
SimulationResults fullModelLinearResults =
linearSimulator.executeSimulation(pFullPatternSimulationJob);
fullModelLinearResults.setLabelPrefix("linear");
fullModelLinearResults.registerForDrawing(colors.fullDeformed, false);
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.registerForDrawing(colors.reducedDeformed);
SimulationResults reducedModelLinearResults =
linearSimulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelLinearResults.setLabelPrefix("linear");
reducedModelLinearResults.registerForDrawing(colors.reducedDeformed,
false);
polyscope::options::programName =
fullPatternSimulationJobs[0]->pMesh->getLabel();
polyscope::view::resetCameraToHomeView();
polyscope::show();
// Save a screensh
const std::string screenshotFilename =
"/home/iason/Coding/Projects/Approximating shapes with flat "
"patterns/RodModelOptimizationForPatterns/Results/Images/" +
fullPatternSimulationJobs[0]->pMesh->getLabel() + "_" +
pFullPatternSimulationJob->getLabel();
polyscope::screenshot(screenshotFilename, false);
fullModelResults.unregister();
reducedModelResults.unregister();
reducedModelLinearResults.unregister();
fullModelLinearResults.unregister();
// double error = computeError(
// reducedModelResults.displacements,fullModelResults.displacements,
// );
// std::cout << "Error of simulation scenario "
// <<
// simula simulationScenarioStrings[simulationScenarioIndex]
// << " is "
// << error << std::endl;
}
}
#endif // POLYSCOPE_DEFINED
void saveMeshFiles() const {
const int numberOfSimulationJobs = fullPatternSimulationJobs.size();
assert(numberOfSimulationJobs != 0 &&
fullPatternSimulationJobs.size() ==
reducedPatternSimulationJobs.size());
fullPatternSimulationJobs[0]->pMesh->savePly(
"FullPattern_undeformed.ply");
reducedPatternSimulationJobs[0]->pMesh->savePly(
"ReducedPattern_undeformed.ply");
FormFinder simulator;
for (int simulationJobIndex = 0;
simulationJobIndex < numberOfSimulationJobs; simulationJobIndex++) {
// Drawing of full pattern results
const std::shared_ptr<SimulationJob> &pFullPatternSimulationJob =
fullPatternSimulationJobs[simulationJobIndex];
SimulationResults fullModelResults =
simulator.executeSimulation(pFullPatternSimulationJob);
fullModelResults.saveDeformedModel();
// Drawing of reduced pattern results
const std::shared_ptr<SimulationJob> &pReducedPatternSimulationJob =
reducedPatternSimulationJobs[simulationJobIndex];
SimulationResults reducedModelResults =
simulator.executeSimulation(pReducedPatternSimulationJob);
reducedModelResults.saveDeformedModel();
}
}
void
writeHeaderTo(const ReducedModelOptimizer::Settings &settings_optimization,
csvFile &os) {
os << "Total raw Obj value";
os << "Total Obj value";
for (int simulationScenarioIndex = 0;
simulationScenarioIndex <
SimulationScenario::NumberOfSimulationScenarios;
simulationScenarioIndex++) {
os << "Obj value " + simulationScenarioStrings[simulationScenarioIndex];
}
for (const ReducedModelOptimizer::xRange &range :
settings_optimization.xRanges) {
os << range.label;
}
os << "Time(s)";
os << "#Crashes";
}
void
writeResultsTo(const ReducedModelOptimizer::Settings &settings_optimization,
csvFile &os) const {
os << rawObjectiveValue;
os << objectiveValue;
for (double scenarioObjValue : objectiveValuePerSimulationScenario) {
os << scenarioObjValue;
}
for (const double &optimalX : x) {
os << optimalX;
}
for (int unusedXVarCounter = 0;
unusedXVarCounter < settings_optimization.xRanges.size() - x.size();
unusedXVarCounter++) {
os << "-";
}
os << time;
if (numberOfSimulationCrashes == 0) {
os << "-";
} else {
os << numberOfSimulationCrashes;
}
}
};
inline static const std::string simulationScenarioStrings[] = {
"Axial", "Shear", "Bending", "Dome", "Saddle"};
Results optimize(const Settings &xRanges,
const std::vector<SimulationScenario> &simulationScenarios =
std::vector<SimulationScenario>());
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);
@ -348,9 +46,8 @@ public:
SimulationJob SimulationJob
getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
void initializePatterns( void initializePatterns(PatternGeometry &fullPattern, PatternGeometry &reducedPatterm,
PatternGeometry &fullPattern, PatternGeometry &reducedPatterm, const std::unordered_set<size_t> &reducedModelExcludedEges, const int &optimizationParameters);
const std::unordered_set<size_t> &reducedModelExcludedEges);
static void runSimulation(const std::string &filename, static void runSimulation(const std::string &filename,
std::vector<double> &x); std::vector<double> &x);
@ -383,7 +80,7 @@ public:
&fullPatternSimulationJobs, &fullPatternSimulationJobs,
const std::vector<std::shared_ptr<SimulationJob>> const std::vector<std::shared_ptr<SimulationJob>>
&reducedPatternSimulationJobs, &reducedPatternSimulationJobs,
const std::vector<SimulationScenario> &simulationScenarios, const std::vector<ReducedModelOptimization::SimulationScenario> &simulationScenarios,
const std::unordered_map<ReducedPatternVertexIndex, const std::unordered_map<ReducedPatternVertexIndex,
FullPatternVertexIndex> FullPatternVertexIndex>
&reducedToFullInterfaceViMap); &reducedToFullInterfaceViMap);
@ -409,12 +106,13 @@ public:
&reducedToFullInterfaceViMap, &reducedToFullInterfaceViMap,
const double &normalizationFactor); const double &normalizationFactor);
private: static double objective(double E, double A, double innerHexagonSize, double innerHexagonRotationAngle);
private:
static void computeDesiredReducedModelDisplacements( static void computeDesiredReducedModelDisplacements(
const SimulationResults &fullModelResults, const SimulationResults &fullModelResults,
const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap, const std::unordered_map<size_t, size_t> &displacementsReducedToFullMap,
Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); Eigen::MatrixX3d &optimalDisplacementsOfReducedModel);
static Results runOptimization(const Settings &settings); static ReducedModelOptimization::Results runOptimization(const ReducedModelOptimization::Settings &settings);
std::vector<std::shared_ptr<SimulationJob>> std::vector<std::shared_ptr<SimulationJob>>
createScenarios(const std::shared_ptr<SimulationMesh> &pMesh); createScenarios(const std::shared_ptr<SimulationMesh> &pMesh);
void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern, void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern,
@ -422,10 +120,11 @@ private:
void createSimulationMeshes(PatternGeometry &fullModel, void createSimulationMeshes(PatternGeometry &fullModel,
PatternGeometry &reducedModel); PatternGeometry &reducedModel);
static void static void
initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh); initializeOptimizationParameters(const std::shared_ptr<SimulationMesh> &mesh, const int &optimizationParamters);
static double objective(long n, const double *x); static double objective(long n, const double *x);
FormFinder simulator; FormFinder simulator;
void computeObjectiveValueNormalizationFactors(); void computeObjectiveValueNormalizationFactors();
}; };
void updateMesh(long n, const double *x);
#endif // REDUCEDMODELOPTIMIZER_HPP #endif // REDUCEDMODELOPTIMIZER_HPP