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

Optimizing the pattern using E,A,J,I2,I3,hexSize,hexAngle

This commit is contained in:
iasonmanolas 2021-03-16 11:57:27 +02:00
parent 711e72551d
commit 07018be421
4 changed files with 51 additions and 129 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
src/linearsimulationmodel.cpp
View File

@ -1,21 +0,0 @@
#include "linearsimulationmodel.hpp"
#include <Eigen/Core>
#include <filesystem>
//#include <igl/list_to_matrix.h>
#include <iostream>
#include <vcg/complex/algorithms/create/platonic.h>
#include <vcg/complex/algorithms/update/normal.h>

19
src/main.cpp
View File

@ -26,27 +26,30 @@ int main(int argc, char *argv[]) {
}
// Populate the pattern pair to be optimized
const int interfaceNodeIndex=3;
////Full pattern
const std::string filepath_fullPattern = argv[1];
PatternGeometry fullPattern(filepath_fullPattern);
fullPattern.setLabel(
std::filesystem::path(filepath_fullPattern).stem().string());
fullPattern.scale(0.03);
fullPattern.scale(0.03,interfaceNodeIndex);
////Reduced pattern
const std::string filepath_reducedPattern = argv[2];
PatternGeometry reducedPattern(filepath_reducedPattern);
reducedPattern.setLabel(
std::filesystem::path(filepath_reducedPattern).stem().string());
reducedPattern.scale(0.03);
reducedPattern.scale(0.03,interfaceNodeIndex);
// Set the optization settings
ReducedModelOptimizer::xRange beamWidth{"B", 0.5, 1.5};
ReducedModelOptimizer::xRange beamDimensionsRatio{"bOverh", 0.7, 1.3};
ReducedModelOptimizer::xRange beamE{"E", 0.1, 1.9};
ReducedModelOptimizer::xRange beamA{"A", 0.25, 2.25};
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};
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};
ReducedModelOptimizer::xRange beamJ{"J", beamI2.min+beamI3.min, beamI2.max+beamI3.max};
ReducedModelOptimizer::xRange innerHexagonSize{"HexSize", 0.1, 0.8};
ReducedModelOptimizer::xRange innerHexagonAngle{"HexAngle", -29.5, 29.5};
ReducedModelOptimizer::Settings settings_optimization;
settings_optimization.xRanges = {beamWidth, beamDimensionsRatio, beamE,
settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3,
innerHexagonSize, innerHexagonAngle};
const bool input_numberOfFunctionCallsDefined = argc >= 4;
settings_optimization.numberOfFunctionCalls =
@ -61,6 +64,7 @@ int main(int argc, char *argv[]) {
fullPattern.getLabel() + "@" + reducedPattern.getLabel();
const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]);
ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
optimizer.initializePatterns(fullPattern, reducedPattern, {});
ReducedModelOptimizer::Results optimizationResults =
@ -102,7 +106,8 @@ int main(int argc, char *argv[]) {
settings_optimization.writeSettingsTo(csv_results);
csv_results << endrow;
// optimizationResults.draw();
#ifdef POLYSCOPE_DEFINED
optimizationResults.draw();
#endif
return 0;
}

138
src/reducedmodeloptimizer.cpp
View File

@ -17,12 +17,12 @@ struct GlobalOptimizationVariables {
reducedToFullInterfaceViMap;
matplot::line_handle gPlotHandle;
std::vector<double> gObjectiveValueHistory;
Eigen::VectorXd g_initialParameters;
Eigen::VectorXd initialParameters;
std::vector<ReducedModelOptimizer::SimulationScenario>
simulationScenarioIndices;
std::vector<VectorType> g_innerHexagonVectors{6, VectorType(0, 0, 0)};
double innerHexagonInitialRotationAngle{30};
double g_innerHexagonInitialPos{0};
double innerHexagonInitialPos{0};
double minY{DBL_MAX};
std::vector<double> minX;
int numOfSimulationCrashes{false};
@ -85,58 +85,35 @@ void updateMesh(long n, const double *x) {
std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh =
global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]
->pMesh;
// const Element &elem = g_reducedPatternSimulationJob[0]->mesh->elements[0];
// std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
// "
// << elem.firstBendingConstFactor << " "
// << elem.secondBendingConstFactor << std::endl;
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];
// if (g_reducedPatternExludedEdges.contains(ei)) {
// continue;
// }
// e.properties.E = g_initialParameters * x[ei];
// e.properties.E = g_initialParameters(0) * x[0];
// e.properties.G = g_initialParameters(1) * x[1];
e.setDimensions(
RectangularBeamDimensions(global.g_initialParameters(0) * x[0],
global.g_initialParameters(0) * x[0] /
(global.g_initialParameters(1) * x[1])));
RectangularBeamDimensions(beamWidth,
beamHeight));
e.setMaterial(ElementMaterial(e.material.poissonsRatio,
global.g_initialParameters(2) * x[2]));
// e.properties.A = g_initialParameters(0) * x[0];
// e.properties.J = g_initialParameters(1) * x[1];
// e.properties.I2 = g_initialParameters(2) * x[2];
// e.properties.I3 = g_initialParameters(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;
E));
e.A = A;
e.J = J;
e.I2 = I2;
e.I3 = I3;
}
// std::cout << elem.axialConstFactor << " " << elem.torsionConstFactor << "
// "
// << elem.firstBendingConstFactor << " "
// << elem.secondBendingConstFactor << std::endl;
// const Element &e = pReducedPatternSimulationMesh->elements[0];
// std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " "
// << e.firstBendingConstFactor << " " <<
// e.secondBendingConstFactor
// << std::endl;
assert(pReducedPatternSimulationMesh->EN() == 12);
auto R = vcg::RotationMatrix(
ReducedModelOptimizer::patternPlaneNormal,
vcg::math::ToRad(x[4] - global.innerHexagonInitialRotationAngle));
// for (VertexIndex vi = 0; vi < pReducedPatternSimulationMesh->VN();
// vi += 2) {
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[3];
R * global.g_innerHexagonVectors[rotationCounter] * x[5];
}
pReducedPatternSimulationMesh->reset();
@ -150,10 +127,10 @@ double ReducedModelOptimizer::objective(double b, double r, double E) {
return ReducedModelOptimizer::objective(x.size(), x.data());
}
double ReducedModelOptimizer::objective(double b, double h, double E,
double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
double innerHexagonSize,
double innerHexagonRotationAngle) {
std::vector<double> x{b, h, E, innerHexagonSize, innerHexagonRotationAngle};
std::vector<double> x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle};
return ReducedModelOptimizer::objective(x.size(), x.data());
}
@ -470,9 +447,9 @@ void ReducedModelOptimizer::initializePatterns(
}
const double innerHexagonInitialPos_x =
copyReducedPattern.vert[0].cP()[0] / global.g_innerHexagonVectors[0][0];
const double innerHexagonInitialPos_y =
copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
global.g_innerHexagonInitialPos = innerHexagonInitialPos_x;
// const double innerHexagonInitialPos_y =
// copyReducedPattern.vert[0].cP()[1] / global.g_innerHexagonVectors[0][1];
global.innerHexagonInitialPos = innerHexagonInitialPos_x;
global.innerHexagonInitialRotationAngle =
30; /* NOTE: Here I assume that the CW reduced pattern is given as input.
This is not very generic */
@ -483,8 +460,8 @@ void ReducedModelOptimizer::initializePatterns(
void ReducedModelOptimizer::initializeOptimizationParameters(
const std::shared_ptr<SimulationMesh> &mesh) {
global.numberOfOptimizationParameters = 5;
global.g_initialParameters.resize(global.numberOfOptimizationParameters);
global.numberOfOptimizationParameters = 7;
global.initialParameters.resize(global.numberOfOptimizationParameters);
// Save save the beam stiffnesses
// for (size_t ei = 0; ei < pReducedModelElementalMesh->EN(); ei++) {
// Element &e = pReducedModelElementalMesh->elements[ei];
@ -495,25 +472,14 @@ void ReducedModelOptimizer::initializeOptimizationParameters(
// e.firstBendingConstFactor *= stiffnessFactor;
// e.secondBendingConstFactor *= stiffnessFactor;
// }
const double initialB = std::sqrt(mesh->elements[0].A);
const double initialRatio = 1;
global.g_initialParameters(0) = initialB;
global.g_initialParameters(1) = initialRatio;
global.g_initialParameters(2) = mesh->elements[0].material.youngsModulus;
global.g_initialParameters(3) = global.g_innerHexagonInitialPos;
global.initialParameters(0) = mesh->elements[0].material.youngsModulus;
global.initialParameters(1) = mesh->elements[0].A;
global.initialParameters(2) = mesh->elements[0].J;
global.initialParameters(3) = mesh->elements[0].I2;
global.initialParameters(4) = mesh->elements[0].I3;
global.initialParameters(5) = global.innerHexagonInitialPos;
global.innerHexagonInitialRotationAngle = 30;
global.g_initialParameters(4) = global.innerHexagonInitialRotationAngle;
// g_initialParameters =
// m_pReducedPatternSimulationMesh->elements[0].properties.E;
// for (size_t ei = 0; ei < m_pReducedPatternSimulationMesh->EN(); ei++) {
// }
// g_initialParameters(0) = mesh->elements[0].properties.E;
// g_initialParameters(1) = mesh->elements[0].properties.G;
// g_initialParameters(0) = mesh->elements[0].properties.A;
// g_initialParameters(1) = mesh->elements[0].properties.J;
// g_initialParameters(2) = mesh->elements[0].properties.I2;
// g_initialParameters(3) = mesh->elements[0].properties.I3;
// }
global.initialParameters(6) = global.innerHexagonInitialRotationAngle;
}
void ReducedModelOptimizer::computeReducedModelSimulationJob(
@ -683,7 +649,7 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
auto start = std::chrono::system_clock::now();
dlib::function_evaluation result;
double (*objF)(double, double, double, double, double) = &objective;
double (*objF)(double, double, double, double, double,double,double) = &objective;
result = dlib::find_min_global(
objF, xMin, xMax,
dlib::max_function_calls(settings.numberOfFunctionCalls),
@ -700,20 +666,9 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
<< std::endl;
}
// Compute raw objective value
if (global.optimizationSettings.normalizationStrategy !=
Settings::NormalizationStrategy::NonNormalized) {
auto temp = global.optimizationSettings.normalizationStrategy;
global.optimizationSettings.normalizationStrategy =
Settings::NormalizationStrategy::NonNormalized;
results.rawObjectiveValue = objective(results.x.size(), results.x.data());
global.optimizationSettings.normalizationStrategy = temp;
} else {
results.rawObjectiveValue = results.objectiveValue;
}
// Compute obj value per simulation scenario
// Compute obj value per simulation scenario and the raw objective value
results.rawObjectiveValue=0;
updateMesh(results.x.size(), results.x.data());
results.objectiveValuePerSimulationScenario.resize(
NumberOfSimulationScenarios);
@ -730,18 +685,11 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) {
global.fullPatternDisplacements[simulationScenarioIndex],
global.reducedToFullInterfaceViMap,
global.objectiveNormalizationValues[simulationScenarioIndex]);
results.rawObjectiveValue+=computeRawError(reducedModelResults.displacements,global.fullPatternDisplacements[simulationScenarioIndex],global.reducedToFullInterfaceViMap);
results.objectiveValuePerSimulationScenario[simulationScenarioIndex] =
error;
}
// if (result.y !=
// std::accumulate(results.objectiveValuePerSimulationScenario.begin(),
// results.objectiveValuePerSimulationScenario.end(), 0))
// {
// std::cerr
// << "Sum of per scenario objectives is not equal to result objective"
// << std::endl;
// }
results.time = elapsed.count() / 1000.0;
const bool printDebugInfo = false;
if (printDebugInfo) {
@ -948,15 +896,8 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() {
Settings::NormalizationStrategy::Epsilon) {
const double epsilon =
global.optimizationSettings.normalizationParameter;
if (epsilon > fullPatternDisplacementNormSum[simulationScenarioIndex]) {
// std::cout << "Epsilon used in "
// <<
// simulationScenarioStrings[simulationScenarioIndex]
// << std::endl;
}
global.objectiveNormalizationValues[simulationScenarioIndex] = std::max(
fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon);
// displacementNormSum;
}
}
}
@ -1017,9 +958,6 @@ ReducedModelOptimizer::Results ReducedModelOptimizer::optimize(
optResults.reducedPatternSimulationJobs.push_back(
global.reducedPatternSimulationJobs[simulationScenarioIndex]);
}
#ifdef POLYSCOPE_DEFINED
optResults.draw();
#endif // POLYSCOPE_DEFINED
return optResults;
}

2
src/reducedmodeloptimizer.hpp
View File

@ -355,7 +355,7 @@ public:
static void runSimulation(const std::string &filename,
std::vector<double> &x);
static double objective(double x0, double x1, double x2, double x3,
static double objective(double x2, double A, double J, double I2, double I3, double x3,
double innerHexagonRotationAngle);
static double objective(double b, double r, double E);