From 07018be421607ec33fc348893f3ebfca91c74925 Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 16 Mar 2021 11:57:27 +0200 Subject: [PATCH] Optimizing the pattern using E,A,J,I2,I3,hexSize,hexAngle --- src/linearsimulationmodel.cpp | 21 ------ src/main.cpp | 19 +++-- src/reducedmodeloptimizer.cpp | 138 ++++++++++------------------------ src/reducedmodeloptimizer.hpp | 2 +- 4 files changed, 51 insertions(+), 129 deletions(-) diff --git a/src/linearsimulationmodel.cpp b/src/linearsimulationmodel.cpp index 6e663c6..e69de29 100644 --- a/src/linearsimulationmodel.cpp +++ b/src/linearsimulationmodel.cpp @@ -1,21 +0,0 @@ -#include "linearsimulationmodel.hpp" -#include -#include -//#include -#include -#include -#include - - - - - - - - - - - - - - diff --git a/src/main.cpp b/src/main.cpp index 33606fe..a8eaf8e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -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 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; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 6cba086..ee6815a 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -17,12 +17,12 @@ struct GlobalOptimizationVariables { reducedToFullInterfaceViMap; matplot::line_handle gPlotHandle; std::vector gObjectiveValueHistory; - Eigen::VectorXd g_initialParameters; + Eigen::VectorXd initialParameters; std::vector simulationScenarioIndices; std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; double innerHexagonInitialRotationAngle{30}; - double g_innerHexagonInitialPos{0}; + double innerHexagonInitialPos{0}; double minY{DBL_MAX}; std::vector minX; int numOfSimulationCrashes{false}; @@ -85,58 +85,35 @@ void updateMesh(long n, const double *x) { std::shared_ptr &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 x{b, h, E, innerHexagonSize, innerHexagonRotationAngle}; + std::vector 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 &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; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index fd5514c..400db52 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -355,7 +355,7 @@ public: static void runSimulation(const std::string &filename, std::vector &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);