diff --git a/src/main.cpp b/src/main.cpp index 5f10e44..c717503 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -40,21 +40,21 @@ int main(int argc, char *argv[]) { reducedPattern.scale(0.03, interfaceNodeIndex); // Set the optization settings - ReducedModelOptimization::xRange beamE{"E", 0.001, 1000}; - ReducedModelOptimization::xRange beamA{"A", 0.001, 1000}; - ReducedModelOptimization::xRange beamI2{"I2", 0.001, 1000}; - ReducedModelOptimization::xRange beamI3{"I3", 0.001, 1000}; - ReducedModelOptimization::xRange beamJ{"J", 0.001, 1000}; - ReducedModelOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95}; - ReducedModelOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0}; - ReducedModelOptimization::Settings settings_optimization; + ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000}; + ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000}; + ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000}; + ReducedPatternOptimization::xRange beamI3{"I3", 0.001, 1000}; + ReducedPatternOptimization::xRange beamJ{"J", 0.001, 1000}; + ReducedPatternOptimization::xRange innerHexagonSize{"HexSize", 0.05, 0.95}; + ReducedPatternOptimization::xRange innerHexagonAngle{"HexAngle", -30.0, 30.0}; + ReducedPatternOptimization::Settings settings_optimization; settings_optimization.xRanges = {beamE,beamA,beamJ,beamI2,beamI3, innerHexagonSize, innerHexagonAngle}; const bool input_numberOfFunctionCallsDefined = argc >= 4; settings_optimization.numberOfFunctionCalls = input_numberOfFunctionCallsDefined ? std::atoi(argv[3]) : 100; settings_optimization.normalizationStrategy - = ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon; + = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; settings_optimization.normalizationParameter = 0.0003; settings_optimization.solutionAccuracy = 0.001; @@ -66,7 +66,8 @@ int main(int argc, char *argv[]) { assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]); ReducedModelOptimizer optimizer(numberOfNodesPerSlot); optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size()); - ReducedModelOptimization::Results optimizationResults = optimizer.optimize(settings_optimization); + ReducedPatternOptimization::Results optimizationResults = optimizer.optimize( + settings_optimization); // Export results const bool input_resultDirectoryDefined = argc >= 5; @@ -105,7 +106,9 @@ int main(int argc, char *argv[]) { csv_results << endrow; #ifdef POLYSCOPE_DEFINED + optimizationResults.saveMeshFiles(); optimizationResults.draw(); #endif + return 0; } diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index ede5cae..4673b45 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -4,10 +4,8 @@ #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" #include -#include -#include -using namespace ReducedModelOptimization; +using namespace ReducedPatternOptimization; struct GlobalOptimizationVariables { // std::vector> fullPatternDisplacements; @@ -22,15 +20,12 @@ struct GlobalOptimizationVariables { std::vector objectiveValueHistory; Eigen::VectorXd initialParameters; std::vector simulationScenarioIndices; - std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; - double innerHexagonInitialRotationAngle{30}; - double innerHexagonInitialPos{0}; double minY{DBL_MAX}; std::vector minX; int numOfSimulationCrashes{false}; int numberOfFunctionCalls{0}; int numberOfOptimizationParameters{5}; - ReducedModelOptimization::Settings optimizationSettings; + ReducedPatternOptimization::Settings optimizationSettings; vcg::Triangle3 baseTriangle; } global; @@ -41,13 +36,13 @@ double ReducedModelOptimizer::computeDisplacementError( &reducedToFullInterfaceViMap, const double &normalizationFactor) { - const double rawError = computeRawDisplacementError(fullPatternDisplacements, - reducedPatternDisplacements, - reducedToFullInterfaceViMap); + const double rawError = computeRawTranslationalError(fullPatternDisplacements, + reducedPatternDisplacements, + reducedToFullInterfaceViMap); return rawError / normalizationFactor; } -double ReducedModelOptimizer::computeRawDisplacementError( +double ReducedModelOptimizer::computeRawTranslationalError( const std::vector &fullPatternDisplacements, const std::vector &reducedPatternDisplacements, const std::unordered_map @@ -56,18 +51,14 @@ double ReducedModelOptimizer::computeRawDisplacementError( double error = 0; for (const auto reducedFullViPair : reducedToFullInterfaceViMap) { const VertexIndex reducedModelVi = reducedFullViPair.first; - Eigen::Vector3d reducedVertexDisplacement( - reducedPatternDisplacements[reducedModelVi][0], - reducedPatternDisplacements[reducedModelVi][1], - reducedPatternDisplacements[reducedModelVi][2]); const VertexIndex fullModelVi = reducedFullViPair.second; - Eigen::Vector3d fullVertexDisplacement( - fullPatternDisplacements[fullModelVi][0], - fullPatternDisplacements[fullModelVi][1], - fullPatternDisplacements[fullModelVi][2]); - Eigen::Vector3d errorVector = - fullVertexDisplacement - reducedVertexDisplacement; - error += errorVector.norm(); + const Eigen::Vector3d fullPatternVertexDiplacement = fullPatternDisplacements[fullModelVi] + .getTranslation(); + const Eigen::Vector3d reducedPatternVertexDiplacement + = reducedPatternDisplacements[reducedModelVi].getTranslation(); + const double vertexError = (fullPatternVertexDiplacement - reducedPatternVertexDiplacement) + .norm(); + error += vertexError; } return error; @@ -121,7 +112,7 @@ double ReducedModelOptimizer::computeError( simulationResults_reducedPattern.rotationalDisplacementQuaternion, reducedToFullInterfaceViMap, normalizationFactor_rotationalDisplacement); - return translationalError + rotationalError; + return 1.5 * translationalError + 0.5 * rotationalError; } double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, @@ -226,10 +217,10 @@ double ReducedModelOptimizer::objective(long n, const double *x) { global.minX.assign(x, x + n); } #ifdef POLYSCOPE_DEFINED -// if (++global.numberOfFunctionCalls % 100 == 0) { -// std::cout << "Number of function calls:" << global.numberOfFunctionCalls -// << std::endl; -// } + ++global.numberOfFunctionCalls; + if (global.numberOfFunctionCalls % 100 == 0) { + std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; + } #endif // compute error and return it // global.objectiveValueHistory.push_back(totalError); @@ -517,9 +508,6 @@ void ReducedModelOptimizer::initializeOptimizationParameters( global.initialParameters(2) = mesh->elements[0].J; global.initialParameters(3) = mesh->elements[0].I2; global.initialParameters(4) = mesh->elements[0].I3; - global.initialParameters(optimizationParamters-2) = global.innerHexagonInitialPos; - global.innerHexagonInitialRotationAngle = 30; - global.initialParameters(optimizationParamters-1) = global.innerHexagonInitialRotationAngle; } void ReducedModelOptimizer::computeReducedModelSimulationJob( @@ -674,89 +662,105 @@ void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( } } -ReducedModelOptimization::Results -ReducedModelOptimizer::runOptimization(const Settings &settings) { +ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( + const dlib::function_evaluation &optimizationResult_dlib, const Settings &settings) +{ + ReducedPatternOptimization::Results results; - global.objectiveValueHistory.clear(); - dlib::matrix xMin(global.numberOfOptimizationParameters); - dlib::matrix xMax(global.numberOfOptimizationParameters); - for (int i = 0; i < global.numberOfOptimizationParameters; i++) { - xMin(i) = settings.xRanges[i].min; - xMax(i) = settings.xRanges[i].max; - } - - auto start = std::chrono::system_clock::now(); - dlib::function_evaluation result; - double (*objF)(double, double, double, double, double,double,double) = &objective; - result = dlib::find_min_global( - objF, xMin, xMax, - dlib::max_function_calls(settings.numberOfFunctionCalls), - std::chrono::hours(24 * 365 * 290), settings.solutionAccuracy); - auto end = std::chrono::system_clock::now(); - auto elapsed = - std::chrono::duration_cast(end - start); - ReducedModelOptimization::Results results; + results.baseTriangle = global.baseTriangle; + //Number of crashes results.numberOfSimulationCrashes = global.numOfSimulationCrashes; + //Value of optimal objective Y + results.objectiveValue.total = optimizationResult_dlib.y; + //Optimal X values results.optimalXNameValuePairs.resize(settings.xRanges.size()); - for(int xVariableIndex=0;xVariableIndex optimalX(settings.xRanges.size()); + for (int xVariableIndex = 0; xVariableIndex < settings.xRanges.size(); xVariableIndex++) { + if (xVariableIndex < 5) { + results.optimalXNameValuePairs[xVariableIndex] + = std::make_pair(settings.xRanges[xVariableIndex].label, + global.minX[xVariableIndex] + * global.initialParameters(xVariableIndex)); + } else { + //Hex size and angle are pure values (not multipliers of the initial values) + results.optimalXNameValuePairs[xVariableIndex] + = std::make_pair(settings.xRanges[xVariableIndex].label, + global.minX[xVariableIndex]); + } + + assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex)); + optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex); } // Compute obj value per simulation scenario and the raw objective value - results.rawObjectiveValue=0; - std::vector optimalX(results.optimalXNameValuePairs.size()); - for (int xVariableIndex = 0; xVariableIndex < global.numberOfOptimizationParameters; - xVariableIndex++) { - optimalX[xVariableIndex] = global.minX[xVariableIndex]; - assert(global.minX[xVariableIndex] == result.x(xVariableIndex)); - } updateMesh(optimalX.size(), optimalX.data()); - results.objectiveValuePerSimulationScenario.resize(totalNumberOfSimulationScenarios); + // results.objectiveValue.totalPerSimulationScenario.resize(totalNumberOfSimulationScenarios); LinearSimulationModel simulator; - double totalObjectiveValue = 0; - double totalRawObjectiveValue = 0; + results.objectiveValue.totalRaw = 0; + results.objectiveValue.perSimulationScenario_translational.resize( + totalNumberOfSimulationScenarios); + results.objectiveValue.perSimulationScenario_rawTranslational.resize( + totalNumberOfSimulationScenarios); + results.objectiveValue.perSimulationScenario_rotational.resize(totalNumberOfSimulationScenarios); + results.objectiveValue.perSimulationScenario_rawRotational.resize( + totalNumberOfSimulationScenarios); + results.objectiveValue.perSimulationScenario_total.resize(totalNumberOfSimulationScenarios); for (int simulationScenarioIndex : global.simulationScenarioIndices) { SimulationResults reducedModelResults = simulator.executeSimulation( global.reducedPatternSimulationJobs[simulationScenarioIndex]); - totalObjectiveValue += computeError( + results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex] = computeError( global.fullPatternResults[simulationScenarioIndex], reducedModelResults, global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex], global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); - const double rawTranslationalError = computeRawDisplacementError( + //Raw translational + const double rawTranslationalError = computeRawTranslationalError( global.fullPatternResults[simulationScenarioIndex].displacements, reducedModelResults.displacements, global.reducedToFullInterfaceViMap); + results.objectiveValue.perSimulationScenario_rawTranslational[simulationScenarioIndex] + = rawTranslationalError; + //Raw rotational const double rawRotationalError = computeRawRotationalError( global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, reducedModelResults.rotationalDisplacementQuaternion, global.reducedToFullInterfaceViMap); + results.objectiveValue.perSimulationScenario_rawRotational[simulationScenarioIndex] + = rawRotationalError; - results.rawObjectiveValue += rawTranslationalError + rawRotationalError; - - const double normalizedTranslationalError - = rawTranslationalError - / global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; - const double test_normalizedTranslationError = computeDisplacementError( + //Normalized translational + const double normalizedTranslationalError = computeDisplacementError( global.fullPatternResults[simulationScenarioIndex].displacements, reducedModelResults.displacements, global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); - const double normalizedRotationalError - = rawRotationalError - / global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]; - const double test_normalizedRotationalError = computeRotationalError( + results.objectiveValue.perSimulationScenario_translational[simulationScenarioIndex] + = 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]); - assert(test_normalizedTranslationError == normalizedTranslationalError); - assert(test_normalizedRotationalError == normalizedRotationalError); - std::cout << "Simulation scenario:" << baseSimulationScenarioNames[simulationScenarioIndex] + results.objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex] + = 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); + std::cout << "Simulation scenario:" + << global.reducedPatternSimulationJobs[simulationScenarioIndex]->getLabel() << std::endl; std::cout << "raw translational error:" << rawTranslationalError << std::endl; std::cout << "translation normalization value:" @@ -768,180 +772,307 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { << global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] << std::endl; std::cout << "Rotational error:" << normalizedRotationalError << std::endl; - results.objectiveValuePerSimulationScenario[simulationScenarioIndex] - = normalizedTranslationalError + normalizedRotationalError; + // results.objectiveValuePerSimulationScenario[simulationScenarioIndex] + // = normalizedTranslationalError + normalizedRotationalError; std::cout << "Objective value:" - << results.objectiveValuePerSimulationScenario[simulationScenarioIndex] + << results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex] << std::endl; - // totalObjectiveValue += results.objectiveValuePerSimulationScenario[simulationScenarioIndex]; - totalRawObjectiveValue += rawTranslationalError + rawRotationalError; + results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError; std::cout << std::endl; } - assert(result.y == global.minY); - results.objectiveValue = result.y; - auto functionReturnValue = objective(optimalX.size(), optimalX.data()); - std::cout << "Total (functionReturn)" << functionReturnValue << std::endl; - std::cout << "Total:" << totalObjectiveValue << std::endl; -#ifdef POLYSCOPE_DEFINED - std::cout << "Total optimal objective value:" << results.objectiveValue << std::endl; - // std::cout << "Total raw objective value:" << global.minY << std::endl; - if (global.minY != result.y) { - std::cerr << "Global min objective is not equal to result objective" << std::endl; - } -#endif - results.time = elapsed.count() / 1000.0; const bool printDebugInfo = false; if (printDebugInfo) { std::cout << "Finished optimizing." << endl; - // std::cout << "Solution x:" << endl; - // std::cout << result.x << endl; - std::cout << "Objective value:" << global.minY << endl; + std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl; + assert(global.minY == optimizationResult_dlib.y); + } + + 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); } return results; } +ReducedPatternOptimization::Results ReducedModelOptimizer::runOptimization(const Settings &settings) +{ + global.objectiveValueHistory.clear(); + dlib::matrix xMin(global.numberOfOptimizationParameters); + dlib::matrix xMax(global.numberOfOptimizationParameters); + for (int i = 0; i < global.numberOfOptimizationParameters; i++) { + xMin(i) = settings.xRanges[i].min; + xMax(i) = settings.xRanges[i].max; + } + + auto start = std::chrono::system_clock::now(); + dlib::function_evaluation result_dlib; + double (*objF)(double, double, double, double, double,double,double) = &objective; + result_dlib = dlib::find_min_global(objF, + xMin, + xMax, + dlib::max_function_calls(settings.numberOfFunctionCalls), + std::chrono::hours(24 * 365 * 290), + settings.solutionAccuracy); + auto end = std::chrono::system_clock::now(); + auto elapsed = + std::chrono::duration_cast(end - start); + ReducedPatternOptimization::Results results = getResults(result_dlib, settings); + results.time = elapsed.count() / 1000.0; + return results; +} + +//TODO: create a function that takes as arguments the magnitude range, a lambda that generates the simulation job and the base sim scenario index std::vector> -ReducedModelOptimizer::createScenarios( - const std::shared_ptr &pMesh) { +ReducedModelOptimizer::createFullPatternSimulationScenarios( + const std::shared_ptr &pMesh) +{ std::vector> scenarios; scenarios.resize(totalNumberOfSimulationScenarios); std::unordered_map> fixedVertices; std::unordered_map nodalForces; //// Axial - int simulationScenarioIndex = BaseSimulationScenario::Axial; - for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); - viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); - viPairIt++) { - if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) { - CoordType forceDirection(1, 0, 0); - const auto viPair = *viPairIt; - nodalForces[viPair.first] - = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) - * 500; - fixedVertices[viPair.second] = - std::unordered_set{0, 1, 2, 3, 4, 5}; + const double maxForceMagnitude_axial = 500; + const double minForceMagnitude_axial = -500; + const int numberOfSimulationScenarios_axial + = simulationScenariosResolution[BaseSimulationScenario::Axial]; + const double forceMagnitudeStep_axial = numberOfSimulationScenarios_axial == 1 + ? maxForceMagnitude_axial + : (maxForceMagnitude_axial - minForceMagnitude_axial) + / (numberOfSimulationScenarios_axial - 1); + const int baseSimulationScenarioIndexOffset_axial + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + BaseSimulationScenario::Axial, + 0); + + for (int axialSimulationScenarioIndex = 0; + axialSimulationScenarioIndex < numberOfSimulationScenarios_axial; + axialSimulationScenarioIndex++) { + const double forceMagnitude = (forceMagnitudeStep_axial * axialSimulationScenarioIndex + + minForceMagnitude_axial); + for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); + viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); + viPairIt++) { + if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) { + CoordType forceDirection(1, 0, 0); + const auto viPair = *viPairIt; + nodalForces[viPair.first] + = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) + * forceMagnitude; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } } + + scenarios[baseSimulationScenarioIndexOffset_axial + axialSimulationScenarioIndex] + = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_" + + std::to_string(axialSimulationScenarioIndex), + fixedVertices, + nodalForces, + {})); } - scenarios[simulationScenarioIndex] = std::make_shared( - SimulationJob(pMesh, - "Axial_" + std::to_string(simulationScenarioIndex), - fixedVertices, - nodalForces, - {})); //// Shear - simulationScenarioIndex = BaseSimulationScenario::Shear; - fixedVertices.clear(); - nodalForces.clear(); - // NewMethod - for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); - viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); - viPairIt++) { - if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) { - CoordType forceDirection(0, 1, 0); - const auto viPair = *viPairIt; - nodalForces[viPair.first] - = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) * 40; - fixedVertices[viPair.second] = - std::unordered_set{0, 1, 2, 3, 4, 5}; + const double maxForceMagnitude_shear = 50; + const double minForceMagnitude_shear = -50; + const int numberOfSimulationScenarios_shear + = simulationScenariosResolution[BaseSimulationScenario::Shear]; + const double forceMagnitudeStep_shear = numberOfSimulationScenarios_shear == 1 + ? maxForceMagnitude_shear + : (maxForceMagnitude_shear - minForceMagnitude_shear) + / (numberOfSimulationScenarios_shear - 1); + const int baseSimulationScenarioIndexOffset_shear + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + BaseSimulationScenario::Shear, + 0); + for (int shearSimulationScenarioIndex = 0; + shearSimulationScenarioIndex < numberOfSimulationScenarios_shear; + shearSimulationScenarioIndex++) { + fixedVertices.clear(); + nodalForces.clear(); + const double forceMagnitude = (forceMagnitudeStep_shear * shearSimulationScenarioIndex + + minForceMagnitude_shear); + for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); + viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); + viPairIt++) { + if (viPairIt != m_fullPatternOppositeInterfaceViPairs.begin()) { + CoordType forceDirection(0, 1, 0); + const auto viPair = *viPairIt; + nodalForces[viPair.first] + = Vector6d({forceDirection[0], forceDirection[1], forceDirection[2], 0, 0, 0}) + * forceMagnitude; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } } + scenarios[baseSimulationScenarioIndexOffset_shear + shearSimulationScenarioIndex] + = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_" + + std::to_string(shearSimulationScenarioIndex), + fixedVertices, + nodalForces, + {})); } - scenarios[simulationScenarioIndex] = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[simulationScenarioIndex], - fixedVertices, - nodalForces, - {})); //// Bending - simulationScenarioIndex = BaseSimulationScenario::Bending; - fixedVertices.clear(); - nodalForces.clear(); - for (const auto &viPair : m_fullPatternOppositeInterfaceViPairs) { - nodalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * 0.005; - fixedVertices[viPair.second] = - std::unordered_set{0, 1, 2, 3, 4, 5}; + const double maxForceMagnitude_bending = 0.005; + const double minForceMagnitude_bending = -0.005; + const int numberOfSimulationScenarios_bending + = simulationScenariosResolution[BaseSimulationScenario::Bending]; + const double forceMagnitudeStep_bending = numberOfSimulationScenarios_bending == 1 + ? maxForceMagnitude_bending + : (maxForceMagnitude_bending + - minForceMagnitude_bending) + / (numberOfSimulationScenarios_bending - 1); + const int baseSimulationScenarioIndexOffset_bending + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + BaseSimulationScenario::Bending, + 0); + for (int bendingSimulationScenarioIndex = 0; + bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending; + bendingSimulationScenarioIndex++) { + fixedVertices.clear(); + nodalForces.clear(); + const double forceMagnitude = (forceMagnitudeStep_bending * bendingSimulationScenarioIndex + + minForceMagnitude_bending); + for (const auto &viPair : m_fullPatternOppositeInterfaceViPairs) { + nodalForces[viPair.first] = Vector6d({0, 0, 1, 0, 0, 0}) * forceMagnitude; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2, 3, 4, 5}; + } + scenarios[baseSimulationScenarioIndexOffset_bending + bendingSimulationScenarioIndex] + = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_" + + std::to_string(bendingSimulationScenarioIndex), + fixedVertices, + nodalForces, + {})); } - scenarios[simulationScenarioIndex] = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[simulationScenarioIndex], - fixedVertices, - nodalForces, - {})); //// Dome - simulationScenarioIndex = BaseSimulationScenario::Dome; - fixedVertices.clear(); - nodalForces.clear(); - std::unordered_map nodalForcedDisplacements; - for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); - viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); - viPairIt++) { - const auto viPair = *viPairIt; - if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) { - CoordType interfaceV = (pMesh->vert[viPair.first].cP() - - pMesh->vert[viPair.second].cP()); - nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceV[0], - -interfaceV[1], - 0) - * 0.025; - nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceV[0], - interfaceV[1], - 0) - * 0.025; - // 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 { - fixedVertices[viPair.first] = std::unordered_set{2}; - fixedVertices[viPair.second] = std::unordered_set{2}; - 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}) * 0.005; - nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.005; + const double maxForceMagnitude_dome = 0.025; + const double minForceMagnitude_dome = -0.025; + const int numberOfSimulationScenarios_dome + = simulationScenariosResolution[BaseSimulationScenario::Dome]; + const double forceMagnitudeStep_dome = numberOfSimulationScenarios_dome == 1 + ? maxForceMagnitude_dome + : (maxForceMagnitude_dome - minForceMagnitude_dome) + / (numberOfSimulationScenarios_dome - 1); + const int baseSimulationScenarioIndexOffset_dome + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + BaseSimulationScenario::Dome, + 0); + for (int domeSimulationScenarioIndex = 0; + domeSimulationScenarioIndex < numberOfSimulationScenarios_dome; + domeSimulationScenarioIndex++) { + fixedVertices.clear(); + nodalForces.clear(); + std::unordered_map nodalForcedDisplacements; + const double forceMagnitude = (forceMagnitudeStep_dome * domeSimulationScenarioIndex + + minForceMagnitude_dome); + for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); + viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); + viPairIt++) { + const auto viPair = *viPairIt; + CoordType interfaceVector = (pMesh->vert[viPair.first].cP() + - pMesh->vert[viPair.second].cP()); + VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0)) + * interfaceVector.Normalize(); + if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) { + nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceVector[0], + -interfaceVector[1], + 0) + * std::abs(forceMagnitude); + nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0], + interfaceVector[1], + 0) + * std::abs(forceMagnitude); + // 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 { + nodalForces[viPair.first] + = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) + * forceMagnitude / 5; + nodalForces[viPair.second] + = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) + * forceMagnitude / 5; + fixedVertices[viPair.first] = std::unordered_set{2}; + fixedVertices[viPair.second] = std::unordered_set{2}; + } } + scenarios[baseSimulationScenarioIndexOffset_dome + domeSimulationScenarioIndex] + = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_" + + std::to_string(domeSimulationScenarioIndex), + fixedVertices, + nodalForces, + nodalForcedDisplacements)); } - scenarios[simulationScenarioIndex] = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[simulationScenarioIndex], - fixedVertices, - nodalForces, - nodalForcedDisplacements)); //// Saddle - simulationScenarioIndex = BaseSimulationScenario::Saddle; - fixedVertices.clear(); - nodalForces.clear(); - for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); - viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); - viPairIt++) { - const auto &viPair = *viPairIt; - CoordType v = - (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) ^ - CoordType(0, 0, -1).Normalize(); - if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) { - nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.002; - nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.002; - } else { - fixedVertices[viPair.first] = std::unordered_set{2}; - fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; + const double maxForceMagnitude_saddle = 0.005; + const double minForceMagnitude_saddle = -0.005; + const int numberOfSimulationScenarios_saddle + = simulationScenariosResolution[BaseSimulationScenario::Saddle]; + const double forceMagnitudeStep_saddle = numberOfSimulationScenarios_saddle == 1 + ? maxForceMagnitude_saddle + : (maxForceMagnitude_saddle + - minForceMagnitude_saddle) + / (numberOfSimulationScenarios_saddle - 1); + const int baseSimulationScenarioIndexOffset_saddle + = std::accumulate(simulationScenariosResolution.begin(), + simulationScenariosResolution.begin() + BaseSimulationScenario::Saddle, + 0); + for (int saddleSimulationScenarioIndex = 0; + saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle; + saddleSimulationScenarioIndex++) { + fixedVertices.clear(); + nodalForces.clear(); + const double forceMagnitude = (forceMagnitudeStep_saddle * saddleSimulationScenarioIndex + + minForceMagnitude_saddle); + for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); + viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); + viPairIt++) { + const auto &viPair = *viPairIt; + CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) + ^ CoordType(0, 0, -1).Normalize(); + if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) { + nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude; + nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude; + } else { + fixedVertices[viPair.first] = std::unordered_set{2}; + fixedVertices[viPair.second] = std::unordered_set{0, 1, 2}; - nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.001; - nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.001; + nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude + / 2; + nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude + / 2; + } } + scenarios[baseSimulationScenarioIndexOffset_saddle + saddleSimulationScenarioIndex] + = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_" + + std::to_string(saddleSimulationScenarioIndex), + fixedVertices, + nodalForces, + {})); } - scenarios[simulationScenarioIndex] = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[simulationScenarioIndex], - fixedVertices, - nodalForces, - {})); return scenarios; } @@ -1050,31 +1181,32 @@ Results ReducedModelOptimizer::optimize( global.numOfSimulationCrashes = 0; global.numberOfFunctionCalls = 0; global.optimizationSettings = optimizationSettings; - global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh); + global.fullPatternSimulationJobs = createFullPatternSimulationScenarios( + m_pFullPatternSimulationMesh); // polyscope::removeAllStructures(); DRMSimulationModel::Settings simulationSettings; simulationSettings.shouldDraw = false; // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( - // ReducedModelOptimization::Colors::fullInitial); + // ReducedPatternOptimization::Colors::fullInitial); // LinearSimulationModel linearSimulator; for (int simulationScenarioIndex : global.simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = global.fullPatternSimulationJobs[simulationScenarioIndex]; SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, simulationSettings); - // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( - // pFullPatternSimulationJob); - // fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); + // fullPatternResults.registerForDrawing(ReducedPatternOptimization::Colors::fullDeformed, // true, // true); - // fullPatternResults_linear.labelPrefix += "_linear"; - // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, - // true, - // true); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + // true, + // true); // polyscope::show(); // fullPatternResults.unregister(); - // fullPatternResults_linear.unregister(); + // fullPatternResults_linear.unregister(); global.fullPatternResults[simulationScenarioIndex] = fullPatternResults; SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; @@ -1083,6 +1215,7 @@ Results ReducedModelOptimizer::optimize( reducedPatternSimulationJob); global.reducedPatternSimulationJobs[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); + // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; } // global.fullPatternSimulationJobs[0]->pMesh->unregister(); @@ -1091,12 +1224,6 @@ Results ReducedModelOptimizer::optimize( computeObjectiveValueNormalizationFactors(); // } Results optResults = runOptimization(optimizationSettings); - for (int simulationScenarioIndex : global.simulationScenarioIndices) { - optResults.fullPatternSimulationJobs.push_back( - global.fullPatternSimulationJobs[simulationScenarioIndex]); - optResults.reducedPatternSimulationJobs.push_back( - global.reducedPatternSimulationJobs[simulationScenarioIndex]); - } return optResults; } diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index 1da7286..bcac219 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -9,6 +9,8 @@ #include "reducedmodeloptimizer_structs.hpp" #include "simulationmesh.hpp" #include +#include +#include #ifdef POLYSCOPE_DEFINED #include "polyscope/color_management.h" @@ -28,17 +30,17 @@ class ReducedModelOptimizer std::unordered_map> slotToNode; public: - constexpr static std::array simulationScenariosResolution = {1, 1, 1, 1, 1}; + constexpr static std::array simulationScenariosResolution = {2, 2, 6, 6, 6}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.end(), 0); inline static int fanSize{6}; inline static VectorType patternPlaneNormal{0, 0, 1}; - ReducedModelOptimization::Results optimize( - const ReducedModelOptimization::Settings &xRanges, - const std::vector &simulationScenarios - = std::vector()); + ReducedPatternOptimization::Results optimize( + const ReducedPatternOptimization::Settings &xRanges, + const std::vector &simulationScenarios + = std::vector()); double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); @@ -64,7 +66,7 @@ public: double innerHexagonRotationAngle); static double objective(double b, double r, double E); - static std::vector> createScenarios( + static std::vector> createFullPatternSimulationScenarios( const std::shared_ptr &pMesh, const std::unordered_map &fullPatternOppositeInterfaceViMap); @@ -85,7 +87,7 @@ public: static void visualizeResults( const std::vector> &fullPatternSimulationJobs, const std::vector> &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, + const std::vector &simulationScenarios, const std::unordered_map &reducedToFullInterfaceViMap); static void registerResultsForDrawing( @@ -94,7 +96,7 @@ public: const std::unordered_map &reducedToFullInterfaceViMap); - static double computeRawDisplacementError( + static double computeRawTranslationalError( const std::vector &fullPatternDisplacements, const std::vector &reducedPatternDisplacements, const std::unordered_map @@ -136,9 +138,9 @@ private: const SimulationResults &fullModelResults, const std::unordered_map &displacementsReducedToFullMap, Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static ReducedModelOptimization::Results runOptimization( - const ReducedModelOptimization::Settings &settings); - std::vector> createScenarios( + static ReducedPatternOptimization::Results runOptimization( + const ReducedPatternOptimization::Settings &settings); + std::vector> createFullPatternSimulationScenarios( const std::shared_ptr &pMesh); void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); @@ -148,6 +150,9 @@ private: static double objective(long n, const double *x); DRMSimulationModel simulator; void computeObjectiveValueNormalizationFactors(); + static ReducedPatternOptimization::Results getResults( + const dlib::function_evaluation &optimizationResult_dlib, + const ReducedPatternOptimization::Settings &settings); }; void updateMesh(long n, const double *x); #endif // REDUCEDMODELOPTIMIZER_HPP