diff --git a/CMakeLists.txt b/CMakeLists.txt index ed65807..253e82a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,7 +29,6 @@ if(${CMAKE_BUILD_TYPE} STREQUAL "Release") set(USE_POLYSCOPE FALSE) set(MYSOURCES_STATIC_LINK TRUE) else() - add_compile_definitions(POLYSCOPE_DEFINED) set(USE_POLYSCOPE TRUE) add_compile_definitions(POLYSCOPE_DEFINED) set(MYSOURCES_STATIC_LINK FALSE) diff --git a/src/main.cpp b/src/main.cpp index 7097a73..f4a9f19 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -47,8 +47,8 @@ int main(int argc, char *argv[]) { ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95}; ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0}; ReducedPatternOptimization::Settings settings_optimization; - settings_optimization.xRanges - = {/*beamE,*/ beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle}; + settings_optimization.parameterRanges + = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle}; const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) @@ -66,7 +66,7 @@ int main(int argc, char *argv[]) { settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]); std::string xConcatNames; - for (const auto &x : settings_optimization.xRanges) { + for (const auto &x : settings_optimization.parameterRanges) { xConcatNames.append(x.label + "_"); } xConcatNames.pop_back(); @@ -117,7 +117,7 @@ int main(int argc, char *argv[]) { const std::vector numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1}; assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]); ReducedModelOptimizer optimizer(numberOfNodesPerSlot); - optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges); + optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.parameterRanges); optimizer.optimize(settings_optimization, optimizationResults); optimizationResults.label = optimizationName; optimizationResults.baseTriangleFullPattern.copy(fullPattern); diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 74d0edc..c64d26b 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -8,22 +8,22 @@ using namespace ReducedPatternOptimization; -struct GlobalOptimizationVariables { +struct GlobalOptimizationVariables +{ // std::vector> fullPatternDisplacements; std::vector fullPatternResults; std::vector translationalDisplacementNormalizationValues; std::vector rotationalDisplacementNormalizationValues; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; - std::unordered_map - reducedToFullInterfaceViMap; + std::unordered_map reducedToFullInterfaceViMap; std::vector> fullPatternInterfaceViPairs; matplot::line_handle gPlotHandle; std::vector objectiveValueHistory_iteration; std::vector objectiveValueHistory; std::vector plotColors; - Eigen::VectorXd initialParameters; + std::array parametersInitialValue; std::vector simulationScenarioIndices; double minY{DBL_MAX}; std::vector minX; @@ -42,9 +42,10 @@ struct GlobalOptimizationVariables { double desiredMaxDisplacementValue; double desiredMaxRotationAngle; std::string currentScenarioName; - std::vector &pReducedPatternSimulationMesh, - const double &newValue)>> - updateReducedPatternFunctions_material; + std::array &pReducedPatternSimulationMesh)>, + 7> + functions_updateReducedPatternParameter; } global; @@ -147,15 +148,20 @@ double ReducedModelOptimizer::computeError( // return ReducedModelOptimizer::objective(x.size(), x.data()); //} +double ReducedModelOptimizer::objective(const double &xValue) +{ + dlib::matrix x{xValue}; + return objective(x); +} double ReducedModelOptimizer::objective(const dlib::matrix &x) { - // std::cout.precision(17); - // for (int i = 0; i < x.size(); i++) { - // std::cout << x(i) << " "; - // } - // std::cout << std::endl; + std::cout.precision(17); + for (int i = 0; i < x.size(); i++) { + std::cout << x(i) << " "; + } + std::cout << std::endl; - // std::cout << x[n - 2] << " " << x[n - 1] << std::endl; + // std::cout << x(x.size() - 2) << " " << x(x.size() - 1) << std::endl; // const Element &e = // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << // e.axialConstFactor << " " << e.torsionConstFactor << " " @@ -253,13 +259,13 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) } } // std::cout << error << std::endl; -// global.objectiveValueHistory.push_back(totalError); -// global.plotColors.push_back(10); + // global.objectiveValueHistory.push_back(totalError); + // global.plotColors.push_back(10); ++global.numberOfFunctionCalls; if (totalError < global.minY) { global.minY = totalError; - global.objectiveValueHistory.push_back(totalError); - global.objectiveValueHistory_iteration.push_back(global.numberOfFunctionCalls); + // global.objectiveValueHistory.push_back(totalError); + // global.objectiveValueHistory_iteration.push_back(global.numberOfFunctionCalls); // std::cout << "New best:" << totalError << std::endl; // // global.minX.assign(x.begin(), x.begin() + n); // std::cout.precision(17); @@ -294,33 +300,33 @@ double ReducedModelOptimizer::objective(const dlib::matrix &x) } void ReducedModelOptimizer::createSimulationMeshes( - PatternGeometry &fullModel, PatternGeometry &reducedModel, + PatternGeometry &fullModel, + PatternGeometry &reducedModel, std::shared_ptr &pFullPatternSimulationMesh, - std::shared_ptr &pReducedPatternSimulationMesh) { + std::shared_ptr &pReducedPatternSimulationMesh) +{ if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { std::cerr << "Error: A rectangular cross section is expected." << std::endl; terminate(); } - const double ym = 1 * 1e9; // Full pattern pFullPatternSimulationMesh = std::make_shared(fullModel); - pFullPatternSimulationMesh->setBeamCrossSection( - CrossSectionType{0.002, 0.002}); - pFullPatternSimulationMesh->setBeamMaterial(0.3, ym); + pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); + pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); // Reduced pattern - pReducedPatternSimulationMesh = - std::make_shared(reducedModel); - pReducedPatternSimulationMesh->setBeamCrossSection( - CrossSectionType{0.002, 0.002}); - pReducedPatternSimulationMesh->setBeamMaterial(0.3, ym); + pReducedPatternSimulationMesh = std::make_shared(reducedModel); + pReducedPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); + pReducedPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); } -void ReducedModelOptimizer::createSimulationMeshes( - PatternGeometry &fullModel, PatternGeometry &reducedModel) { - ReducedModelOptimizer::createSimulationMeshes( - fullModel, reducedModel, m_pFullPatternSimulationMesh, - m_pReducedPatternSimulationMesh); +void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel, + PatternGeometry &reducedModel) +{ + ReducedModelOptimizer::createSimulationMeshes(fullModel, + reducedModel, + m_pFullPatternSimulationMesh, + m_pReducedPatternSimulationMesh); } void ReducedModelOptimizer::computeMaps( @@ -336,24 +342,22 @@ void ReducedModelOptimizer::computeMaps( { // Compute the offset between the interface nodes const size_t interfaceSlotIndex = 4; // bottom edge - assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && - slotToNode.find(interfaceSlotIndex)->second.size() == 1); + 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 size_t baseTriangleInterfaceVi = *(slotToNode.find(interfaceSlotIndex)->second.begin()); vcg::tri::Allocator::PointerUpdater pu_fullModel; fullPattern.deleteDanglingVertices(pu_fullModel); - const size_t fullModelBaseTriangleInterfaceVi = - pu_fullModel.remap.empty() ? baseTriangleInterfaceVi - : pu_fullModel.remap[baseTriangleInterfaceVi]; + const size_t fullModelBaseTriangleInterfaceVi = pu_fullModel.remap.empty() + ? baseTriangleInterfaceVi + : pu_fullModel.remap[baseTriangleInterfaceVi]; const size_t fullModelBaseTriangleVN = fullPattern.VN(); fullPattern.createFan(); - const size_t duplicateVerticesPerFanPair = - fullModelBaseTriangleVN - fullPattern.VN() / 6; - const size_t fullPatternInterfaceVertexOffset = - fullModelBaseTriangleVN - duplicateVerticesPerFanPair; + const size_t duplicateVerticesPerFanPair = fullModelBaseTriangleVN - fullPattern.VN() / 6; + const size_t fullPatternInterfaceVertexOffset = fullModelBaseTriangleVN + - duplicateVerticesPerFanPair; // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << // std::endl; @@ -372,8 +376,7 @@ void ReducedModelOptimizer::computeMaps( reducedToFullInterfaceViMap.clear(); vcg::tri::Allocator::PointerUpdater pu_reducedModel; reducedPattern.deleteDanglingVertices(pu_reducedModel); - const size_t reducedModelBaseTriangleInterfaceVi = - pu_reducedModel.remap[baseTriangleInterfaceVi]; + const size_t reducedModelBaseTriangleInterfaceVi = pu_reducedModel.remap[baseTriangleInterfaceVi]; const size_t reducedModelInterfaceVertexOffset = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/; Results::applyOptimizationResults_innerHexagon(initialHexagonSize, @@ -382,10 +385,9 @@ void ReducedModelOptimizer::computeMaps( reducedPattern); reducedPattern.createFan({0}); //TODO: should be an input parameter from main for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { - reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + - reducedModelBaseTriangleInterfaceVi] = - fullModelBaseTriangleInterfaceVi + - fanIndex * fullPatternInterfaceVertexOffset; + reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + + reducedModelBaseTriangleInterfaceVi] + = fullModelBaseTriangleInterfaceVi + fanIndex * fullPatternInterfaceVertexOffset; } fullToReducedInterfaceViMap.clear(); constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap); @@ -420,8 +422,7 @@ void ReducedModelOptimizer::computeMaps( // ->setEnabled(true); // polyscope::show(); - std::vector nodeColorsOpposite(fullPattern.VN(), - glm::vec3(0, 0, 0)); + std::vector nodeColorsOpposite(fullPattern.VN(), glm::vec3(0, 0, 0)); for (const std::pair oppositeVerts : fullPatternOppositeInterfaceViPairs) { auto color = polyscope::getNextUniqueColor(); nodeColorsOpposite[oppositeVerts.first] = color; @@ -435,24 +436,19 @@ void ReducedModelOptimizer::computeMaps( std::vector nodeColorsReducedToFull_reduced(reducedPattern.VN(), glm::vec3(0, 0, 0)); - std::vector nodeColorsReducedToFull_full(fullPattern.VN(), - glm::vec3(0, 0, 0)); + std::vector nodeColorsReducedToFull_full(fullPattern.VN(), glm::vec3(0, 0, 0)); for (size_t vi = 0; vi < reducedPattern.VN(); vi++) { if (global.reducedToFullInterfaceViMap.contains(vi)) { - auto color = polyscope::getNextUniqueColor(); nodeColorsReducedToFull_reduced[vi] = color; - nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = - color; + nodeColorsReducedToFull_full[global.reducedToFullInterfaceViMap[vi]] = color; } } polyscope::getCurveNetwork(reducedPattern.getLabel()) - ->addNodeColorQuantity("reducedToFull_reduced", - nodeColorsReducedToFull_reduced) + ->addNodeColorQuantity("reducedToFull_reduced", nodeColorsReducedToFull_reduced) ->setEnabled(true); polyscope::getCurveNetwork(fullPattern.getLabel()) - ->addNodeColorQuantity("reducedToFull_full", - nodeColorsReducedToFull_full) + ->addNodeColorQuantity("reducedToFull_full", nodeColorsReducedToFull_full) ->setEnabled(true); polyscope::show(); } @@ -514,47 +510,88 @@ void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern, void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() { - function_updateReducedPattern_geometry = [&](const dlib::matrix &x, - std::shared_ptr - &pReducedPatternSimulationMesh) { - const int n = x.size(); - assert(n >= 2); + // //TODO:Maybe I should split theta from R + // function_updateReducedPattern_geometry = [&](const dlib::matrix &x, + // std::shared_ptr + // &pReducedPatternSimulationMesh) { + // const int n = x.size(); + // assert(n >= 2); - // CoordType center_barycentric(1, 0, 0); - // CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); - // CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric) - // * x[n - 2]); + // // CoordType center_barycentric(1, 0, 0); + // // CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); + // // CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric) + // // * x[n - 2]); - CoordType movableVertex_barycentric(1 - x(n - 2), x(n - 2) / 2, x(n - 2) / 2); - CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0) - * movableVertex_barycentric[0] - + global.baseTriangle.cP(1) - * movableVertex_barycentric[1] - + global.baseTriangle.cP(2) - * movableVertex_barycentric[2]; - baseTriangleMovableVertexPosition - = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, - vcg::math::ToRad(x(n - 1))) - * baseTriangleMovableVertexPosition; + // CoordType movableVertex_barycentric(1 - x(n - 2), x(n - 2) / 2, x(n - 2) / 2); + // CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0) + // * movableVertex_barycentric[0] + // + global.baseTriangle.cP(1) + // * movableVertex_barycentric[1] + // + global.baseTriangle.cP(2) + // * movableVertex_barycentric[2]; + // baseTriangleMovableVertexPosition + // = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, + // vcg::math::ToRad(x(n - 1))) + // * baseTriangleMovableVertexPosition; - for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; - rotationCounter++) { - pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + // for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; + // rotationCounter++) { + // pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + // = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, + // vcg::math::ToRad(60.0 * rotationCounter)) + // * baseTriangleMovableVertexPosition; + // } + // }; + + int i = 0; + i++; + global.functions_updateReducedPatternParameter[R] = + [](const double &newR, std::shared_ptr &pReducedPatternSimulationMesh) { + const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR, + newR / 2, + newR / 2); + const CoordType hexagonBaseTriangleVertexPosition + = global.baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0] + + global.baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1] + + global.baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2]; + + for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; + rotationCounter++) { + pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, + vcg::math::ToRad(60.0 * rotationCounter)) + * hexagonBaseTriangleVertexPosition; + } + }; + + global.functions_updateReducedPatternParameter[Theta] = + [](const double &newTheta, std::shared_ptr &pReducedPatternSimulationMesh) { + const CoordType baseTriangleHexagonVertexPosition + = pReducedPatternSimulationMesh->vert[0].cP(); + + const CoordType thetaRotatedHexagonBaseTriangleVertexPosition = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, - vcg::math::ToRad(60.0 * rotationCounter)) - * baseTriangleMovableVertexPosition; - } - }; + vcg::math::ToRad(newTheta)) + * baseTriangleHexagonVertexPosition; - function_updateReducedPattern_material_E = - [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newE) { + for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; + rotationCounter++) { + pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, + vcg::math::ToRad(60.0 * rotationCounter)) + * thetaRotatedHexagonBaseTriangleVertexPosition; + } + }; + + global.functions_updateReducedPatternParameter[E] = + [](const double &newE, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE)); } }; - function_updateReducedPattern_material_A = - [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newA) { + global.functions_updateReducedPatternParameter[A] = + [](const double &newA, std::shared_ptr &pReducedPatternSimulationMesh) { const double beamWidth = std::sqrt(newA); const double beamHeight = beamWidth; for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { @@ -562,30 +599,30 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); } }; - function_updateReducedPattern_material_I = - [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI) { - for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { - Element &e = pReducedPatternSimulationMesh->elements[ei]; - e.inertia.I2 = newI; - e.inertia.I3 = newI; - } - }; - function_updateReducedPattern_material_I2 = - [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI2) { + // global.functions_updateReducedPatternParameter[I] = + // [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI) { + // for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + // Element &e = pReducedPatternSimulationMesh->elements[ei]; + // e.inertia.I2 = newI; + // e.inertia.I3 = newI; + // } + // }; + global.functions_updateReducedPatternParameter[I2] = + [](const double &newI2, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.inertia.I2 = newI2; } }; - function_updateReducedPattern_material_I3 = - [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newI3) { + global.functions_updateReducedPatternParameter[I3] = + [](const double &newI3, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.inertia.I3 = newI3; } }; - function_updateReducedPattern_material_J = - [&](std::shared_ptr &pReducedPatternSimulationMesh, const double &newJ) { + global.functions_updateReducedPatternParameter[J] = + [](const double &newJ, std::shared_ptr &pReducedPatternSimulationMesh) { for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { Element &e = pReducedPatternSimulationMesh->elements[ei]; e.inertia.J = newJ; @@ -596,38 +633,68 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() void ReducedModelOptimizer::initializeOptimizationParameters( const std::shared_ptr &mesh, const std::vector &optimizationParameters) { - global.numberOfOptimizationParameters = optimizationParameters.size(); - global.initialParameters.resize(global.numberOfOptimizationParameters); + global.numberOfOptimizationParameters = NumberOfOptimizationParameters; for (int optimizationParameterIndex = 0; optimizationParameterIndex < optimizationParameters.size(); optimizationParameterIndex++) { - const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex]; - if (xOptimizationParameter.label == "E") { - global.initialParameters(optimizationParameterIndex) = mesh->elements[0] - .material.youngsModulus; - global.updateReducedPatternFunctions_material.push_back( - function_updateReducedPattern_material_E); - } else if (xOptimizationParameter.label == "A") { - global.initialParameters(optimizationParameterIndex) = mesh->elements[0].A; - global.updateReducedPatternFunctions_material.push_back( - function_updateReducedPattern_material_A); - } else if (xOptimizationParameter.label == "I") { - global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2; - global.updateReducedPatternFunctions_material.push_back( - function_updateReducedPattern_material_I); - } else if (xOptimizationParameter.label == "I2") { - global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I2; - global.updateReducedPatternFunctions_material.push_back( - function_updateReducedPattern_material_I2); - } else if (xOptimizationParameter.label == "I3") { - global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.I3; - global.updateReducedPatternFunctions_material.push_back( - function_updateReducedPattern_material_I3); - } else if (xOptimizationParameter.label == "J") { - global.initialParameters(optimizationParameterIndex) = mesh->elements[0].inertia.J; - global.updateReducedPatternFunctions_material.push_back( - function_updateReducedPattern_material_J); + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationParameters; + optimizationParameterIndex++) { + // const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex]; + switch (optimizationParameterIndex) { + case E: + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].material.youngsModulus; + break; + case A: + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].material.youngsModulus; + break; + case I2: + global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] + .inertia.I2; + break; + case I3: + global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] + .inertia.I3; + break; + case J: + global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] + .inertia.J; + break; + case R: + global.parametersInitialValue[optimizationParameterIndex] = 0; + break; + case Theta: + global.parametersInitialValue[optimizationParameterIndex] = 0; + break; + } + + // if (xOptimizationParameter.label == "E") { + // global.updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_E); + // } else if (xOptimizationParameter.label == "A") { + // global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].A; + // global.updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_A); + // } else if (xOptimizationParameter.label == "I") { + // global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].inertia.I2; + // global.updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_I); + // } else if (xOptimizationParameter.label == "I2") { + // global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0].inertia.I2; + // global.updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_I2); + // } else if (xOptimizationParameter.label == "I3") { + // global.parametersInitialValue(optimizationParameterIndex) = mesh->elements[0].inertia.I3; + // global.updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_I3); + // } else if (xOptimizationParameter.label == "J") { + // global.parametersInitialValue(optimizationParameterIndex) = mesh->elements[0].inertia.J; + // global.updateReducedPatternFunctions_material.push_back( + // function_updateReducedPattern_material_J); + // } } //NOTE:Assuming that I2,I3 and J are be passed after E and A because the update of E and A changes I2,I3,J. // case "HexSize": @@ -640,25 +707,25 @@ void ReducedModelOptimizer::initializeOptimizationParameters( // break; } - if (global.updateReducedPatternFunctions_material.size() != 0) { - function_updateReducedPattern_material = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - for (int optimizationParameterIndex = 0; - optimizationParameterIndex - < global.updateReducedPatternFunctions_material.size(); - optimizationParameterIndex++) { - global.updateReducedPatternFunctions_material[optimizationParameterIndex]( - pReducedPatternSimulationMesh, - global.initialParameters(optimizationParameterIndex) - * x(optimizationParameterIndex)); - } - }; - } else { - function_updateReducedPattern_material = - [](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) {}; - } + // if (global.updateReducedPatternFunctions_material.size() != 0) { + // function_updateReducedPattern_material = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // for (int optimizationParameterIndex = 0; + // optimizationParameterIndex + // < global.updateReducedPatternFunctions_material.size(); + // optimizationParameterIndex++) { + // global.updateReducedPatternFunctions_material[optimizationParameterIndex]( + // pReducedPatternSimulationMesh, + // global.parametersInitialValue(optimizationParameterIndex) + // * x(optimizationParameterIndex)); + // } + // }; + // } else { + // function_updateReducedPattern_material = + // [](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) {}; + // } } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -667,17 +734,14 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( SimulationJob &simulationJobOfReducedModel) { assert(simulationJobOfReducedModel.pMesh->VN() != 0); - std::unordered_map> - reducedModelFixedVertices; - for (auto fullModelFixedVertex : - simulationJobOfFullModel.constrainedVertices) { + std::unordered_map> reducedModelFixedVertices; + for (auto fullModelFixedVertex : simulationJobOfFullModel.constrainedVertices) { reducedModelFixedVertices[fullToReducedMap.at(fullModelFixedVertex.first)] = fullModelFixedVertex.second; } std::unordered_map reducedModelNodalForces; - for (auto fullModelNodalForce : - simulationJobOfFullModel.nodalExternalForces) { + for (auto fullModelNodalForce : simulationJobOfFullModel.nodalExternalForces) { reducedModelNodalForces[fullToReducedMap.at(fullModelNodalForce.first)] = fullModelNodalForce.second; } @@ -790,22 +854,20 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, - Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) { - - assert(optimalDisplacementsOfReducedModel.rows() != 0 && - optimalDisplacementsOfReducedModel.cols() == 3); - optimalDisplacementsOfReducedModel.setZero( - optimalDisplacementsOfReducedModel.rows(), - optimalDisplacementsOfReducedModel.cols()); + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel) +{ + assert(optimalDisplacementsOfReducedModel.rows() != 0 + && optimalDisplacementsOfReducedModel.cols() == 3); + optimalDisplacementsOfReducedModel.setZero(optimalDisplacementsOfReducedModel.rows(), + optimalDisplacementsOfReducedModel.cols()); for (auto reducedFullViPair : displacementsReducedToFullMap) { const VertexIndex fullModelVi = reducedFullViPair.second; - const Vector6d fullModelViDisplacements = - fullModelResults.displacements[fullModelVi]; - optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) = - Eigen::Vector3d(fullModelViDisplacements[0], - fullModelViDisplacements[1], - fullModelViDisplacements[2]); + const Vector6d fullModelViDisplacements = fullModelResults.displacements[fullModelVi]; + optimalDisplacementsOfReducedModel.row(reducedFullViPair.first) + = Eigen::Vector3d(fullModelViDisplacements[0], + fullModelViDisplacements[1], + fullModelViDisplacements[2]); } } @@ -818,18 +880,19 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza //Value of optimal objective Y results.objectiveValue.total = optimizationResult_dlib.y; //Optimal X values - results.optimalXNameValuePairs.resize(settings.xRanges.size()); - std::vector optimalX(settings.xRanges.size()); - for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) { - if (xVariableIndex < settings.xRanges.size() - 2) { + results.optimalXNameValuePairs.resize(settings.parameterRanges.size()); + std::vector optimalX(settings.parameterRanges.size()); + for (int xVariableIndex = 0; xVariableIndex < settings.parameterRanges.size(); + xVariableIndex++) { + if (xVariableIndex < settings.parameterRanges.size() - 2) { results.optimalXNameValuePairs[xVariableIndex] - = std::make_pair(settings.xRanges[xVariableIndex].label, + = std::make_pair(settings.parameterRanges[xVariableIndex].label, optimizationResult_dlib.x(xVariableIndex) - * global.initialParameters(xVariableIndex)); + * global.parametersInitialValue[xVariableIndex]); } else { //Hex size and angle are pure values (not multipliers of the initial values) results.optimalXNameValuePairs[xVariableIndex] - = std::make_pair(settings.xRanges[xVariableIndex].label, + = std::make_pair(settings.parameterRanges[xVariableIndex].label, optimizationResult_dlib.x(xVariableIndex)); } @@ -843,6 +906,7 @@ void ReducedModelOptimizer::getResults(const dlib::function_evaluation &optimiza std::cout << std::endl; #endif + results.fullPatternYoungsModulus = youngsModulus; // Compute obj value per simulation scenario and the raw objective value // updateMeshFunction(optimalX); std::shared_ptr &pReducedPatternSimulationMesh @@ -1013,18 +1077,18 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( desiredBaseSimulationScenarioIndices); #ifdef POLYSCOPE_DEFINED - nlohmann::json json; - json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes; + nlohmann::json json; + json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes; - std::filesystem::create_directories(forceMagnitudesDirectoryPath); - std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); - jsonFile << json; + std::filesystem::create_directories(forceMagnitudesDirectoryPath); + std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); + jsonFile << json; #endif - assert(fullPatternSimulationScenarioMaxMagnitudes.size() - == desiredBaseSimulationScenarioIndices.size()); + assert(fullPatternSimulationScenarioMaxMagnitudes.size() + == desiredBaseSimulationScenarioIndices.size()); - return fullPatternSimulationScenarioMaxMagnitudes; + return fullPatternSimulationScenarioMaxMagnitudes; } void ReducedModelOptimizer::runOptimization(const Settings &settings, @@ -1036,167 +1100,66 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, global.objectiveValueHistory_iteration.reserve(settings.numberOfFunctionCalls / 100); #if POLYSCOPE_DEFINED - global.plotColors.reserve(settings.numberOfFunctionCalls); +// global.plotColors.reserve(settings.numberOfFunctionCalls); #endif double (*objF)(const dlib::matrix &) = &objective; + enum OptimizationParameterComparisonScenarioIndex { + AllVar, + noYM, + SplittedMatGeo, + SplittedYM_before, + SplittedYM_after, + NumberOfScenarios + }; + const std::vector>> scenarioParameters = + [&]() { + std::vector>> scenarioParameters( + NumberOfScenarios); + scenarioParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}}; + scenarioParameters[noYM] = {{A, I2, I3, J, R, Theta}}; + scenarioParameters[SplittedMatGeo] = {{E, A, I2, I3, J}, {R, Theta}}; + scenarioParameters[SplittedYM_before] = {{E}, {A, I2, I3, J, R, Theta}}; + scenarioParameters[SplittedYM_after] = {{A, I2, I3, J, R, Theta}, {E}}; + return scenarioParameters; + }(); + constexpr OptimizationParameterComparisonScenarioIndex scenario = AllVar; + const std::vector> scenarioParameterGroups + = scenarioParameters[scenario]; + + //TODO:Ensure that the reduced pattern mesh has the initial parameter values before starting the optimization dlib::function_evaluation optimalResult; - const auto hexAngleParameterIt = std::find_if(settings.xRanges.begin(), - settings.xRanges.end(), - [](const xRange &x) { - return x.label == "HexAngle"; - }); - const bool hasHexAngleParameter = hexAngleParameterIt != settings.xRanges.end(); - const auto hexSizeParameterIt = std::find_if(settings.xRanges.begin(), - settings.xRanges.end(), - [](const xRange &x) { - return x.label == "HexSize"; - }); - const bool hasHexSizeParameter = hexSizeParameterIt != settings.xRanges.end(); - const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter; - assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter)); - const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0; - const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters - - numberOfGeometryOptimizationParameters; - const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters - && numberOfMaterialOptimizationParameters != 0; - if (settings.splitGeometryMaterialOptimization && hasGeometryAndMaterialParameters) { - //Geometry optimization of the reduced pattern - dlib::matrix xGeometryMin(numberOfGeometryOptimizationParameters); - dlib::matrix xGeometryMax(numberOfGeometryOptimizationParameters); - // const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(), - // hexAngleParameterIt); - // const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(), - // hexSizeParameterIt); - xGeometryMin(0) = hexSizeParameterIt->min; - xGeometryMax(0) = hexSizeParameterIt->max; - xGeometryMin(1) = hexAngleParameterIt->min; - xGeometryMax(1) = hexAngleParameterIt->max; + for (const std::vector ¶meterGroup : scenarioParameterGroups) { + //Set min max values + dlib::matrix xMin(parameterGroup.size()); + dlib::matrix xMax(parameterGroup.size()); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; - //Set reduced pattern update functions to be used during optimization - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - const int optimizationFunctionCalls_geometry - = global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters - ? settings.numberOfFunctionCalls - : static_cast(settings.numberOfFunctionCalls / 3.0); - dlib::function_evaluation result_dlib_geometry - = dlib::find_min_global(objF, - xGeometryMin, - xGeometryMax, - dlib::max_function_calls(optimizationFunctionCalls_geometry), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); - //Material optimization of the reduced pattern - // std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl; - // std::cout << "opt size:" << global.minX[0] << std::endl; - dlib::function_evaluation result_dlib_material; - if (numberOfMaterialOptimizationParameters != 0) { - dlib::matrix xMaterialMin(numberOfMaterialOptimizationParameters); - dlib::matrix xMaterialMax(numberOfMaterialOptimizationParameters); - for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) { - xMaterialMin(i) = settings.xRanges[i].min; - xMaterialMax(i) = settings.xRanges[i].max; + xMin(xIndex) = settings.parameterRanges[parameterIndex].min; + xMax(xIndex) = settings.parameterRanges[parameterIndex].max; + } + //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup + function_updateReducedPattern = [&](const dlib::matrix &x, + std::shared_ptr &pMesh) { + for (int xIndex = 0; xIndex < x.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + // const double parameterInitialValue=optimizationSettings.parameterRanges[parameterIndex].initialValue; + const double parameterInitialValue = global.parametersInitialValue[parameterIndex]; + const double parameterNewValue = [&]() { + if (parameterIndex == R + || parameterIndex + == Theta) { //NOTE:Here I explore the geometry parameters linearly + return x(xIndex) + parameterInitialValue; + } + //and the material parameters exponentially(?).TODO: Check what happens if I make all linear + return x(xIndex) * parameterInitialValue; + }(); + global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, + pMesh); } - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - function_updateReducedPattern_geometry( - result_dlib_geometry.x, - pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this - pReducedPatternSimulationMesh->reset(); - }; - result_dlib_material = dlib::find_min_global(objF, - xMaterialMin, - xMaterialMax, - dlib::max_function_calls(static_cast( - 2.0 * settings.numberOfFunctionCalls - / 3.0)), - std::chrono::hours(24 * 365 * 290), - settings.solverAccuracy); - } - // constexpr bool useBOBYQA = false; - // if (useBOBYQA) { - // const size_t npt = 2 * global.numberOfOptimizationParameters; - // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; - // const double rhobeg = 0.1; - // // const double rhobeg = 10; - // const double rhoend = rhobeg * 1e-6; - // // const size_t maxFun = 10 * (x.size() ^ 2); - // const size_t bobyqa_maxFunctionCalls = 200; - // dlib::find_min_bobyqa(objF, - // result_dlib.x, - // npt, - // xMaterialMin, - // xMaterialMax, - // rhobeg, - // rhoend, - // bobyqa_maxFunctionCalls); - // } + }; - optimalResult.x.set_size(global.numberOfOptimizationParameters); - std::copy(result_dlib_material.x.begin(), - result_dlib_material.x.end(), - optimalResult.x.begin()); - // debug_x.begin()); - std::copy(result_dlib_geometry.x.begin(), - result_dlib_geometry.x.end(), - optimalResult.x.begin() + numberOfMaterialOptimizationParameters); - // std::cout << "opt x:"; - // for (const auto optx : optimalResult.x) { - // std::cout << optx << " "; - // } - // std::cout << std::endl; - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - // const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel(); - // global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre"); - // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial); - // polyscope::show(); - // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); - optimalResult.y = objective(optimalResult.x); - } else { - dlib::matrix xMin(global.numberOfOptimizationParameters); - dlib::matrix xMax(global.numberOfOptimizationParameters); - for (int i = 0; i < global.numberOfOptimizationParameters; i++) { - xMin(i) = settings.xRanges[i].min; - xMax(i) = settings.xRanges[i].max; - } - if (numberOfGeometryOptimizationParameters != 0 - && numberOfMaterialOptimizationParameters != 0) { - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - } else if (numberOfGeometryOptimizationParameters == 0) { - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - - } else { - function_updateReducedPattern = - [&](const dlib::matrix &x, - std::shared_ptr &pReducedPatternSimulationMesh) { - function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); - pReducedPatternSimulationMesh->reset(); - }; - } optimalResult = dlib::find_min_global(objF, xMin, xMax, @@ -1204,7 +1167,360 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, settings.numberOfFunctionCalls), std::chrono::hours(24 * 365 * 290), settings.solverAccuracy); + function_updateReducedPattern( + optimalResult.x, + global.reducedPatternSimulationJobs[0] + ->pMesh); //TODO: Check if its ok to update only the changed parameters } + + // optimalResult.x.set_size(global.numberOfOptimizationParameters); + // const auto hexAngleParameterIt = std::find_if(settings.parameterRanges.begin(), + // settings.parameterRanges.end(), + // [](const xRange &x) { + // return x.label == ParameterLabels::theta; + // }); + // const bool hasHexAngleParameter = hexAngleParameterIt != settings.parameterRanges.end(); + // const auto hexSizeParameterIt = std::find_if(settings.parameterRanges.begin(), + // settings.parameterRanges.end(), + // [](const xRange &x) { + // return x.label == ParameterLabels::R; + // }); + // const bool hasHexSizeParameter = hexSizeParameterIt != settings.parameterRanges.end(); + // const bool hasBothGeometricalParameters = hasHexAngleParameter && hasHexSizeParameter; + // assert(hasBothGeometricalParameters || (!hasHexAngleParameter && !hasHexSizeParameter)); + // const int numberOfGeometryOptimizationParameters = hasBothGeometricalParameters ? 2 : 0; + // const int numberOfMaterialOptimizationParameters = global.numberOfOptimizationParameters + // - numberOfGeometryOptimizationParameters; + // const bool hasGeometryAndMaterialParameters = hasBothGeometricalParameters + // && numberOfMaterialOptimizationParameters != 0; + // constexpr bool splitYMOptimization = false; + // const auto parameterIt_E = std::find_if(settings.parameterRanges.begin(), + // settings.parameterRanges.end(), + // [](const xRange &x) { + // return x.label == ParameterLabels::E; + // }); + // //E exists as a parameter + // if (splitYMOptimization && parameterIt_E == settings.parameterRanges.end()) { + // std::cerr << "YM is not a parameter." << std::endl; + // std::terminate(); + // } + + // std::unordered_map parameterLabelGlobalIndexMap; + // { + // for (int globalParameterIndex = 0; + // globalParameterIndex < global.numberOfOptimizationParameters; + // globalParameterIndex++) { + // parameterLabelGlobalIndexMap[settings.parameterRanges[globalParameterIndex].label] + // = globalParameterIndex; + // } + // } + // const int parameterIndexGlobal_E = parameterLabelGlobalIndexMap.contains(ParameterLabels::E) + // ? parameterLabelGlobalIndexMap[ParameterLabels::E] + // : -1; + // const int parameterIndexGlobal_A = parameterLabelGlobalIndexMap.contains(ParameterLabels::A) + // ? parameterLabelGlobalIndexMap[ParameterLabels::A] + // : -1; + // const int parameterIndexGlobal_I2 = parameterLabelGlobalIndexMap.contains(ParameterLabels::I2) + // ? parameterLabelGlobalIndexMap[ParameterLabels::I2] + // : -1; + // const int parameterIndexGlobal_I3 = parameterLabelGlobalIndexMap.contains(ParameterLabels::I3) + // ? parameterLabelGlobalIndexMap[ParameterLabels::I3] + // : -1; + // const int parameterIndexGlobal_J = parameterLabelGlobalIndexMap.contains(ParameterLabels::J) + // ? parameterLabelGlobalIndexMap[ParameterLabels::J] + // : -1; + // const int parameterIndexGlobal_R = parameterLabelGlobalIndexMap.contains(ParameterLabels::R) + // ? parameterLabelGlobalIndexMap[ParameterLabels::R] + // : -1; + // const int parameterIndexGlobal_theta = parameterLabelGlobalIndexMap.contains( + // ParameterLabels::theta) + // ? parameterLabelGlobalIndexMap[ParameterLabels::theta] + // : -1; + + // if (splitYMOptimization) { + // std::unordered_map parameterLabelLocalIndexMap; + // //Optimize all variables besides the youngs modulus + // const int parameterIndex_E = std::distance(settings.parameterRanges.begin(), parameterIt_E); + // dlib::matrix xMin(global.numberOfOptimizationParameters - 1); + // dlib::matrix xMax(global.numberOfOptimizationParameters - 1); + // { + // int xIndex = 0; + // for (int globalParameterIndex = 0; + // globalParameterIndex < global.numberOfOptimizationParameters; + // globalParameterIndex++) { + // if (globalParameterIndex == parameterIndexGlobal_E) { + // continue; + // } + // parameterLabelLocalIndexMap[settings.parameterRanges[globalParameterIndex].label] + // = xIndex; + // xMin(xIndex) = settings.parameterRanges[globalParameterIndex].min; + // xMax(xIndex) = settings.parameterRanges[globalParameterIndex].max; + // xIndex++; + // } + // } + // const int parameterLocalIndex_A = parameterLabelLocalIndexMap[ParameterLabels::A]; + // const int parameterLocalIndex_I2 = parameterLabelLocalIndexMap[ParameterLabels::I2]; + // const int parameterLocalIndex_I3 = parameterLabelLocalIndexMap[ParameterLabels::I3]; + // const int parameterLocalIndex_J = parameterLabelLocalIndexMap[ParameterLabels::J]; + // const int parameterLocalIndex_R = parameterLabelLocalIndexMap[ParameterLabels::R]; + // const int parameterLocalIndex_theta = parameterLabelLocalIndexMap[ParameterLabels::theta]; + // constexpr std::array<7, bool> optimizeAllBesidesE(false); + + // function_updateReducedPattern = + // [](const dlib::matrix &x, + // const std::vector &initialValues, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_material_A(pReducedPatternSimulationMesh, + // initialValues[parameterIndexGlobal_A] + // * x(parameterLocalIndex_A)); + // function_updateReducedPattern_material_I2(pReducedPatternSimulationMesh, + // initialValues[parameterIndexGlobal_I2] + // * x(parameterLocalIndex_I2)); + // function_updateReducedPattern_material_I3(pReducedPatternSimulationMesh, + // initialValues[parameterIndexGlobal_I3] + // * x(parameterLocalIndex_I3)); + // function_updateReducedPattern_material_J(pReducedPatternSimulationMesh, + // initialValues[parameterIndexGlobal_J] + // * x(parameterLocalIndex_J)); + // function_updateReducedPattern_geometry({x(parameterLocalIndex_R), + // x(parameterLocalIndex_theta)}, + // pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + + // const dlib::function_evaluation result_dlib_noYM + // = dlib::find_min_global(objF, + // xMin, + // xMax, + // dlib::max_function_calls(settings.numberOfFunctionCalls), + // std::chrono::hours(24 * 365 * 290), + // settings.solverAccuracy); + // function_updateReducedPattern(result_dlib_noYM.x, + // global.reducedPatternSimulationJobs[0]->pMesh); + // //Optimize for the youngs modulus only + // dlib::matrix xEMin(1); + // dlib::matrix xEMax(1); + // xEMin(0) = parameterIt_E->min; + // xEMax(0) = parameterIt_E->max; + + // function_updateReducedPattern = [&](const dlib::matrix &x, + // const std::vector &initialValues, + // std::shared_ptr + // &pReducedPatternSimulationMesh) { + // // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + // function_updateReducedPattern_material_E(pReducedPatternSimulationMesh, + // x(0) * initialValues[parameterIndexGlobal_E]); + // function_updateReducedPattern_geometry( + // {result_dlib_noYM.x(parameterLocalIndex_R), + // result_dlib_noYM.x(parameterLocalIndex_theta)}, + // pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this + // pReducedPatternSimulationMesh->reset(); + // }; + // double optimalE = 1; + // double (*objective_singleVariable)(const double &) = &objective; + // const double optimalObjectiveValue = dlib::find_min_single_variable(objective_singleVariable, + // optimalE, + // parameterIt_E->min, + // parameterIt_E->max, + // 1e-2, + // 100); + // //Gather results of the two optimizations + // for (int parameterIndex = 0; parameterIndex < numberOfMaterialOptimizationParameters; + // parameterIndex++) { + // if (parameterIndex == parameterIndex_E) { + // optimalResult.x(parameterIndex) = optimalE; + // continue; + // } + // assert(parameterLabelLocalIndexMap.contains( + // settings.parameterRanges[parameterIndex].label)); + // const int firstOptimizationParameterIndex + // = parameterLabelLocalIndexMap[settings.parameterRanges[parameterIndex].label]; + // optimalResult.x(parameterIndex) = result_dlib_noYM.x(firstOptimizationParameterIndex); + // } + // function_updateReducedPattern = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + + // } else if (settings.splitGeometryMaterialOptimization && hasGeometryAndMaterialParameters) { + // //Geometry optimization of the reduced pattern + // dlib::matrix xGeometryMin(numberOfGeometryOptimizationParameters); + // dlib::matrix xGeometryMax(numberOfGeometryOptimizationParameters); + // // const int hexAngleParameterIndex = std::distance(settings.xRanges.begin(), + // // hexAngleParameterIt); + // // const int hexSizeParameterIndex = std::distance(settings.xRanges.begin(), + // // hexSizeParameterIt); + // xGeometryMin(0) = hexSizeParameterIt->min; + // xGeometryMax(0) = hexSizeParameterIt->max; + // xGeometryMin(1) = hexAngleParameterIt->min; + // xGeometryMax(1) = hexAngleParameterIt->max; + + // //Set reduced pattern update functions to be used during optimization + // function_updateReducedPattern = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + // const int optimizationFunctionCalls_geometry + // = global.numberOfOptimizationParameters == numberOfGeometryOptimizationParameters + // ? settings.numberOfFunctionCalls + // : static_cast(settings.numberOfFunctionCalls / 3.0); + // dlib::function_evaluation result_dlib_geometry + // = dlib::find_min_global(objF, + // xGeometryMin, + // xGeometryMax, + // dlib::max_function_calls(optimizationFunctionCalls_geometry), + // std::chrono::hours(24 * 365 * 290), + // settings.solverAccuracy); + // function_updateReducedPattern(result_dlib_geometry.x, + // global.reducedPatternSimulationJobs[0]->pMesh); + // //Material optimization of the reduced pattern + // // std::cout << "opt size:" << result_dlib_geometry.x(0) << std::endl; + // // std::cout << "opt size:" << global.minX[0] << std::endl; + // dlib::function_evaluation result_dlib_material; + // if (numberOfMaterialOptimizationParameters != 0) { + // dlib::matrix xMaterialMin(numberOfMaterialOptimizationParameters); + // dlib::matrix xMaterialMax(numberOfMaterialOptimizationParameters); + // for (int i = 0; i < numberOfMaterialOptimizationParameters; i++) { + // xMaterialMin(i) = settings.parameterRanges[i].min; + // xMaterialMax(i) = settings.parameterRanges[i].max; + // } + // function_updateReducedPattern = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + // function_updateReducedPattern_geometry( + // result_dlib_geometry.x, + // pReducedPatternSimulationMesh); //FIXME: I shouldn't need to call this + // pReducedPatternSimulationMesh->reset(); + // }; + // result_dlib_material = dlib::find_min_global(objF, + // xMaterialMin, + // xMaterialMax, + // dlib::max_function_calls(static_cast( + // 2.0 * settings.numberOfFunctionCalls + // / 3.0)), + // std::chrono::hours(24 * 365 * 290), + // settings.solverAccuracy); + // } + // // constexpr bool useBOBYQA = false; + // // if (useBOBYQA) { + // // const size_t npt = 2 * global.numberOfOptimizationParameters; + // // // ((n + 2) + ((n + 1) * (n + 2) / 2)) / 2; + // // const double rhobeg = 0.1; + // // // const double rhobeg = 10; + // // const double rhoend = rhobeg * 1e-6; + // // // const size_t maxFun = 10 * (x.size() ^ 2); + // // const size_t bobyqa_maxFunctionCalls = 200; + // // dlib::find_min_bobyqa(objF, + // // result_dlib.x, + // // npt, + // // xMaterialMin, + // // xMaterialMax, + // // rhobeg, + // // rhoend, + // // bobyqa_maxFunctionCalls); + // // } + + // std::copy(result_dlib_material.x.begin(), + // result_dlib_material.x.end(), + // optimalResult.x.begin()); + // // debug_x.begin()); + // std::copy(result_dlib_geometry.x.begin(), + // result_dlib_geometry.x.end(), + // optimalResult.x.begin() + numberOfMaterialOptimizationParameters); + // // std::cout << "opt x:"; + // // for (const auto optx : optimalResult.x) { + // // std::cout << optx << " "; + // // } + // // std::cout << std::endl; + // function_updateReducedPattern = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + // // const auto meshLabel = global.reducedPatternSimulationJobs[0]->pMesh->getLabel(); + // // global.reducedPatternSimulationJobs[0]->pMesh->setLabel("pre"); + // // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(Colors::reducedInitial); + // // polyscope::show(); + // // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); + // } else { + // dlib::matrix xMin(global.numberOfOptimizationParameters); + // dlib::matrix xMax(global.numberOfOptimizationParameters); + // for (int i = 0; i < global.numberOfOptimizationParameters; i++) { + // xMin(i) = settings.parameterRanges[i].min; + // xMax(i) = settings.parameterRanges[i].max; + // } + // if (numberOfGeometryOptimizationParameters != 0 + // && numberOfMaterialOptimizationParameters != 0) { + // function_updateReducedPattern = [&](const dlib::matrix &x, + // std::shared_ptr + // &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // // function_updateReducedPattern_material_E(pReducedPatternSimulationMesh, + // // global.initialParameters( + // // parameterIndexGlobal_E) + // // * x(parameterIndexGlobal_E)); + // // function_updateReducedPattern_material_A(pReducedPatternSimulationMesh, + // // global.initialParameters( + // // parameterIndexGlobal_A) + // // * x(parameterIndexGlobal_A)); + // // function_updateReducedPattern_material_I2(pReducedPatternSimulationMesh, + // // global.initialParameters( + // // parameterIndexGlobal_I2) + // // * x(parameterIndexGlobal_I2)); + // // function_updateReducedPattern_material_I3(pReducedPatternSimulationMesh, + // // global.initialParameters( + // // parameterIndexGlobal_I3) + // // * x(parameterIndexGlobal_I3)); + // // function_updateReducedPattern_material_J(pReducedPatternSimulationMesh, + // // global.initialParameters( + // // parameterIndexGlobal_J) + // // * x(parameterIndexGlobal_J)); + // // function_updateReducedPattern_geometry({x(parameterIndexGlobal_R), + // // x(parameterIndexGlobal_theta)}, + // // pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + // } /* else if (numberOfGeometryOptimizationParameters == 0) { + // function_updateReducedPattern = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_material(x, pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + + // } else { + // function_updateReducedPattern = + // [&](const dlib::matrix &x, + // std::shared_ptr &pReducedPatternSimulationMesh) { + // function_updateReducedPattern_geometry(x, pReducedPatternSimulationMesh); + // pReducedPatternSimulationMesh->reset(); + // }; + // }*/ + // optimalResult = dlib::find_min_global(objF, + // xMin, + // xMax, + // dlib::max_function_calls( + // settings.numberOfFunctionCalls), + // std::chrono::hours(24 * 365 * 290), + // settings.solverAccuracy); + // } + // function_updateReducedPattern(optimalResult.x, global.reducedPatternSimulationJobs[0]->pMesh); + // const auto debug_currObj = objective(optimalResult.x); + // if (std::abs(debug_currObj - optimalResult.y) > 1e-5) { + // std::cout << "wrong optimal result.Previous was:" << optimalResult.y + // << " .Current is:" << debug_currObj << std::endl; + // } + // optimalResult.y = debug_currObj; getResults(optimalResult, settings, results); } @@ -1225,7 +1541,6 @@ void ReducedModelOptimizer::constructAxialSimulationScenario( job.constrainedVertices[viPairIt->second] = std::unordered_set{0, 1, 2, 3, 4, 5}; } } - } void ReducedModelOptimizer::constructShearSimulationScenario( @@ -1581,75 +1896,74 @@ std::vector> ReducedModelOptimizer::createFullPat return scenarios; } -void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { - - // Compute the sum of the displacement norms - std::vector fullPatternTranslationalDisplacementNormSum( - totalNumberOfSimulationScenarios); - std::vector fullPatternAngularDistance(totalNumberOfSimulationScenarios); - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - double translationalDisplacementNormSum = 0; - for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { - const int fullPatternVi = interfaceViPair.second; - //If the full pattern vertex is translationally constrained dont take it into account - if (global.fullPatternSimulationJobs[simulationScenarioIndex] - ->constrainedVertices.contains(fullPatternVi)) { - const std::unordered_set constrainedDof - = global.fullPatternSimulationJobs[simulationScenarioIndex] - ->constrainedVertices.at(fullPatternVi); - if (constrainedDof.contains(0) && constrainedDof.contains(1) - && constrainedDof.contains(2)) { - continue; - } +void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() +{ + // Compute the sum of the displacement norms + std::vector fullPatternTranslationalDisplacementNormSum( + totalNumberOfSimulationScenarios); + std::vector fullPatternAngularDistance(totalNumberOfSimulationScenarios); + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + double translationalDisplacementNormSum = 0; + for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + const int fullPatternVi = interfaceViPair.second; + //If the full pattern vertex is translationally constrained dont take it into account + if (global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.contains(fullPatternVi)) { + const std::unordered_set constrainedDof + = global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.at(fullPatternVi); + if (constrainedDof.contains(0) && constrainedDof.contains(1) + && constrainedDof.contains(2)) { + continue; } - - const Vector6d &vertexDisplacement = global - .fullPatternResults[simulationScenarioIndex] - .displacements[fullPatternVi]; - translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); - } - double angularDistanceSum = 0; - for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { - const int fullPatternVi = interfaceViPair.second; - //If the full pattern vertex is rotationally constrained dont take it into account - if (global.fullPatternSimulationJobs[simulationScenarioIndex] - ->constrainedVertices.contains(fullPatternVi)) { - const std::unordered_set constrainedDof - = global.fullPatternSimulationJobs[simulationScenarioIndex] - ->constrainedVertices.at(fullPatternVi); - if (constrainedDof.contains(3) && constrainedDof.contains(5) - && constrainedDof.contains(4)) { - continue; - } - } - angularDistanceSum += global.fullPatternResults[simulationScenarioIndex] - .rotationalDisplacementQuaternion[fullPatternVi] - .angularDistance(Eigen::Quaterniond::Identity()); } - fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] - = translationalDisplacementNormSum; - fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; + const Vector6d &vertexDisplacement = global.fullPatternResults[simulationScenarioIndex] + .displacements[fullPatternVi]; + translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); + } + double angularDistanceSum = 0; + for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + const int fullPatternVi = interfaceViPair.second; + //If the full pattern vertex is rotationally constrained dont take it into account + if (global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.contains(fullPatternVi)) { + const std::unordered_set constrainedDof + = global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.at(fullPatternVi); + if (constrainedDof.contains(3) && constrainedDof.contains(5) + && constrainedDof.contains(4)) { + continue; + } + } + angularDistanceSum += global.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion[fullPatternVi] + .angularDistance(Eigen::Quaterniond::Identity()); } - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - if (global.optimizationSettings.normalizationStrategy == - Settings::NormalizationStrategy::Epsilon) { - const double epsilon_translationalDisplacement - = global.optimizationSettings.translationNormalizationParameter; - global.translationalDisplacementNormalizationValues[simulationScenarioIndex] - = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], - epsilon_translationalDisplacement); - const double epsilon_rotationalDisplacement = global.optimizationSettings - .rotationNormalizationParameter; - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] - = std::max(fullPatternAngularDistance[simulationScenarioIndex], - epsilon_rotationalDisplacement); - } else { - global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; - global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; - } + fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] + = translationalDisplacementNormSum; + fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; + } + + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + if (global.optimizationSettings.normalizationStrategy + == Settings::NormalizationStrategy::Epsilon) { + const double epsilon_translationalDisplacement = global.optimizationSettings + .translationNormalizationParameter; + global.translationalDisplacementNormalizationValues[simulationScenarioIndex] + = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], + epsilon_translationalDisplacement); + const double epsilon_rotationalDisplacement = global.optimizationSettings + .rotationNormalizationParameter; + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] + = std::max(fullPatternAngularDistance[simulationScenarioIndex], + epsilon_rotationalDisplacement); + } else { + global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; } + } } void ReducedModelOptimizer::optimize( @@ -1692,8 +2006,8 @@ void ReducedModelOptimizer::optimize( results.baseTriangle = global.baseTriangle; DRMSimulationModel::Settings simulationSettings; -// simulationSettings.maxDRMIterations = 200000; -// simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; + simulationSettings.maxDRMIterations = 200000; + simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; simulationSettings.viscousDampingFactor = 5e-3; simulationSettings.useKineticDamping = true; @@ -1715,8 +2029,11 @@ void ReducedModelOptimizer::optimize( const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, - simulationSettings); + // SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, + // simulationSettings); + LinearSimulationModel linearSimulator; + SimulationResults fullPatternResults = linearSimulator.executeSimulation( + pFullPatternSimulationJob); // if (!fullPatternResults.converged) { // DRMSimulationModel::Settings simulationSettings_secondRound; // simulationSettings_secondRound.viscousDampingFactor = 2e-3; @@ -1797,10 +2114,10 @@ void ReducedModelOptimizer::optimize( global.fullPatternSimulationJobs[0]->pMesh->unregister(); } #endif - // if (global.optimizationSettings.normalizationStrategy - // != Settings::NormalizationStrategy::NonNormalized) { - computeObjectiveValueNormalizationFactors(); - // } + if (global.optimizationSettings.normalizationStrategy + == Settings::NormalizationStrategy::Epsilon) { + computeObjectiveValueNormalizationFactors(); + } #ifdef POLYSCOPE_DEFINED std::cout << "Running reduced model optimization" << std::endl; #endif diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 454d59f..8bcd6a2 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -36,8 +36,22 @@ class ReducedModelOptimizer constructBaseScenarioFunctions; std::vector scenarioIsSymmetrical; int fullPatternNumberOfEdges; + constexpr static double youngsModulus{1 * 1e9}; public: + constexpr static struct ParameterLabels + { + inline const static std::string E = {"E"}; + inline const static std::string A = {"A"}; + inline const static std::string I2 = {"I2"}; + inline const static std::string I3 = {"I3"}; + inline const static std::string J = {"J"}; + inline const static std::string theta = {"Theta"}; + inline const static std::string R = {"R"}; + } parameterLabels; + + enum OptimizationParameterIndex { E, A, I2, I3, J, R, Theta, NumberOfOptimizationParameters }; + constexpr static std::array simulationScenariosResolution = {11, 11, 20, 20, 20}; // constexpr static std::array simulationScenariosResolution = {3, 3, 3, 3, 3}; inline static int totalNumberOfSimulationScenarios @@ -74,15 +88,6 @@ public: static void runSimulation(const std::string &filename, std::vector &x); - 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); - static std::vector> createFullPatternSimulationJobs( const std::shared_ptr &pMesh, const std::unordered_map &fullPatternOppositeInterfaceViMap); @@ -179,31 +184,31 @@ public: &oppositeInterfaceViPairs, SimulationJob &job); static std::function &x, - std::shared_ptr &m)> + std::shared_ptr &pReducedPatternSimulationMesh)> function_updateReducedPattern; static std::function &x, - std::shared_ptr &m)> + std::shared_ptr &pReducedPatternSimulationMesh)> function_updateReducedPattern_material; - static std::function &pReducedPatternSimulationMesh, - const double &newE)> + static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_E; - static std::function &pReducedPatternSimulationMesh, - const double &newA)> + static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_A; - static std::function &pReducedPatternSimulationMesh, - const double &newI)> + static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_I; - static std::function &pReducedPatternSimulationMesh, - const double &newI2)> + static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_I2; - static std::function &pReducedPatternSimulationMesh, - const double &newI3)> + static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_I3; - static std::function &pReducedPatternSimulationMesh, - const double &newJ)> + static std::function &pReducedPatternSimulationMesh)> function_updateReducedPattern_material_J; static std::function &x, - std::shared_ptr &m)> + std::shared_ptr &pReducedPatternSimulationMesh)> function_updateReducedPattern_geometry; private: @@ -223,7 +228,6 @@ private: const std::shared_ptr &mesh, const std::vector &optimizationParamters); - static double objective(long n, const double *x); DRMSimulationModel simulator; void computeObjectiveValueNormalizationFactors(); static void getResults(const dlib::function_evaluation &optimizationResult_dlib, @@ -242,30 +246,32 @@ private: &desiredBaseSimulationScenarioIndices); static double objective(const dlib::matrix &x); static void initializeUpdateReducedPatternFunctions(); + static double objective(const double &xValue); }; -inline std::function &x, std::shared_ptr &m)> - ReducedModelOptimizer::function_updateReducedPattern; inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern_material; -inline std::function &pReducedPatternSimulationMesh, - const double &newE)> +inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_E; -inline std::function &pReducedPatternSimulationMesh, - const double &newA)> +inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_A; -inline std::function &pReducedPatternSimulationMesh, - const double &newI)> +inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_I; -inline std::function &pReducedPatternSimulationMesh, - const double &newI2)> +inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_I2; -inline std::function &pReducedPatternSimulationMesh, - const double &newI3)> +inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_I3; -inline std::function &pReducedPatternSimulationMesh, - const double &newJ)> +inline std::function &pReducedPatternSimulationMesh)> ReducedModelOptimizer::function_updateReducedPattern_material_J; inline std::function &x, std::shared_ptr &m)> ReducedModelOptimizer::function_updateReducedPattern_geometry; +inline std::function &x, + std::shared_ptr &m)> + ReducedModelOptimizer::function_updateReducedPattern; #endif // REDUCEDMODELOPTIMIZER_HPP