diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp new file mode 100644 index 0000000..08267f5 --- /dev/null +++ b/reducedmodeloptimizer.cpp @@ -0,0 +1,2073 @@ +#include "reducedmodeloptimizer.hpp" +#include "hexagonremesher.hpp" +#include "linearsimulationmodel.hpp" +#include "simulationhistoryplotter.hpp" +#include "trianglepatterngeometry.hpp" +#include "trianglepattterntopology.hpp" +#include +#include +#include +#include +#include +#include +#include + +//#define USE_SCENARIO_WEIGHTS + +using namespace ReducedModelOptimization; + +struct GlobalOptimizationVariables +{ + std::vector fullPatternResults; + std::vector translationalDisplacementNormalizationValues; + std::vector rotationalDisplacementNormalizationValues; + std::vector> fullPatternSimulationJobs; + std::vector> reducedPatternSimulationJobs; + std::unordered_map reducedToFullInterfaceViMap; + std::vector> + fullPatternInterfaceViPairs; + matplot::line_handle gPlotHandle; + std::vector objectiveValueHistory_iteration; + std::vector objectiveValueHistory; + std::vector plotColors; + std::array + parametersInitialValue; + std::array + optimizationInitialValue; + std::vector simulationScenarioIndices; + double minY{DBL_MAX}; + std::vector minX; + int numOfSimulationCrashes{false}; + int numberOfFunctionCalls{0}; + ReducedModelOptimization::Settings optimizationSettings; + vcg::Triangle3 baseTriangle; + //Variables for finding the full pattern simulation forces + std::shared_ptr pFullPatternSimulationMesh; + std::function> &, + SimulationJob &)> + constructScenarioFunction; + FullPatternVertexIndex interfaceViForComputingScenarioError; + double desiredMaxDisplacementValue; + double desiredMaxRotationAngle; + std::string currentScenarioName; + std::array &pReducedPatternSimulationMesh)>, + 7> + functions_updateReducedPatternParameter; + std::vector xMin; + std::vector xMax; + std::vector scenarioWeights; + std::vector objectiveWeights; +} global; + +double ReducedModelOptimizer::computeDisplacementError( + const std::vector &fullPatternDisplacements, + const std::vector &reducedPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor) +{ + const double rawError = computeRawTranslationalError(fullPatternDisplacements, + reducedPatternDisplacements, + reducedToFullInterfaceViMap); + // std::cout << "raw trans error:" << rawError << std::endl; + // std::cout << "raw trans error:" << normalizationFactor << std::endl; + return rawError / normalizationFactor; +} + +double ReducedModelOptimizer::computeRawTranslationalError( + const std::vector &fullPatternDisplacements, + const std::vector &reducedPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap) +{ + double error = 0; + for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { + const VertexIndex reducedModelVi = reducedFullViPair.first; + const VertexIndex fullModelVi = reducedFullViPair.second; + const Eigen::Vector3d fullPatternVertexDiplacement = fullPatternDisplacements[fullModelVi] + .getTranslation(); + const Eigen::Vector3d reducedPatternVertexDiplacement + = reducedPatternDisplacements[reducedModelVi].getTranslation(); + const double vertexError = (fullPatternVertexDiplacement - reducedPatternVertexDiplacement) + .norm(); + error += vertexError; + } + + return error; +} + +double ReducedModelOptimizer::computeRawRotationalError( + const std::vector> &rotatedQuaternion_fullPattern, + const std::vector> &rotatedQuaternion_reducedPattern, + const std::unordered_map + &reducedToFullInterfaceViMap) +{ + double rawRotationalError = 0; + for (const auto &reducedToFullInterfaceViPair : reducedToFullInterfaceViMap) { + const double vertexRotationalError = std::abs( + rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance( + rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first])); + rawRotationalError += vertexRotationalError; + } + + return rawRotationalError; +} + +double ReducedModelOptimizer::computeRotationalError( + const std::vector> &rotatedQuaternion_fullPattern, + const std::vector> &rotatedQuaternion_reducedPattern, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor) +{ + const double rawRotationalError = computeRawRotationalError(rotatedQuaternion_fullPattern, + rotatedQuaternion_reducedPattern, + reducedToFullInterfaceViMap); + return rawRotationalError / normalizationFactor; +} + +double ReducedModelOptimizer::computeError( + const SimulationResults &simulationResults_fullPattern, + const SimulationResults &simulationResults_reducedPattern, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor_translationalDisplacement, + const double &normalizationFactor_rotationalDisplacement, + const double &scenarioWeight, + const Settings::ObjectiveWeights &objectiveWeights) +{ + const double translationalError + = computeDisplacementError(simulationResults_fullPattern.displacements, + simulationResults_reducedPattern.displacements, + reducedToFullInterfaceViMap, + normalizationFactor_translationalDisplacement); + // const double translationalError + // = computeRawTranslationalError(simulationResults_fullPattern.displacements, + // simulationResults_reducedPattern.displacements, + // reducedToFullInterfaceViMap); + + // std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl; + // std::cout << "trans error:" << translationalError << std::endl; + const double rotationalError + = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, + simulationResults_reducedPattern.rotationalDisplacementQuaternion, + reducedToFullInterfaceViMap, + normalizationFactor_rotationalDisplacement); + // const double rotationalError + // = computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, + // simulationResults_reducedPattern.rotationalDisplacementQuaternion, + // reducedToFullInterfaceViMap); + // std::cout << "rot error:" << rotationalError<< std::endl; + const double simulationError = objectiveWeights.translational * translationalError + + objectiveWeights.rotational * rotationalError; + assert(!std::isnan(simulationError)); + return simulationError * simulationError * scenarioWeight * scenarioWeight; +} + +//double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, +// double innerHexagonSize, +// double innerHexagonRotationAngle) { +// std::vector x{E,A,J,I2,I3, innerHexagonSize, innerHexagonRotationAngle}; +// return ReducedModelOptimizer::objective(x.size(), x.data()); +//} + +#ifdef USE_ENSMALLEN +struct EnsmallenOptimizationObjective +{ + static double Evaluate(const arma::mat &x_arma) + { + std::vector x(x_arma.begin(), x_arma.end()); + for (int xi = 0; xi < x.size(); xi++) { + if (x[xi] > global.xMax[xi]) { + //#ifdef POLYSCOPE_DEFINED + // std::cout << "Out of range" << std::endl; + //#endif + return std::numeric_limits::max(); + // return 1e10; + // x[xi] = global.xMax[xi]; + } else if (x[xi] < global.xMin[xi]) { + //#ifdef POLYSCOPE_DEFINED + // std::cout << "Out of range" << std::endl; + //#endif + return std::numeric_limits::max(); + // return 1e10; + // x[xi] = global.xMin[xi]; + } + } + return ReducedModelOptimizer::objective(x); + } +}; +#endif + +#ifdef DLIB_DEFINED +double ReducedModelOptimizer::objective(const dlib::matrix &x) +{ + return objective(std::vector(x.begin(), x.end())); +} +#endif + +double ReducedModelOptimizer::objective(const double &xValue) +{ + return objective(std::vector({xValue})); +} + +double ReducedModelOptimizer::objective(const std::vector &x) +{ + // std::cout.precision(17); + // for (int i = 0; i < x.size(); i++) { + // std::cout << x[i] << " "; + // } + // std::cout << 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 << " " + // << e.firstBendingConstFactor << " " << + // e.secondBendingConstFactor + // << std::endl; + // const int n = x.size(); + std::shared_ptr &pReducedPatternSimulationMesh + = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; + function_updateReducedPattern(x, pReducedPatternSimulationMesh); + // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); + // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); + // polyscope::show(); + // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); + + // run simulations + std::vector simulationErrorsPerScenario(totalNumberOfSimulationScenarios, 0); + LinearSimulationModel simulator; + simulator.setStructure(pReducedPatternSimulationMesh); + // simulator.initialize(); + // FormFinder simulator; + + std::for_each( + std::execution::par_unseq, + global.simulationScenarioIndices.begin(), + global.simulationScenarioIndices.end(), + [&](const int &simulationScenarioIndex) { + // for (const int simulationScenarioIndex : global.simulationScenarioIndices) { + const std::shared_ptr &reducedJob + = global.reducedPatternSimulationJobs[simulationScenarioIndex]; + //#ifdef POLYSCOPE_DEFINED + // std::cout << reducedJob->getLabel() << ":" << std::endl; + //#endif + SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); + // std::string filename; + if (!reducedModelResults.converged) { + simulationErrorsPerScenario[simulationScenarioIndex] + = std::numeric_limits::max(); + global.numOfSimulationCrashes++; +#ifdef POLYSCOPE_DEFINED + std::cout << "Failed simulation" << std::endl; +#endif + } else { +// const double simulationScenarioError_PE = std::abs( +// (reducedModelResults.internalPotentialEnergy +// - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy) +// / global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); +// else { +#ifdef POLYSCOPE_DEFINED +#ifdef USE_SCENARIO_WEIGHTS + const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex]; +#else + const double scenarioWeight = 1; +#endif +#else + const double scenarioWeight = 1; +#endif + const double simulationScenarioError_displacements = computeError( + global.fullPatternResults[simulationScenarioIndex], + reducedModelResults, + global.reducedToFullInterfaceViMap, + global.translationalDisplacementNormalizationValues[simulationScenarioIndex], + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + scenarioWeight, + global.objectiveWeights[simulationScenarioIndex]); + + simulationErrorsPerScenario[simulationScenarioIndex] + = simulationScenarioError_displacements /*+ simulationScenarioError_PE*/; + // } + // #ifdef POLYSCOPE_DEFINED + // reducedJob->pMesh->registerForDrawing(Colors::reducedInitial); + // reducedModelResults.registerForDrawing(Colors::reducedDeformed); + // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing( + // Colors::fullDeformed); + // polyscope::show(); + // reducedModelResults.unregister(); + // global.pFullPatternSimulationMesh->unregister(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); + // #endif + } + }); + // double totalError = std::accumulate(simulationErrorsPerScenario.begin(), + // simulationErrorsPerScenario.end(), + // 0); + const double totalError = std::reduce(std::execution::par_unseq, + simulationErrorsPerScenario.begin(), + simulationErrorsPerScenario.end()); + // std::cout << totalError << std::endl; + // 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); +#ifdef POLYSCOPE_DEFINED + std::cout << "New best:" << totalError << std::endl; + // // global.minX.assign(x.begin(), x.begin() + n); + std::cout.precision(17); + for (int i = 0; i < x.size(); i++) { + std::cout << x[i] << " "; + } + std::cout << std::endl; +#endif + // global.objectiveValueHistoryY.push_back(std::log(totalError)); + // global.objectiveValueHistoryX.push_back(global.numberOfFunctionCalls); + // global.plotColors.push_back(0.1); + // auto xPlot = matplot::linspace(0, + // global.objectiveValueHistoryY.size(), + // global.objectiveValueHistoryY.size()); + // global.gPlotHandle = matplot::scatter(global.objectiveValueHistoryX, + // global.objectiveValueHistoryY, + // 4, + // global.plotColors); + // matplot::show(); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Objective value", + // global.objectiveValueHistoryY); + } +#ifdef POLYSCOPE_DEFINED + if (global.optimizationSettings.numberOfFunctionCalls >= 100 + && global.numberOfFunctionCalls % (global.optimizationSettings.numberOfFunctionCalls / 100) + == 0) { + std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; + } +#endif + // compute error and return it + return totalError; +} + +void ReducedModelOptimizer::createSimulationMeshes( + PatternGeometry &fullModel, + PatternGeometry &reducedModel, + std::shared_ptr &pFullPatternSimulationMesh, + std::shared_ptr &pReducedPatternSimulationMesh) +{ + if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) { + std::cerr << "Error: A rectangular cross section is expected." << std::endl; + std::terminate(); + } + // Full pattern + pFullPatternSimulationMesh = std::make_shared(fullModel); + pFullPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); + pFullPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); + pFullPatternSimulationMesh->reset(); + + // Reduced pattern + pReducedPatternSimulationMesh = std::make_shared(reducedModel); + pReducedPatternSimulationMesh->setBeamCrossSection(CrossSectionType{0.002, 0.002}); + pReducedPatternSimulationMesh->setBeamMaterial(0.3, youngsModulus); + pReducedPatternSimulationMesh->reset(); +} + +void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel, + PatternGeometry &reducedModel) +{ + ReducedModelOptimizer::createSimulationMeshes(fullModel, + reducedModel, + m_pFullPatternSimulationMesh, + m_pReducedPatternSimulationMesh); +} + +void ReducedModelOptimizer::computeMaps( + const std::unordered_map> &slotToNode, + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + std::unordered_map + &reducedToFullInterfaceViMap, + std::unordered_map + &fullToReducedInterfaceViMap, + std::vector> + &fullPatternOppositeInterfaceViPairs) +{ + // 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); + // 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()); + + 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 fullModelBaseTriangleVN = fullPattern.VN(); + fullPattern.createFan(); + const size_t duplicateVerticesPerFanPair = fullModelBaseTriangleVN - fullPattern.VN() / 6; + const size_t fullPatternInterfaceVertexOffset = fullModelBaseTriangleVN + - duplicateVerticesPerFanPair; + // std::cout << "Dups in fan pair:" << duplicateVerticesPerFanPair << + // std::endl; + + // Save excluded edges TODO:this changes the global object. Should this be a + // function parameter? + // global.reducedPatternExludedEdges.clear(); + // const size_t reducedBaseTriangleNumberOfEdges = reducedPattern.EN(); + // for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { + // for (const size_t ei : reducedModelExcludedEdges) { + // global.reducedPatternExludedEdges.insert( + // fanIndex * reducedBaseTriangleNumberOfEdges + ei); + // } + // } + + // Construct reduced->full and full->reduced interface vi map + reducedToFullInterfaceViMap.clear(); + vcg::tri::Allocator::PointerUpdater pu_reducedModel; + reducedPattern.deleteDanglingVertices(pu_reducedModel); + const size_t reducedModelBaseTriangleInterfaceVi = pu_reducedModel.remap[baseTriangleInterfaceVi]; + const size_t reducedModelInterfaceVertexOffset + = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/; + Results::applyOptimizationResults_innerHexagon(initialHexagonSize, + 0, + global.baseTriangle, + 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; + } + fullToReducedInterfaceViMap.clear(); + constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap); + + // Create opposite vertex map + fullPatternOppositeInterfaceViPairs.clear(); + fullPatternOppositeInterfaceViPairs.reserve(fanSize / 2); + // for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { + for (int fanIndex = 0; fanIndex < fanSize / 2; fanIndex++) { + const size_t vi0 = fullModelBaseTriangleInterfaceVi + + fanIndex * fullPatternInterfaceVertexOffset; + const size_t vi1 = vi0 + (fanSize / 2) * fullPatternInterfaceVertexOffset; + assert(vi0 < fullPattern.VN() && vi1 < fullPattern.VN()); + fullPatternOppositeInterfaceViPairs.emplace_back(std::make_pair(vi0, vi1)); + } + global.fullPatternInterfaceViPairs = fullPatternOppositeInterfaceViPairs; + +#if POLYSCOPE_DEFINED + const bool debugMapping = false; + if (debugMapping) { + reducedPattern.registerForDrawing(); + + // std::vector colors_reducedPatternExcludedEdges( + // reducedPattern.EN(), glm::vec3(0, 0, 0)); + // for (const size_t ei : global.reducedPatternExludedEdges) { + // colors_reducedPatternExcludedEdges[ei] = glm::vec3(1, 0, 0); + // } + // const std::string label = reducedPattern.getLabel(); + // polyscope::getCurveNetwork(label) + // ->addEdgeColorQuantity("Excluded edges", + // colors_reducedPatternExcludedEdges) + // ->setEnabled(true); + // polyscope::show(); + + std::vector nodeColorsOpposite(fullPattern.VN(), glm::vec3(0, 0, 0)); + for (const std::pair oppositeVerts : fullPatternOppositeInterfaceViPairs) { + auto color = polyscope::getNextUniqueColor(); + nodeColorsOpposite[oppositeVerts.first] = color; + nodeColorsOpposite[oppositeVerts.second] = color; + } + fullPattern.registerForDrawing(); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("oppositeMap", nodeColorsOpposite) + ->setEnabled(true); + polyscope::show(); + + std::vector nodeColorsReducedToFull_reduced(reducedPattern.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; + } + } + polyscope::getCurveNetwork(reducedPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_reduced", nodeColorsReducedToFull_reduced) + ->setEnabled(true); + polyscope::getCurveNetwork(fullPattern.getLabel()) + ->addNodeColorQuantity("reducedToFull_full", nodeColorsReducedToFull_full) + ->setEnabled(true); + polyscope::show(); + } +#endif +} + +void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, + PatternGeometry &reducedPattern) +{ + ReducedModelOptimizer::computeMaps(slotToNode, + fullPattern, + reducedPattern, + global.reducedToFullInterfaceViMap, + m_fullToReducedInterfaceViMap, + m_fullPatternOppositeInterfaceViPairs); +} + +ReducedModelOptimizer::ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot) +{ + FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); + FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); + interfaceNodeIndex = numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]; + constructBaseScenarioFunctions.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios); + scenarioIsSymmetrical.resize(BaseSimulationScenario::NumberOfBaseSimulationScenarios); + constructBaseScenarioFunctions[BaseSimulationScenario::Axial] = &constructAxialSimulationScenario; + scenarioIsSymmetrical[BaseSimulationScenario::Axial] = false; + constructBaseScenarioFunctions[BaseSimulationScenario::Shear] = &constructShearSimulationScenario; + scenarioIsSymmetrical[BaseSimulationScenario::Shear] = false; + constructBaseScenarioFunctions[BaseSimulationScenario::Bending] + = &constructBendingSimulationScenario; + scenarioIsSymmetrical[BaseSimulationScenario::Bending] = true; + constructBaseScenarioFunctions[BaseSimulationScenario::Dome] = &constructDomeSimulationScenario; + scenarioIsSymmetrical[BaseSimulationScenario::Dome] = true; + constructBaseScenarioFunctions[BaseSimulationScenario::Saddle] + = &constructSaddleSimulationScenario; + scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true; +} + +void ReducedModelOptimizer::initializePatterns( + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + const std::array &optimizationParameters) +{ + assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN()); + +#if POLYSCOPE_DEFINED + polyscope::removeAllStructures(); +#endif + // Create copies of the input models + PatternGeometry copyFullPattern; + PatternGeometry copyReducedPattern; + copyFullPattern.copy(fullPattern); + copyReducedPattern.copy(reducedPattern); + global.baseTriangle = copyReducedPattern.getBaseTriangle(); + + computeMaps(copyFullPattern, copyReducedPattern); + createSimulationMeshes(copyFullPattern, copyReducedPattern); + initializeUpdateReducedPatternFunctions(); + initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); +} + +void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern, + const ReducedModelOptimization::Settings &optimizationSettings, + ReducedModelOptimization::Results &optimizationResults) +{ + //scale pattern and reduced model + fullPattern.scale(optimizationSettings.targetBaseTriangleSize, interfaceNodeIndex); + ReducedModel reducedModel; + reducedModel.scale(optimizationSettings.targetBaseTriangleSize, interfaceNodeIndex); + auto start = std::chrono::system_clock::now(); + optimizationResults.label = fullPattern.getLabel(); + optimizationResults.baseTriangleFullPattern.copy(fullPattern); + optimizationResults.settings = optimizationSettings; + initializePatterns(fullPattern, reducedModel, optimizationSettings.variablesRanges); + optimize(optimizationSettings, optimizationResults); + auto end = std::chrono::system_clock::now(); + auto elapsed = std::chrono::duration_cast(end - start); + optimizationResults.time = elapsed.count() / 1000.0; +} + +void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() +{ + 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(PatternGeometry::DefaultNormal, + 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(PatternGeometry::DefaultNormal, vcg::math::ToRad(newTheta)) + * baseTriangleHexagonVertexPosition; + + for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; + rotationCounter++) { + pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + = vcg::RotationMatrix(PatternGeometry::DefaultNormal, + 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)); + } + }; + global.functions_updateReducedPatternParameter[A] = + [](const double &newA, std::shared_ptr &pReducedPatternSimulationMesh) { + assert(pReducedPatternSimulationMesh != nullptr); + const double beamWidth = std::sqrt(newA); + const double beamHeight = beamWidth; + for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) { + Element &e = pReducedPatternSimulationMesh->elements[ei]; + e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); + } + }; + + 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.dimensions.inertia.I2 = newI2; + e.updateRigidity(); + } + }; + 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.dimensions.inertia.I3 = newI3; + e.updateRigidity(); + } + }; + 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.dimensions.inertia.J = newJ; + e.updateRigidity(); + } + }; +} + +void ReducedModelOptimizer::initializeOptimizationParameters( + const std::shared_ptr &mesh, + const std::array &optimizationParameters) +{ + for (int optimizationParameterIndex = 0; + optimizationParameterIndex < optimizationParameters.size(); + optimizationParameterIndex++) { + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationVariables; + optimizationParameterIndex++) { + // const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex]; + switch (optimizationParameterIndex) { + case E: + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].material.youngsModulus; + global.optimizationInitialValue[optimizationParameterIndex] = 1; + break; + case A: + global.parametersInitialValue[optimizationParameterIndex] = mesh->elements[0] + .dimensions.A; + global.optimizationInitialValue[optimizationParameterIndex] = 1; + break; + case I2: + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.inertia.I2; + global.optimizationInitialValue[optimizationParameterIndex] = 1; + break; + case I3: + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.inertia.I3; + global.optimizationInitialValue[optimizationParameterIndex] = 1; + break; + case J: + global.parametersInitialValue[optimizationParameterIndex] + = mesh->elements[0].dimensions.inertia.J; + global.optimizationInitialValue[optimizationParameterIndex] = 1; + break; + case R: + global.parametersInitialValue[optimizationParameterIndex] = 0; + global.optimizationInitialValue[optimizationParameterIndex] = 0.5; + break; + case Theta: + global.parametersInitialValue[optimizationParameterIndex] = 0; + global.optimizationInitialValue[optimizationParameterIndex] = 0; + break; + } + } + } +} + +void ReducedModelOptimizer::computeReducedModelSimulationJob( + const SimulationJob &simulationJobOfFullModel, + const std::unordered_map &fullToReducedMap, + SimulationJob &simulationJobOfReducedModel) +{ + assert(simulationJobOfReducedModel.pMesh->VN() != 0); + 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) { + reducedModelNodalForces[fullToReducedMap.at(fullModelNodalForce.first)] + = fullModelNodalForce.second; + } + + std::unordered_map reducedNodalForcedDisplacements; + for (auto fullForcedDisplacement : simulationJobOfFullModel.nodalForcedDisplacements) { + reducedNodalForcedDisplacements[fullToReducedMap.at(fullForcedDisplacement.first)] + = fullForcedDisplacement.second; + } + + simulationJobOfReducedModel.constrainedVertices = reducedModelFixedVertices; + simulationJobOfReducedModel.nodalExternalForces = reducedModelNodalForces; + simulationJobOfReducedModel.label = simulationJobOfFullModel.getLabel(); + simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements; +} + +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()); + + 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]); + } +} + +void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjective, + const Settings &settings, + ReducedModelOptimization::Results &results) +{ + std::shared_ptr &pReducedPatternSimulationMesh + = global.reducedPatternSimulationJobs[global.simulationScenarioIndices[0]]->pMesh; + //Number of crashes + // results.numberOfSimulationCrashes = global.numOfSimulationCrashes; + //Value of optimal objective Y + results.objectiveValue.total = optimalObjective.y; + // results.objectiveValue.total = 0; + if (optimalObjective.y != global.minY) { + std::cout << "Different optimal y:" << optimalObjective.y << " " << global.minY + << std::endl; + } + //Optimal X values + results.optimalXNameValuePairs.resize(NumberOfOptimizationVariables); + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationVariables; + optimizationParameterIndex++) { + if (optimizationParameterIndex != OptimizationParameterIndex::R + && optimizationParameterIndex != OptimizationParameterIndex::Theta) { + results.optimalXNameValuePairs[optimizationParameterIndex] + = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label, + optimalObjective.x[optimizationParameterIndex] + * global.parametersInitialValue[optimizationParameterIndex]); + } else { + //Hex size and angle are pure values (not multipliers of the initial values) + results.optimalXNameValuePairs[optimizationParameterIndex] + = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label, + optimalObjective.x[optimizationParameterIndex]); + } + + global.functions_updateReducedPatternParameter[optimizationParameterIndex]( + results.optimalXNameValuePairs[optimizationParameterIndex].second, + pReducedPatternSimulationMesh); +#ifdef POLYSCOPE_DEFINED + std::cout << results.optimalXNameValuePairs[optimizationParameterIndex].first << ":" + << results.optimalXNameValuePairs[optimizationParameterIndex].second << " "; +#endif + } +#ifdef POLYSCOPE_DEFINED + std::cout << std::endl; +#endif + pReducedPatternSimulationMesh->reset(); + + results.fullPatternYoungsModulus = youngsModulus; + // Compute obj value per simulation scenario and the raw objective value + // updateMeshFunction(optimalX); + // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); + //TODO:use push_back it will make the code more readable + LinearSimulationModel simulator; + results.objectiveValue.totalRaw = 0; + results.objectiveValue.perSimulationScenario_translational.resize( + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_rawTranslational.resize( + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_rotational.resize( + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_rawRotational.resize( + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_total.resize( + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_total_unweighted.resize( + global.simulationScenarioIndices.size()); + //#ifdef POLYSCOPE_DEFINED + // global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed); + //#endif + + results.perScenario_fullPatternPotentialEnergy.resize(global.simulationScenarioIndices.size()); + for (int i = 0; i < global.simulationScenarioIndices.size(); i++) { + const int simulationScenarioIndex = global.simulationScenarioIndices[i]; + SimulationResults reducedModelResults = simulator.executeSimulation( + global.reducedPatternSimulationJobs[simulationScenarioIndex]); + +#ifdef POLYSCOPE_DEFINED +#ifdef USE_SCENARIO_WEIGHTS + const double scenarioWeight = global.scenarioWeights[simulationScenarioIndex]; +#else + const double scenarioWeight = 1; +#endif +#else + const double scenarioWeight = 1; +#endif + results.objectiveValue.perSimulationScenario_total[i] = computeError( + global.fullPatternResults[simulationScenarioIndex], + reducedModelResults, + global.reducedToFullInterfaceViMap, + global.translationalDisplacementNormalizationValues[simulationScenarioIndex], + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + scenarioWeight, + global.objectiveWeights[simulationScenarioIndex]); + + results.objectiveValue.perSimulationScenario_total_unweighted[i] = computeError( + global.fullPatternResults[simulationScenarioIndex], + reducedModelResults, + global.reducedToFullInterfaceViMap, + global.translationalDisplacementNormalizationValues[simulationScenarioIndex], + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex], + 1, + global.objectiveWeights[simulationScenarioIndex]); + + // results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i]; + //Raw translational + const double rawTranslationalError = computeRawTranslationalError( + global.fullPatternResults[simulationScenarioIndex].displacements, + reducedModelResults.displacements, + global.reducedToFullInterfaceViMap); + results.objectiveValue.perSimulationScenario_rawTranslational[i] = rawTranslationalError; + //Raw rotational + const double rawRotationalError = computeRawRotationalError( + global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, + reducedModelResults.rotationalDisplacementQuaternion, + global.reducedToFullInterfaceViMap); + results.objectiveValue.perSimulationScenario_rawRotational[i] = rawRotationalError; + + //Normalized translational + const double normalizedTranslationalError = computeDisplacementError( + global.fullPatternResults[simulationScenarioIndex].displacements, + reducedModelResults.displacements, + global.reducedToFullInterfaceViMap, + global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); + results.objectiveValue.perSimulationScenario_translational[i] = normalizedTranslationalError; + // const double test_normalizedTranslationError = computeDisplacementError( + // global.fullPatternResults[simulationScenarioIndex].displacements, + // reducedModelResults.displacements, + // global.reducedToFullInterfaceViMap, + // global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); + //Normalized rotational + const double normalizedRotationalError = computeRotationalError( + global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, + reducedModelResults.rotationalDisplacementQuaternion, + global.reducedToFullInterfaceViMap, + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + results.objectiveValue.perSimulationScenario_rotational[i] = normalizedRotationalError; + // const double test_normalizedRotationalError = computeRotationalError( + // global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, + // reducedModelResults.rotationalDisplacementQuaternion, + // global.reducedToFullInterfaceViMap, + // global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); + // assert(test_normalizedTranslationError == normalizedTranslationalError); + // assert(test_normalizedRotationalError == normalizedRotationalError); + results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError; + results.perScenario_fullPatternPotentialEnergy[i] + = global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy; +#ifdef POLYSCOPE_DEFINED + std::cout << "Simulation scenario:" + << global.reducedPatternSimulationJobs[simulationScenarioIndex]->getLabel() + << std::endl; + std::cout << "raw translational error:" << rawTranslationalError << std::endl; + std::cout << "translation normalization value:" + << global.translationalDisplacementNormalizationValues[simulationScenarioIndex] + << std::endl; + std::cout << "raw Rotational error:" << rawRotationalError << std::endl; + std::cout << "rotational normalization value:" + << global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] + << std::endl; + std::cout << "Translational error:" << normalizedTranslationalError << std::endl; + std::cout << "Rotational error:" << normalizedRotationalError << std::endl; + // results.objectiveValuePerSimulationScenario[simulationScenarioIndex] + // = normalizedTranslationalError + normalizedRotationalError; + std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i] + << std::endl; + std::cout << std::endl; + + // reducedModelResults.registerForDrawing(Colors::reducedDeformed); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed); + // polyscope::show(); + // reducedModelResults.unregister(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); + +#endif + } + + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + results.fullPatternSimulationJobs.push_back( + global.fullPatternSimulationJobs[simulationScenarioIndex]); + results.reducedPatternSimulationJobs.push_back( + global.reducedPatternSimulationJobs[simulationScenarioIndex]); + // const std::string temp = global.reducedPatternSimulationJobs[simulationScenarioIndex] + // ->pMesh->getLabel(); + // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel("temp"); + // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->registerForDrawing(); + // global.reducedPatternSimulationJobs[simulationScenarioIndex]->pMesh->setLabel(temp); + } + results.objectiveValueHistory = global.objectiveValueHistory; + results.objectiveValueHistory_iteration = global.objectiveValueHistory_iteration; + // results.draw(); +} + +std::array +ReducedModelOptimizer::computeFullPatternMaxSimulationForces( + const std::vector &desiredBaseSimulationScenario) const +{ + std::array fullPatternMaxSimulationForces; + for (const BaseSimulationScenario &scenario : desiredBaseSimulationScenario) { + const double maxForce = computeFullPatternMaxSimulationForce(scenario); + fullPatternMaxSimulationForces[scenario] = maxForce; + } + + return fullPatternMaxSimulationForces; +} + +std::array +ReducedModelOptimizer::getFullPatternMaxSimulationForces( + const std::vector &desiredBaseSimulationScenarioIndices, + const std::filesystem::path &intermediateResultsDirectoryPath, + const bool &recomputeForceMagnitudes) +{ + std::array fullPatternSimulationScenarioMaxMagnitudes; + //#ifdef POLYSCOPE_DEFINED + const std::filesystem::path forceMagnitudesDirectoryPath( + std::filesystem::path(intermediateResultsDirectoryPath).append("ForceMagnitudes")); + std::filesystem::path patternMaxForceMagnitudesFilePath( + std::filesystem::path(forceMagnitudesDirectoryPath) + .append(m_pFullPatternSimulationMesh->getLabel() + ".json")); + const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( + patternMaxForceMagnitudesFilePath); + if (fullPatternScenarioMagnitudesExist && !recomputeForceMagnitudes) { + nlohmann::json json; + std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); + ifs >> json; + fullPatternSimulationScenarioMaxMagnitudes + = static_cast>(json.at("maxMagn")); + return fullPatternSimulationScenarioMaxMagnitudes; + } + + //#endif + fullPatternSimulationScenarioMaxMagnitudes = computeFullPatternMaxSimulationForces( + desiredBaseSimulationScenarioIndices); + + //#ifdef POLYSCOPE_DEFINED + nlohmann::json json; + json["maxMagn"] = fullPatternSimulationScenarioMaxMagnitudes; + + std::filesystem::create_directories(forceMagnitudesDirectoryPath); + std::ofstream jsonFile(patternMaxForceMagnitudesFilePath.string()); + jsonFile << json; + +#ifdef POLYSCOPE_DEFINED + std::cout << "Saved base scenario max magnitudes to:" << patternMaxForceMagnitudesFilePath + << std::endl; +#endif + + //#endif + assert(fullPatternSimulationScenarioMaxMagnitudes.size() + == desiredBaseSimulationScenarioIndices.size()); + + return fullPatternSimulationScenarioMaxMagnitudes; +} + +void ReducedModelOptimizer::runOptimization(const Settings &settings, + ReducedModelOptimization::Results &results) +{ + global.objectiveValueHistory.clear(); + global.objectiveValueHistory_iteration.clear(); + global.objectiveValueHistory.reserve(settings.numberOfFunctionCalls / 100); + global.objectiveValueHistory_iteration.reserve(settings.numberOfFunctionCalls / 100); + +#if POLYSCOPE_DEFINED +// global.plotColors.reserve(settings.numberOfFunctionCalls); +#endif + assert(!settings.optimizationStrategy.empty()); + const std::vector> &optimizationParametersGroups + = settings.optimizationStrategy; + +#ifndef USE_ENSMALLEN + const int totalNumberOfOptimizationParameters + = std::accumulate(optimizationParametersGroups.begin(), + optimizationParametersGroups.end(), + 0, + [](const int &sum, + const std::vector ¶meterGroup) { + return sum + parameterGroup.size(); + }); +#endif + + FunctionEvaluation optimization_optimalResult; + optimization_optimalResult.x.resize(NumberOfOptimizationVariables, 0); + for (int optimizationParameterIndex = E; + optimizationParameterIndex != NumberOfOptimizationVariables; + optimizationParameterIndex++) { + optimization_optimalResult.x[optimizationParameterIndex] + = global.optimizationInitialValue[optimizationParameterIndex]; + } + for (int parameterGroupIndex=0;parameterGroupIndex ¶meterGroup = + optimizationParametersGroups[parameterGroupIndex]; + FunctionEvaluation parameterGroup_optimalResult; + //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup + function_updateReducedPattern = [&](const std::vector &x, + std::shared_ptr &pMesh) { + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + // const double parameterInitialValue=optimizationSettings.parameterRanges[parameterIndex].initialValue; + const double parameterNewValue = [&]() { + if (parameterIndex == R || parameterIndex == Theta) { + return x[xIndex] /*+ parameterInitialValue*/; + } + //and the material parameters exponentially(?).TODO: Check what happens if I make all linear + const double parameterInitialValue + = global.parametersInitialValue[parameterIndex]; + return x[xIndex] * parameterInitialValue; + }(); + // std::cout << "Optimization parameter:" << parameterIndex << std::endl; + // std::cout << "New value:" << parameterNewValue << std::endl; + global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, + pMesh); + } + pMesh->reset(); //NOTE: I could put this code into each updateParameter function for avoiding unessecary calculations + }; + + std::vector xMin; + std::vector xMax; + xMin.resize(parameterGroup.size()); + xMax.resize(parameterGroup.size()); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + + xMin[xIndex] = settings.variablesRanges[parameterIndex].min; + xMax[xIndex] = settings.variablesRanges[parameterIndex].max; + } +#ifdef USE_ENSMALLEN + arma::mat x(parameterGroup.size(), 1); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + + x(xIndex, 0) = global.optimizationInitialValue[parameterIndex]; + } + + global.xMin = xMin; + global.xMax = xMax; + + // Create simulated annealing optimizer with default options. + // The ens::SA<> type can be replaced with any suitable ensmallen optimizer + // that is able to handle arbitrary functions. + EnsmallenOptimizationObjective optimizationFunction; + //Set min max values + // ens::SA optimizer; + // ens::CNE optimizer; + // ens::DE optimizer; + // ens::SPSA optimizer; + arma::mat xMin_arma(global.xMin); + arma::mat xMax_arma(global.xMax); + // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); + ens::LBestPSO optimizer(100, + xMin_arma, + xMax_arma, + settings.numberOfFunctionCalls, + 500, + settings.solverAccuracy, + 2.05, + 2.35); + // ens::LBestPSO optimizer; + const double minima = optimizer.Optimize(optimizationFunction, x); + // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + // if (x[xIndex] > global.xMax[xIndex]) { + // x[xIndex] = global.xMax[xIndex]; + // } else if (x[xIndex] < global.xMin[xIndex]) { + // x[xIndex] = global.xMin[xIndex]; + // } + // } + +#ifdef POLYSCOPE_DEFINED + std::cout << "optimal x:" + << "\n" + << x << std::endl; + std::cout << "Minima:" << minima << std::endl; + std::cout << "optimal objective:" << optimizationFunction.Evaluate(x) << std::endl; +#endif + parameterGroup_optimalResult.x.clear(); + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(x.begin(), x.end(), parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = minima; +#else + //Set min max values + dlib::matrix xMin_dlib(parameterGroup.size()); + dlib::matrix xMax_dlib(parameterGroup.size()); + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + + xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min; + xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max; + } + const double portionOfTotalFunctionCallsBudget = [&]() { + if (settings.optimizationVariablesGroupsWeights.has_value()) { + assert(settings.optimizationVariablesGroupsWeights->size() + == settings.optimizationStrategy.size()); + return settings.optimizationVariablesGroupsWeights.value()[parameterGroupIndex]; + + } else { + return static_cast(parameterGroup.size()) + / totalNumberOfOptimizationParameters; + } + }(); + const int optimizationNumberOfFunctionCalls = portionOfTotalFunctionCallsBudget + * settings.numberOfFunctionCalls; + const dlib::function_evaluation optimalResult_dlib = [&]() { + if (parameterGroup.size() == 1) { + dlib::function_evaluation result; + double optimalX = global.optimizationInitialValue[parameterGroup[0]]; + double (*singleVariableObjective)(const double &) = &objective; + try { + result.y = dlib::find_min_single_variable(singleVariableObjective, + optimalX, + xMin_dlib(0), + xMax_dlib(0), + 1e-5, + optimizationNumberOfFunctionCalls); + } catch (const std::exception &e) { // reference to the base of a polymorphic object +#ifdef POLYSCOPE_DEFINED + std::cout << e.what() << std::endl; +#endif + } + + result.x.set_size(1); + result.x(0) = optimalX; + return result; + } + double (*objF)(const dlib::matrix &) = &objective; + return dlib::find_min_global(objF, + xMin_dlib, + xMax_dlib, + dlib::max_function_calls(optimizationNumberOfFunctionCalls), + 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 = 10000; + // dlib::matrix x; + // dlib::function_evaluation optimalResult_dlib; + // optimalResult_dlib.x.set_size(parameterGroup.size()); + // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + // const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + // optimalResult_dlib.x(xIndex) = global.optimizationInitialValue[parameterIndex]; + + // std::cout << "xIndex:" << xIndex << std::endl; + // std::cout << "xInit:" << optimalResult_dlib.x(xIndex) << std::endl; + // } + + // optimalResult_dlib.y = dlib::find_min_bobyqa(objF, + // optimalResult_dlib.x, + // npt, + // xMin, + // xMax, + // rhobeg, + // rhoend, + // bobyqa_maxFunctionCalls); + // } + + parameterGroup_optimalResult.x.clear(); + parameterGroup_optimalResult.x.resize(parameterGroup.size()); + std::copy(optimalResult_dlib.x.begin(), + optimalResult_dlib.x.end(), + parameterGroup_optimalResult.x.begin()); + parameterGroup_optimalResult.y = optimalResult_dlib.y; + +#endif + + const auto firstNonNullReducedSimulationJobIt + = std::find_if(global.reducedPatternSimulationJobs.begin(), + global.reducedPatternSimulationJobs.end(), + [](const std::shared_ptr &pJob) { + return pJob != nullptr; + }); + function_updateReducedPattern( + std::vector(parameterGroup_optimalResult.x.begin(), + parameterGroup_optimalResult.x.end()), + (*firstNonNullReducedSimulationJobIt) + ->pMesh); //TODO: Check if its ok to update only the changed parameters +#ifdef POLYSCOPE_DEFINED + std::cout << "Optimal x:" + << "\n"; + for (const auto &x : parameterGroup_optimalResult.x) { + std::cout << x << std::endl; + } + std::cout << "optimal objective:" << objective(parameterGroup_optimalResult.x) << std::endl; + // std::cout << "Minima:" << minima << std::endl; + std::cout << "GLOBAL MIN:" << global.minY << std::endl; + std::cout << "opt res y:" << parameterGroup_optimalResult.y << std::endl; +#endif + +optimization_optimalResult.y=parameterGroup_optimalResult.y; + for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex]; + optimization_optimalResult.x[parameterIndex] = parameterGroup_optimalResult.x[xIndex]; + } + } + getResults(optimization_optimalResult, settings, results); +} + +void ReducedModelOptimizer::constructAxialSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job) +{ + for (auto viPairIt = oppositeInterfaceViPairs.begin(); + viPairIt != oppositeInterfaceViPairs.end(); + viPairIt++) { + if (viPairIt != oppositeInterfaceViPairs.begin()) { + CoordType forceDirection(1, 0, 0); + job.nodalExternalForces[viPairIt->first] + = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) + * forceMagnitude; + job.constrainedVertices[viPairIt->second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } + } +} + +void ReducedModelOptimizer::constructShearSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job) +{ + for (auto viPairIt = oppositeInterfaceViPairs.begin(); + viPairIt != oppositeInterfaceViPairs.end(); + viPairIt++) { + if (viPairIt != oppositeInterfaceViPairs.begin()) { + CoordType forceDirection(0, 1, 0); + const auto viPair = *viPairIt; + job.nodalExternalForces[viPair.first] + = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) + * forceMagnitude; + job.constrainedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } + } +} + +void ReducedModelOptimizer::constructBendingSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job) +{ + for (const auto &viPair : oppositeInterfaceViPairs) { + job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude; + job.constrainedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } +} + +/*NOTE: From the results it seems as if the forced displacements are different in the linear and in the drm + * */ +void ReducedModelOptimizer::constructDomeSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job) +{ + for (auto viPairIt = oppositeInterfaceViPairs.begin(); + viPairIt != oppositeInterfaceViPairs.end(); + viPairIt++) { + const auto viPair = *viPairIt; + CoordType interfaceVector = (job.pMesh->vert[viPair.first].cP() + - job.pMesh->vert[viPair.second].cP()); + VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0)) + * interfaceVector.Normalize(); + if (viPairIt == oppositeInterfaceViPairs.begin()) { + job.nodalForcedDisplacements[viPair.first] + = Eigen::Vector3d(-interfaceVector[0], + -interfaceVector[1], + 0) + * 0.001 + * std::abs( + forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude? + // * std::abs(forceMagnitude / maxForceMagnitude_dome); + job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0], + interfaceVector[1], + 0) + * 0.001 * std::abs(forceMagnitude); + // * std::abs(forceMagnitude / maxForceMagnitude_dome); + // CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + // ^ CoordType(0, 0, -1).Normalize(); + // nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude + // * 0.0001; + // nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude + // * 0.0001; + } else { + job.nodalExternalForces[viPair.first] + = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude + / 3; + job.nodalExternalForces[viPair.second] + = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) + * forceMagnitude / 3; + job.constrainedVertices[viPair.first] = std::unordered_set{2}; + job.constrainedVertices[viPair.second] = std::unordered_set{2}; + } + } +} + +void ReducedModelOptimizer::constructSaddleSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job) +{ + for (auto viPairIt = oppositeInterfaceViPairs.begin(); + viPairIt != oppositeInterfaceViPairs.end(); + viPairIt++) { + const auto &viPair = *viPairIt; + CoordType v = (job.pMesh->vert[viPair.first].cP() - job.pMesh->vert[viPair.second].cP()) + ^ CoordType(0, 0, -1).Normalize(); + if (viPairIt == oppositeInterfaceViPairs.begin()) { + job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) + * forceMagnitude; + job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) + * forceMagnitude; + } else { + job.constrainedVertices[viPair.first] = std::unordered_set{2}; + job.constrainedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + + job.nodalExternalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) + * forceMagnitude / 2; + job.nodalExternalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) + * forceMagnitude / 2; + } + } +} + +double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagnitude) +{ + SimulationJob job; + job.pMesh = global.pFullPatternSimulationMesh; + global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job); + // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, + // global.fullPatternInterfaceViPairs, + // job); + + DRMSimulationModel simulator; + DRMSimulationModel::Settings settings; + // settings.totalExternalForcesNormPercentageTermination = 1e-2; + settings.totalTranslationalKineticEnergyThreshold = 1e-10; + // settings.viscousDampingFactor = 1e-2; + // settings.useKineticDamping = true; + settings.shouldDraw = false; + settings.debugModeStep = 200; + // settings.averageResidualForcesCriterionThreshold = 1e-5; + settings.maxDRMIterations = 200000; + SimulationResults results = simulator.executeSimulation(std::make_shared(job), + settings); + const double &desiredRotationAngle = global.desiredMaxRotationAngle; + double error; + // job.pMesh->setLabel("initial"); + // job.pMesh->registerForDrawing(); + // results.registerForDrawing(); + // polyscope::show(); + std::string saveJobToPath; + if (!results.converged) { + // std::cout << "Force used:" << forceMagnitude << std::endl; + error = std::numeric_limits::max(); + // DRMSimulationModel::Settings debugSimulationSettings; + // debugSimulationSettings.isDebugMode = true; + // debugSimulationSettings.debugModeStep = 1000; + // // debugSimulationSettings.maxDRMIterations = 100000; + // debugSimulationSettings.shouldDraw = true; + // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; + // debugSimulationSettings.shouldCreatePlots = true; + // // debugSimulationSettings.Dtini = 0.06; + // debugSimulationSettings.beVerbose = true; + // debugSimulationSettings.useAverage = true; + // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; + // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; + // auto debugResults = simulator.executeSimulation(std::make_shared(job), + // debugSimulationSettings); + // std::terminate(); + saveJobToPath = "../nonConvergingJobs"; + } else { + error = std::abs( + results.rotationalDisplacementQuaternion[global.interfaceViForComputingScenarioError] + .angularDistance(Eigen::Quaterniond::Identity()) + - desiredRotationAngle); + saveJobToPath = "../convergingJobs"; + } + // std::filesystem::path outputPath(std::filesystem::path(saveJobToPath) + // .append(job.pMesh->getLabel()) + // .append("mag_" + global.currentScenarioName)); + // std::filesystem::create_directories(outputPath); + // job.save(outputPath); + // settings.save(outputPath); + +#ifdef POLYSCOPE_DEFINED + std::cout << "Force:" << forceMagnitude << " Error is:" << vcg::math::ToDeg(error) << std::endl; +#endif + return error; +} + +void search(const std::function objectiveFunction, + const double &targetY, + double &optimalX, + double &error, + const double &epsilon, + const double &maxIterations) +{ + error = std::numeric_limits::max(); + int iterationIndex = 0; + double low = optimalX / 1e-3, high = optimalX * 1e3; + while (error > epsilon && iterationIndex < maxIterations) { + const double y = objectiveFunction(optimalX); + error = std::abs(targetY - y); + if (y > targetY) { + high = optimalX; + } else { + low = optimalX; + } + optimalX = low + (high - low) / 2; + + iterationIndex++; + } + + if (iterationIndex == maxIterations) { + std::cerr << "Max iterations reached." << std::endl; + } +} + +double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMagnitude) +{ + SimulationJob job; + job.pMesh = global.pFullPatternSimulationMesh; + global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job); + // ReducedModelOptimizer::axialSimulationScenario(forceMagnitude, + // global.fullPatternInterfaceViPairs, + // job); + + DRMSimulationModel simulator; + DRMSimulationModel::Settings settings; + // settings.totalResidualForcesNormThreshold = 1e-3; + settings.totalExternalForcesNormPercentageTermination = 1e-2; + settings.totalTranslationalKineticEnergyThreshold = 1e-8; + settings.viscousDampingFactor = 5e-3; + settings.useKineticDamping = true; + // settings.averageResidualForcesCriterionThreshold = 1e-5; + // settings.useAverage = true; + // settings.totalTranslationalKineticEnergyThreshold = 1e-10; + // settings.shouldUseTranslationalKineticEnergyThreshold = true; + // settings.shouldDraw = true; + // settings.isDebugMode = true; + // settings.drawingStep = 200000; + // settings.beVerbose = true; + // settings.debugModeStep = 200; + settings.maxDRMIterations = 200000; + SimulationResults results = simulator.executeSimulation(std::make_shared(job), + settings); + const double &desiredDisplacementValue = global.desiredMaxDisplacementValue; + double error; + if (!results.converged) { + error = std::numeric_limits::max(); + std::filesystem::path outputPath(std::filesystem::path("../nonConvergingJobs") + .append(job.pMesh->getLabel()) + .append("mag_" + global.currentScenarioName)); + std::filesystem::create_directories(outputPath); + job.save(outputPath); + settings.save(outputPath); + // std::terminate(); + } else { + error = std::abs( + results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm() + - desiredDisplacementValue); + } + +#ifdef POLYSCOPE_DEFINED + std::cout << "Force:" << forceMagnitude << " Error is:" << error << std::endl; +#endif + + return error; +} + +#ifdef USE_ENSMALLEN +struct ForceMagnitudeOptimization +{ + std::function objectiveFunction; + ForceMagnitudeOptimization(std::function &f) : objectiveFunction(f) {} + double Evaluate(const arma::mat &x) { return objectiveFunction(x(0, 0)); } +}; +#endif + +double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( + const BaseSimulationScenario &scenario) const +{ + double forceMagnitude = 1; + double minimumError; + bool wasSuccessful = false; + global.constructScenarioFunction = constructBaseScenarioFunctions[scenario]; + const double baseTriangleHeight = vcg::Distance(global.baseTriangle.cP(0), + (global.baseTriangle.cP(1) + + global.baseTriangle.cP(2)) + / 2); + std::function objectiveFunction; + double translationalOptimizationEpsilon{baseTriangleHeight * 0.001}; + double objectiveEpsilon = translationalOptimizationEpsilon; + objectiveFunction = &fullPatternMaxSimulationForceTranslationalObjective; + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; + global.desiredMaxDisplacementValue = 0.1 * baseTriangleHeight; + global.currentScenarioName = baseSimulationScenarioNames[scenario]; + double forceMagnitudeEpsilon = 1e-4; + + switch (scenario) { + case Axial: + global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight; + break; + case Shear: + global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight; + break; + case Bending: + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; + global.desiredMaxDisplacementValue = 0.05 * baseTriangleHeight; + break; + case Dome: + global.desiredMaxRotationAngle = vcg::math::ToRad(20.0); + objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective; + forceMagnitudeEpsilon *= 1e-2; + objectiveEpsilon = vcg::math::ToRad(1.0); + forceMagnitude = 0.6; + break; + case Saddle: + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; + global.desiredMaxDisplacementValue = 0.05 * baseTriangleHeight; + break; + } + + constexpr int maxIterations = 1000; + minimumError = dlib::find_min_single_variable(objectiveFunction, + forceMagnitude, + 1e-4, + 1e4, + forceMagnitudeEpsilon, + maxIterations); +#ifdef DLIB_DEFINED +#else + // ens::SA<> optimizer; + // arma::vec lowerBound("0.00001"); + // arma::vec upperBound("10000"); + // ens::LBestPSO optimizer(64, lowerBound, upperBound, maxIterations, 350, forceMagnitudeEpsilon); + // ForceMagnitudeOptimization f(objectiveFunction); // Create function to be optimized. + // arma::mat forceMagnitude_mat({forceMagnitude}); + // minimumError = optimizer.Optimize(f, forceMagnitude_mat); + // std::cout << ReducedModelOptimization::baseSimulationScenarioNames[scenario] << ": " + // << optimalObjective << std::endl; + // forceMagnitude = forceMagnitude_mat(0, 0); + +// search(objectiveFunction, +// global.desiredMaxDisplacementValue, +// forceMagnitude, +// minimumError, +// objectiveEpsilon, +// maxIterations); +#endif + + wasSuccessful = minimumError < objectiveEpsilon; + +#ifdef POLYSCOPE_DEFINED + std::cout << "Max " << ReducedModelOptimization::baseSimulationScenarioNames[scenario] + << " magnitude:" << forceMagnitude << std::endl; + if (!wasSuccessful) { + SimulationJob job; + job.pMesh = global.pFullPatternSimulationMesh; + global.constructScenarioFunction(forceMagnitude, global.fullPatternInterfaceViPairs, job); + std::cerr << ReducedModelOptimization::baseSimulationScenarioNames[scenario] + + " max scenario magnitude was not succefully determined." + << std::endl; + std::filesystem::path outputPath( + std::filesystem::path("../nonConvergingJobs") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append("magFinal_" + + ReducedModelOptimization::baseSimulationScenarioNames[scenario])); + std::filesystem::create_directories(outputPath); + job.save(outputPath); + std::terminate(); + } +#endif + + return forceMagnitude; +} + +void ReducedModelOptimizer::computeScenarioWeights( + const std::vector &baseSimulationScenarios) +{ + std::array baseScenario_equalizationWeights; + baseScenario_equalizationWeights.fill(0); + std::array baseScenario_userWeights; + baseScenario_userWeights.fill(0); + //Fill only used base scenario weights + for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { + baseScenario_equalizationWeights[baseScenario] = static_cast( + simulationScenariosResolution[baseScenario]); + baseScenario_userWeights[baseScenario] = baseScenarioWeights[baseScenario]; + } + + Utilities::normalize(baseScenario_userWeights.begin(), baseScenario_userWeights.end()); + + Utilities::normalize(baseScenario_equalizationWeights.begin(), + baseScenario_equalizationWeights.end()); + std::transform(baseScenario_equalizationWeights.begin(), + baseScenario_equalizationWeights.end(), + baseScenario_equalizationWeights.begin(), + [](const double &d) { + if (d == 0) { + return d; + } + return 1 / d; + }); + Utilities::normalize(baseScenario_equalizationWeights.begin(), + baseScenario_equalizationWeights.end()); + + std::array baseScenarios_weights; + baseScenarios_weights.fill(0); + for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { + baseScenarios_weights[baseScenario] = baseScenario_equalizationWeights[baseScenario] + + 2 * baseScenario_userWeights[baseScenario]; + } + Utilities::normalize(baseScenarios_weights.begin(), baseScenarios_weights.end()); + + global.objectiveWeights.resize(totalNumberOfSimulationScenarios); + global.scenarioWeights.resize( + totalNumberOfSimulationScenarios, + 0); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization + for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) { + const int baseSimulationScenarioIndexOffset + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + baseScenario, + 0); + for (int simulationScenarioIndex = 0; + simulationScenarioIndex < simulationScenariosResolution[baseScenario]; + simulationScenarioIndex++) { + const int globalScenarioIndex = baseSimulationScenarioIndexOffset + + simulationScenarioIndex; + global.scenarioWeights[globalScenarioIndex] = baseScenarios_weights[baseScenario]; + + global.objectiveWeights[globalScenarioIndex] + = global.optimizationSettings.perBaseScenarioObjectiveWeights[baseScenario]; + } + } + + //Dont normalize since we want the base scenarios to be normalized not the "global scenarios" + // Utilities::normalize(global.scenarioWeights.begin(), global.scenarioWeights.end()); +} + +std::vector> ReducedModelOptimizer::createFullPatternSimulationJobs( + const std::shared_ptr &pMesh, + const std::array &baseScenarioMaxForceMagnitudes) const +{ + std::vector> scenarios; + scenarios.resize(totalNumberOfSimulationScenarios); + + SimulationJob job; + job.pMesh = pMesh; + + for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) { + const double maxForceMagnitude = baseScenarioMaxForceMagnitudes[baseScenario]; + const int numberOfSimulationScenarios = simulationScenariosResolution[baseScenario]; + const double minForceMagnitude = scenarioIsSymmetrical[baseScenario] + ? maxForceMagnitude / numberOfSimulationScenarios + : -maxForceMagnitude; + const double forceMagnitudeStep = numberOfSimulationScenarios == 1 + ? maxForceMagnitude + : (maxForceMagnitude - minForceMagnitude) + / (numberOfSimulationScenarios); + const int baseSimulationScenarioIndexOffset + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + baseScenario, + 0); + for (int simulationScenarioIndex = 0; simulationScenarioIndex < numberOfSimulationScenarios; + simulationScenarioIndex++) { + const int scenarioIndex = baseSimulationScenarioIndexOffset + simulationScenarioIndex; + if (baseScenarioMaxForceMagnitudes[baseScenario] == -1) { + scenarios[scenarioIndex] = nullptr; + continue; + } + job.nodalExternalForces.clear(); + job.constrainedVertices.clear(); + job.nodalForcedDisplacements.clear(); + job.label = baseSimulationScenarioNames[baseScenario] + "_" + + std::to_string(simulationScenarioIndex); + + const double forceMagnitude = (forceMagnitudeStep * simulationScenarioIndex + + minForceMagnitude); + constructBaseScenarioFunctions[baseScenario](forceMagnitude, + m_fullPatternOppositeInterfaceViPairs, + job); + + scenarios[scenarioIndex] = std::make_shared(job); + } + } + +#ifdef POLYSCOPE_DEFINED + std::cout << "Computed full pattern scenario magnitudes" << std::endl; +#endif + return scenarios; +} + +void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() +{ + // m_pFullPatternSimulationMesh->registerForDrawing(); + // m_pFullPatternSimulationMesh->save(std::filesystem::current_path().append("initial.ply")); + // 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]; + // std::cout << "vi:" << fullPatternVi << std::endl; + // std::cout << "displacement:" << vertexDisplacement.getTranslation().norm() << std::endl; + + translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); + } + // global.fullPatternResults[simulationScenarioIndex].saveDeformedModel( + // std::filesystem::current_path()); + // global.fullPatternResults[simulationScenarioIndex].registerForDrawing(); + // polyscope::show(); + // global.fullPatternResults[simulationScenarioIndex].unregister(); + 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 += std::abs(global.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion[fullPatternVi] + .angularDistance(Eigen::Quaterniond::Identity())); + } + + fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] + = translationalDisplacementNormSum; + fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; + } + + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + if (global.optimizationSettings.normalizationStrategy + == Settings::NormalizationStrategy::Epsilon) { + const double epsilon_translationalDisplacement = global.optimizationSettings + .translationEpsilon; + global.translationalDisplacementNormalizationValues[simulationScenarioIndex] + = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], + epsilon_translationalDisplacement); + const double epsilon_rotationalDisplacement = global.optimizationSettings + .angularDistanceEpsilon; + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] + = std::max(fullPatternAngularDistance[simulationScenarioIndex], + epsilon_rotationalDisplacement); + } else { + global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; + } + } +} + +void ReducedModelOptimizer::optimize( + const Settings &optimizationSettings, + ReducedModelOptimization::Results &results, + const std::vector &desiredBaseSimulationScenarioIndices) +{ + assert(!optimizationSettings.optimizationStrategy.empty()); + assert(!optimizationSettings.optimizationVariablesGroupsWeights.has_value() + || (optimizationSettings.optimizationStrategy.size() + == optimizationSettings.optimizationVariablesGroupsWeights->size() + && std::accumulate(optimizationSettings.optimizationVariablesGroupsWeights->begin(), + optimizationSettings.optimizationVariablesGroupsWeights->end(), + 0))); + + for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) { + //Increase the size of the vector holding the simulation scenario indices + global.simulationScenarioIndices.resize( + global.simulationScenarioIndices.size() + + simulationScenariosResolution[baseSimulationScenarioIndex]); + //Add the simulation scenarios indices that correspond to this base simulation scenario + std::iota(global.simulationScenarioIndices.end() + - simulationScenariosResolution[baseSimulationScenarioIndex], + global.simulationScenarioIndices.end(), + std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + + baseSimulationScenarioIndex, + 0)); + } + + // global.totalNumberOfSimulationScenarios = totalNumberOfSimulationScenarios; + global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios); + global.fullPatternResults.resize(totalNumberOfSimulationScenarios); + global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios); + global.rotationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios); + global.minY = std::numeric_limits::max(); + global.numOfSimulationCrashes = 0; + global.numberOfFunctionCalls = 0; + global.optimizationSettings = optimizationSettings; + global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; + +#ifdef POLYSCOPE_DEFINED + constexpr bool recomputeMagnitudes = false; +#else + constexpr bool recomputeMagnitudes = true; +#endif + // std::array fullPatternSimulationScenarioMaxMagnitudes + // = getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices, + // intermediateResultsDirectoryPath, + // recomputeMagnitudes); + std::array fullPatternSimulationScenarioMaxMagnitudes + = optimizationSettings.baseScenarioMagnitudes; + global.fullPatternSimulationJobs + = createFullPatternSimulationJobs(m_pFullPatternSimulationMesh, + fullPatternSimulationScenarioMaxMagnitudes); + // polyscope::removeAllStructures(); + + // ComputeScenarioWeights({Axial, Shear, Bending, Dome, Saddle}, + computeScenarioWeights(desiredBaseSimulationScenarioIndices); + + results.baseTriangle = global.baseTriangle; + + DRMSimulationModel::Settings simulationSettings; + simulationSettings.totalExternalForcesNormPercentageTermination = 1e-5; +// simulationSettings.maxDRMIterations = 200000; +// simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-10; +// simulationSettings.viscousDampingFactor = 5e-3; +// simulationSettings.useKineticDamping = true; +// simulationSettings.save(std::filesystem::path(std::string(__FILE__)) +// .parent_path() +// .parent_path() +// .append("DefaultSettings") +// .append("DRMSettings") +// .append("defaultDRMSimulationSettings.json")); + +// simulationSettings.averageResidualForcesCriterionThreshold = 1e-5; +// simulationSettings.viscousDampingFactor = 1e-3; +// simulationSettings.beVerbose = true; +// simulationSettings.shouldDraw = true; +// simulationSettings.isDebugMode = true; +// simulationSettings.debugModeStep = 100000; +#ifdef POLYSCOPE_DEFINED + constexpr bool drawFullPatternSimulationResults = false; + if (drawFullPatternSimulationResults) { + global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( + ReducedModelOptimization::Colors::fullInitial); + } +#endif + results.wasSuccessful = true; + for (int simulationScenarioIndex : global.simulationScenarioIndices) { + const std::shared_ptr &pFullPatternSimulationJob + = global.fullPatternSimulationJobs[simulationScenarioIndex]; + + std::filesystem::path jobResultsDirectoryPath( + std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) + .append("FullPatternResults") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append(pFullPatternSimulationJob->getLabel())); + // .append(pFullPatternSimulationJob->getLabel() + ".json") +#ifdef POLYSCOPE_DEFINED + constexpr bool recomputeFullPatternResults = false; +#else + constexpr bool recomputeFullPatternResults = true; +#endif + SimulationResults fullPatternResults; + if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) { + fullPatternResults.load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), + pFullPatternSimulationJob); + } else { + // const bool fullPatternScenarioMagnitudesExist = std::filesystem::exists( + // patternMaxForceMagnitudesFilePath); + // if (fullPatternScenarioMagnitudesExist) { + // nlohmann::json json; + // std::ifstream ifs(patternMaxForceMagnitudesFilePath.string()); + // ifs >> json; + // fullPatternSimulationScenarioMaxMagnitudes + // = static_cast>>( + // json.at("maxMagn")); + // const bool shouldRecompute = fullPatternSimulationScenarioMaxMagnitudes.size() + // != desiredBaseSimulationScenarioIndices.size(); + // if (!shouldRecompute) { + // return fullPatternSimulationScenarioMaxMagnitudes; + // } + // } + //#ifdef POLYSCOPE_DEFINED + // LinearSimulationModel linearSimulator; + // SimulationResults fullPatternResults = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); + + //#else + fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, + simulationSettings); + fullPatternResults.save(jobResultsDirectoryPath); + } +//#endif +// if (!fullPatternResults.converged) { +// DRMSimulationModel::Settings simulationSettings_secondRound; +// simulationSettings_secondRound.viscousDampingFactor = 2e-3; +// simulationSettings_secondRound.useKineticDamping = true; +// simulationSettings.maxDRMIterations = 200000; +// fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, +// simulationSettings_secondRound); +#ifdef POLYSCOPE_DEFINED +// std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel() +// << " used viscous damping." << std::endl; +#endif + + if (!fullPatternResults.converged) { + results.wasSuccessful = false; + std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() + << " did not converge." << std::endl; +#ifdef POLYSCOPE_DEFINED + // DRMSimulationModel::Settings debugSimulationSettings; + // debugSimulationSettings.debugModeStep = 50; + // // // debugSimulationSettings.maxDRMIterations = 100000; + // debugSimulationSettings.shouldDraw = true; + // // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; + // debugSimulationSettings.shouldCreatePlots = true; + // // // debugSimulationSettings.Dtini = 0.06; + // debugSimulationSettings.beVerbose = true; + // debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5; + // debugSimulationSettings.maxDRMIterations = 100000; + // debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; + // debugSimulationSettings.viscousDampingFactor = 1e-2; + // // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; + // // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; + // auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob, + // debugSimulationSettings); + // debugResults.setLabelPrefix("debugResults"); + // debugResults.registerForDrawing(); + // polyscope::show(); + // debugResults.unregister(); + std::filesystem::path outputPath( + std::filesystem::path("../nonConvergingJobs") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append("final_" + pFullPatternSimulationJob->getLabel())); + std::filesystem::create_directories(outputPath); + pFullPatternSimulationJob->save(outputPath); + simulationSettings.save(outputPath); + std::terminate(); + return; +#endif + // } + } +#ifdef POLYSCOPE_DEFINED + if (drawFullPatternSimulationResults) { + // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); + fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + true); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + // true, + // true); + polyscope::show(); + fullPatternResults.unregister(); + // fullPatternResults_linear.unregister(); + } +#endif + global.fullPatternResults[simulationScenarioIndex] = fullPatternResults; + SimulationJob reducedPatternSimulationJob; + reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; + computeReducedModelSimulationJob(*pFullPatternSimulationJob, + m_fullToReducedInterfaceViMap, + reducedPatternSimulationJob); + global.reducedPatternSimulationJobs[simulationScenarioIndex] + = std::make_shared(reducedPatternSimulationJob); + // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; + } + +#ifdef POLYSCOPE_DEFINED + if (drawFullPatternSimulationResults) { + global.fullPatternSimulationJobs[0]->pMesh->unregister(); + } +#endif + // if (global.optimizationSettings.normalizationStrategy + // == Settings::NormalizationStrategy::Epsilon) { + computeObjectiveValueNormalizationFactors(); +// } +#ifdef POLYSCOPE_DEFINED + std::cout << "Running reduced model optimization" << std::endl; +#endif + runOptimization(optimizationSettings, results); + results.notes = optimizationNotes; +} diff --git a/reducedmodeloptimizer.hpp b/reducedmodeloptimizer.hpp new file mode 100644 index 0000000..93a4885 --- /dev/null +++ b/reducedmodeloptimizer.hpp @@ -0,0 +1,307 @@ +#ifndef REDUCEDMODELOPTIMIZER_HPP +#define REDUCEDMODELOPTIMIZER_HPP + +#include "csvfile.hpp" +#include "drmsimulationmodel.hpp" +#include "edgemesh.hpp" +#include "linearsimulationmodel.hpp" +#ifdef POLYSCOPE_DEFINED +#include "matplot/matplot.h" +#endif +#include "reducedmodel.hpp" +#include "reducedmodeloptimizer_structs.hpp" +#include "simulationmesh.hpp" +#include +#ifdef DLIB_DEFINED +#include +#include +#endif + +#ifdef POLYSCOPE_DEFINED +#include "polyscope/color_management.h" +#endif // POLYSCOPE_DEFINED +using FullPatternVertexIndex = VertexIndex; +using ReducedPatternVertexIndex = VertexIndex; + +class ReducedModelOptimizer +{ + std::shared_ptr m_pReducedPatternSimulationMesh; + std::shared_ptr m_pFullPatternSimulationMesh; + std::unordered_map + m_fullToReducedInterfaceViMap; + std::vector> + m_fullPatternOppositeInterfaceViPairs; + std::unordered_map nodeToSlot; + std::unordered_map> slotToNode; + std::string optimizationNotes; + std::vector> &, + SimulationJob &)>> + constructBaseScenarioFunctions; + std::vector scenarioIsSymmetrical; + int fullPatternNumberOfEdges; + constexpr static double youngsModulus{1 * 1e9}; + std::filesystem::path intermediateResultsDirectoryPath; + std::string fullPatternLabel; + +public: + struct FunctionEvaluation + { + FunctionEvaluation() = default; + FunctionEvaluation(const std::vector &x, double y) : x(x), y(y) {} + + std::vector x; + double y = std::numeric_limits::quiet_NaN(); + }; + + // 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 th= {"Theta"}; + // inline const static std::string R = {"R"}; + // }; + // inline constexpr static ParameterLabels parameterLabels(); + + inline static std::array + parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"}; + constexpr static std::array + simulationScenariosResolution = {11, 11, 20, 20, 20}; + constexpr static std::array + baseScenarioWeights = {1, 1, 2, 3, 2}; + inline static int totalNumberOfSimulationScenarios + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.end(), + 0); + inline static int fanSize{6}; + inline static double initialHexagonSize{0.3}; + int interfaceNodeIndex; + double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; + + ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); + static void computeReducedModelSimulationJob( + const SimulationJob &simulationJobOfFullModel, + const std::unordered_map &fullToReducedMap, + SimulationJob &simulationJobOfReducedModel); + + SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); + + static void runSimulation(const std::string &filename, std::vector &x); + + static std::vector> createFullPatternSimulationJobs( + const std::shared_ptr &pMesh, + const std::unordered_map &fullPatternOppositeInterfaceViMap); + + static void createSimulationMeshes( + PatternGeometry &fullModel, + PatternGeometry &reducedModel, + std::shared_ptr &pFullPatternSimulationMesh, + std::shared_ptr &pReducedPatternSimulationMesh); + static void computeMaps(const std::unordered_map> &slotToNode, + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + std::unordered_map + &reducedToFullInterfaceViMap, + std::unordered_map + &fullToReducedInterfaceViMap, + std::vector> + &fullPatternOppositeInterfaceViMap); + static void visualizeResults( + const std::vector> &fullPatternSimulationJobs, + const std::vector> &reducedPatternSimulationJobs, + const std::vector &simulationScenarios, + const std::unordered_map + &reducedToFullInterfaceViMap); + static void registerResultsForDrawing( + const std::shared_ptr &pFullPatternSimulationJob, + const std::shared_ptr &pReducedPatternSimulationJob, + const std::unordered_map + &reducedToFullInterfaceViMap); + + static double computeRawTranslationalError( + const std::vector &fullPatternDisplacements, + const std::vector &reducedPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap); + + static double computeDisplacementError( + const std::vector &fullPatternDisplacements, + const std::vector &reducedPatternDisplacements, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor); + + static double objective(double E, + double A, + double innerHexagonSize, + double innerHexagonRotationAngle); + static double computeRawRotationalError( + const std::vector &rotatedQuaternion_fullPattern, + const std::vector &rotatedQuaternion_reducedPattern, + const std::unordered_map + &reducedToFullInterfaceViMap); + + static double computeRotationalError(const std::vector &rotatedQuaternion_fullPattern, + const std::vector &rotatedQuaternion_reducedPattern, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor); + static double computeError( + const SimulationResults &simulationResults_fullPattern, + const SimulationResults &simulationResults_reducedPattern, + const std::unordered_map + &reducedToFullInterfaceViMap, + const double &normalizationFactor_translationalDisplacement, + const double &normalizationFactor_rotationalDisplacement, + const double &scenarioWeight, + const ReducedModelOptimization::Settings::ObjectiveWeights &objectiveWeights); + static void constructAxialSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job); + + static void constructShearSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job); + + static void constructBendingSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job); + + static void constructDomeSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job); + + static void constructSaddleSimulationScenario( + const double &forceMagnitude, + const std::vector> + &oppositeInterfaceViPairs, + SimulationJob &job); + static std::function &x, + std::shared_ptr &pReducedPatternSimulationMesh)> + function_updateReducedPattern; + static std::function &pReducedPatternSimulationMesh)> + function_updateReducedPattern_material_E; + static std::function &pReducedPatternSimulationMesh)> + function_updateReducedPattern_material_A; + static std::function &pReducedPatternSimulationMesh)> + function_updateReducedPattern_material_I; + static std::function &pReducedPatternSimulationMesh)> + function_updateReducedPattern_material_I2; + static std::function &pReducedPatternSimulationMesh)> + function_updateReducedPattern_material_I3; + static std::function &pReducedPatternSimulationMesh)> + function_updateReducedPattern_material_J; + static double objective(const std::vector &x); + static void initializeUpdateReducedPatternFunctions(); + static double objective(const double &xValue); + + ReducedModelOptimization::Results optimize( + ConstPatternGeometry &fullPattern, + const ReducedModelOptimization::Settings &optimizationSettings); + + void optimize(ConstPatternGeometry &fullPattern, + const ReducedModelOptimization::Settings &optimizationSettings, + ReducedModelOptimization::Results &optimizationResults); + +private: + void optimize( + const ReducedModelOptimization::Settings &xRanges, + ReducedModelOptimization::Results &results, + const std::vector &simulationScenarios + = std::vector( + {ReducedModelOptimization::Axial, + ReducedModelOptimization::Shear, + ReducedModelOptimization::Bending, + ReducedModelOptimization::Dome, + ReducedModelOptimization::Saddle})); + + void initializePatterns(PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, + const std::array + &optimizationParameters); + static void computeDesiredReducedModelDisplacements( + const SimulationResults &fullModelResults, + const std::unordered_map &displacementsReducedToFullMap, + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); + static void runOptimization(const ReducedModelOptimization::Settings &settings, + ReducedModelOptimization::Results &results); + void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); + void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); + static void initializeOptimizationParameters( + const std::shared_ptr &mesh, + const std::array + &optimizationParamters); + + DRMSimulationModel simulator; + void computeObjectiveValueNormalizationFactors(); + + static void getResults(const FunctionEvaluation &optimalObjective, + const ReducedModelOptimization::Settings &settings, + ReducedModelOptimization::Results &results); + double computeFullPatternMaxSimulationForce( + const ReducedModelOptimization::BaseSimulationScenario &scenario) const; + +#ifdef DLIB_DEFINED + static double objective(const dlib::matrix &x); +#endif + void computeScenarioWeights(const std::vector + &baseSimulationScenarios); + std::array + computeFullPatternMaxSimulationForces( + const std::vector + &desiredBaseSimulationScenario) const; + std::vector> createFullPatternSimulationJobs( + const std::shared_ptr &pMesh, + const std::array + &baseScenarioMaxForceMagnitudes) const; + std::array + getFullPatternMaxSimulationForces( + const std::vector + &desiredBaseSimulationScenarioIndices, + const std::filesystem::path &intermediateResultsDirectoryPath, + const bool &recomputeForceMagnitudes); + std::array + getFullPatternMaxSimulationForces(); +}; +inline std::function &pReducedPatternSimulationMesh)> + ReducedModelOptimizer::function_updateReducedPattern_material_E; +inline std::function &pReducedPatternSimulationMesh)> + ReducedModelOptimizer::function_updateReducedPattern_material_A; +inline std::function &pReducedPatternSimulationMesh)> + ReducedModelOptimizer::function_updateReducedPattern_material_I; +inline std::function &pReducedPatternSimulationMesh)> + ReducedModelOptimizer::function_updateReducedPattern_material_I2; +inline std::function &pReducedPatternSimulationMesh)> + ReducedModelOptimizer::function_updateReducedPattern_material_I3; +inline std::function &pReducedPatternSimulationMesh)> + ReducedModelOptimizer::function_updateReducedPattern_material_J; +inline std::function &x, std::shared_ptr &m)> + ReducedModelOptimizer::function_updateReducedPattern; + +#endif // REDUCEDMODELOPTIMIZER_HPP