diff --git a/TestSet/ReducedPatterns/SingleBar_reduced.ply b/TestSet/ReducedPatterns/SingleBar_reduced.ply deleted file mode 100644 index 146888d..0000000 --- a/TestSet/ReducedPatterns/SingleBar_reduced.ply +++ /dev/null @@ -1,25 +0,0 @@ -ply -format ascii 1.0 -comment VCGLIB generated -element vertex 7 -property double x -property double y -property double z -property uchar red -property uchar green -property uchar blue -property uchar alpha -element face 0 -property list uchar int vertex_indices -element edge 1 -property int vertex1 -property int vertex2 -end_header -0 0 0 255 255 255 255 --0.1666666666666666 -0.2886751345948129 0 255 255 255 255 --0.3333333333333333 -0.5773502691896257 0 255 255 255 255 -0 -0.8660254037844387 0 255 255 255 255 -0.3333333333333333 -0.5773502691896258 0 255 255 255 255 -0.1666666666666666 -0.288675134594813 0 255 255 255 255 -0.1 -0.5773502691896258 0 255 255 255 255 -0 3 diff --git a/src/main.cpp b/src/main.cpp index c717503..b83a77a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,5 +1,5 @@ -#include "drmsimulationmodel.hpp" #include "csvfile.hpp" +#include "drmsimulationmodel.hpp" #include "edgemesh.hpp" #include "reducedmodeloptimizer.hpp" #include "simulationhistoryplotter.hpp" @@ -57,6 +57,8 @@ int main(int argc, char *argv[]) { = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon; settings_optimization.normalizationParameter = 0.0003; settings_optimization.solutionAccuracy = 0.001; + settings_optimization.objectiveWeights.translational = std::atof(argv[4]); + settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]); // Optimize pair const std::string pairName = @@ -66,28 +68,37 @@ int main(int argc, char *argv[]) { assert(interfaceNodeIndex==numberOfNodesPerSlot[0]+numberOfNodesPerSlot[3]); ReducedModelOptimizer optimizer(numberOfNodesPerSlot); optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.xRanges.size()); - ReducedPatternOptimization::Results optimizationResults = optimizer.optimize( - settings_optimization); + ReducedPatternOptimization::Results optimizationResults + // = optimizer.optimize( + // settings_optimization); + = optimizer.optimize(settings_optimization); // Export results - const bool input_resultDirectoryDefined = argc >= 5; - std::string optimizationResultsDirectory = - input_resultDirectoryDefined ? argv[4] : std::filesystem::current_path().append("OptimizationResults"); + const bool input_resultDirectoryDefined = argc >= 6; + std::string optimizationResultsDirectory = input_resultDirectoryDefined + ? argv[5] + : std::filesystem::current_path().append( + "OptimizationResults"); std::string resultsOutputDir; if (optimizationResults.numberOfSimulationCrashes != 0) { - const auto crashedJobsDirPath = - std::filesystem::path(optimizationResultsDirectory) - .append("CrashedJobs") - .append(pairName); - std::filesystem::create_directories(crashedJobsDirPath); - resultsOutputDir = crashedJobsDirPath.string(); + const auto crashedJobsDirPath = std::filesystem::path( + optimizationResultsDirectory.append("CrashedJobs") + .append( + pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_" + + to_string_with_precision(settings_optimization.objectiveWeights.translational) + + ")")); + std::filesystem::create_directories(crashedJobsDirPath); + resultsOutputDir = crashedJobsDirPath.string(); } else { - std::filesystem::path convergedJobsDirPath( - std::filesystem::path(optimizationResultsDirectory) - .append("ConvergedJobs") - .append(pairName)); - std::filesystem::create_directories(convergedJobsDirPath); - resultsOutputDir = convergedJobsDirPath.string(); + std::filesystem::path convergedJobsDirPath( + std::filesystem::path(optimizationResultsDirectory) + .append("ConvergedJobs") + .append( + pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_" + + to_string_with_precision(settings_optimization.objectiveWeights.translational) + + ")")); + std::filesystem::create_directories(convergedJobsDirPath); + resultsOutputDir = convergedJobsDirPath.string(); } optimizationResults.save(resultsOutputDir); // Write results in csv diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 4673b45..31294c7 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -16,6 +16,8 @@ struct GlobalOptimizationVariables { std::vector> reducedPatternSimulationJobs; std::unordered_map reducedToFullInterfaceViMap; + std::vector> + fullPatternInterfaceViPairs; matplot::line_handle gPlotHandle; std::vector objectiveValueHistory; Eigen::VectorXd initialParameters; @@ -27,6 +29,14 @@ struct GlobalOptimizationVariables { int numberOfOptimizationParameters{5}; ReducedPatternOptimization::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; } global; double ReducedModelOptimizer::computeDisplacementError( @@ -112,7 +122,8 @@ double ReducedModelOptimizer::computeError( simulationResults_reducedPattern.rotationalDisplacementQuaternion, reducedToFullInterfaceViMap, normalizationFactor_rotationalDisplacement); - return 1.5 * translationalError + 0.5 * rotationalError; + return global.optimizationSettings.objectiveWeights.translational * translationalError + + global.optimizationSettings.objectiveWeights.rotational * rotationalError; } double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3, @@ -124,11 +135,12 @@ double ReducedModelOptimizer::objective(double E,double A,double J,double I2,dou double ReducedModelOptimizer::objective(long n, const double *x) { // std::cout.precision(17); - // for (int i = 0; i < n; i++) { - // std::cout << x[i] << " "; - // } + // for (int i = 0; i < n; i++) { + // std::cout << x[i] << " "; + // } // std::cout << std::endl; + // std::cout << x[n - 2] << " " << x[n - 1] << std::endl; // const Element &e = // global.reducedPatternSimulationJobs[0]->pMesh->elements[0]; std::cout << // e.axialConstFactor << " " << e.torsionConstFactor << " " @@ -136,6 +148,11 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // e.secondBendingConstFactor // << std::endl; updateMesh(n, x); + // global.reducedPatternSimulationJobs[0]->pMesh->registerForDrawing(); + // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing(); + // polyscope::show(); + // global.reducedPatternSimulationJobs[0]->pMesh->unregister(); + // std::cout << e.axialConstFactor << " " << e.torsionConstFactor << " " // << e.firstBendingConstFactor << " " << // e.secondBendingConstFactor @@ -144,10 +161,11 @@ double ReducedModelOptimizer::objective(long n, const double *x) { // run simulations double totalError = 0; LinearSimulationModel simulator; -// FormFinder simulator; + // FormFinder simulator; for (const int simulationScenarioIndex : global.simulationScenarioIndices) { - SimulationResults reducedModelResults = simulator.executeSimulation( - global.reducedPatternSimulationJobs[simulationScenarioIndex]); + const std::shared_ptr &reducedJob + = global.reducedPatternSimulationJobs[simulationScenarioIndex]; + SimulationResults reducedModelResults = simulator.executeSimulation(reducedJob); // std::string filename; if (!reducedModelResults.converged) { totalError += std::numeric_limits::max(); @@ -218,7 +236,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { } #ifdef POLYSCOPE_DEFINED ++global.numberOfFunctionCalls; - if (global.numberOfFunctionCalls % 100 == 0) { + if (global.numberOfFunctionCalls % 1000 == 0) { std::cout << "Number of function calls:" << global.numberOfFunctionCalls << std::endl; } #endif @@ -276,7 +294,7 @@ void ReducedModelOptimizer::computeMaps( std::unordered_map &fullToReducedInterfaceViMap, std::vector> - &fullPatternOppositeInterfaceViMap) + &fullPatternOppositeInterfaceViPairs) { // Compute the offset between the interface nodes const size_t interfaceSlotIndex = 4; // bottom edge @@ -320,7 +338,11 @@ void ReducedModelOptimizer::computeMaps( pu_reducedModel.remap[baseTriangleInterfaceVi]; const size_t reducedModelInterfaceVertexOffset = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/; - reducedPattern.createFan({1}); //TODO: should be an input parameter from main + 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] = @@ -331,20 +353,21 @@ void ReducedModelOptimizer::computeMaps( constructInverseMap(reducedToFullInterfaceViMap, fullToReducedInterfaceViMap); // Create opposite vertex map - fullPatternOppositeInterfaceViMap.clear(); - fullPatternOppositeInterfaceViMap.reserve(fanSize / 2); + 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()); - fullPatternOppositeInterfaceViMap.emplace_back(std::make_pair(vi0, vi1)); + fullPatternOppositeInterfaceViPairs.emplace_back(std::make_pair(vi0, vi1)); } + global.fullPatternInterfaceViPairs = fullPatternOppositeInterfaceViPairs; +#if POLYSCOPE_DEFINED const bool debugMapping = false; if (debugMapping) { -#if POLYSCOPE_DEFINED reducedPattern.registerForDrawing(); // std::vector colors_reducedPatternExcludedEdges( @@ -361,8 +384,7 @@ void ReducedModelOptimizer::computeMaps( std::vector nodeColorsOpposite(fullPattern.VN(), glm::vec3(0, 0, 0)); - for (const std::pair oppositeVerts : - fullPatternOppositeInterfaceViMap) { + for (const std::pair oppositeVerts : fullPatternOppositeInterfaceViPairs) { auto color = polyscope::getNextUniqueColor(); nodeColorsOpposite[oppositeVerts.first] = color; nodeColorsOpposite[oppositeVerts.second] = color; @@ -395,8 +417,8 @@ void ReducedModelOptimizer::computeMaps( nodeColorsReducedToFull_full) ->setEnabled(true); polyscope::show(); -#endif } +#endif } void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, @@ -462,12 +484,14 @@ const double I3=global.initialParameters(4) * x[4]; e.I3 = I3; } assert(pReducedPatternSimulationMesh->EN() == 12); - assert(n>=2); + assert(n >= 2); - CoordType center_barycentric(1, 0, 0); - CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); - CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric) - * x[n - 2]); + // CoordType center_barycentric(1, 0, 0); + // CoordType interfaceEdgeMiddle_barycentric(0, 0.5, 0.5); + // CoordType movableVertex_barycentric((center_barycentric + interfaceEdgeMiddle_barycentric) + // * x[n - 2]); + + CoordType movableVertex_barycentric(1 - x[n - 2], x[n - 2] / 2, x[n - 2] / 2); CoordType baseTriangleMovableVertexPosition = global.baseTriangle.cP(0) * movableVertex_barycentric[0] + global.baseTriangle.cP(1) @@ -480,7 +504,7 @@ const double I3=global.initialParameters(4) * x[4]; for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { - pReducedPatternSimulationMesh->vert[2 * rotationCounter + 1].P() + pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() = vcg::RotationMatrix(ReducedModelOptimizer::patternPlaneNormal, vcg::math::ToRad(60.0 * rotationCounter)) * baseTriangleMovableVertexPosition; @@ -690,26 +714,33 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( assert(global.minX[xVariableIndex] == optimizationResult_dlib.x(xVariableIndex)); optimalX[xVariableIndex] = optimizationResult_dlib.x(xVariableIndex); + std::cout << results.optimalXNameValuePairs[xVariableIndex].first << ":" + << optimalX[xVariableIndex] << " "; } + std::cout << std::endl; // Compute obj value per simulation scenario and the raw objective value updateMesh(optimalX.size(), optimalX.data()); // 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( - totalNumberOfSimulationScenarios); + global.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_rawTranslational.resize( - totalNumberOfSimulationScenarios); - results.objectiveValue.perSimulationScenario_rotational.resize(totalNumberOfSimulationScenarios); + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_rotational.resize( + global.simulationScenarioIndices.size()); results.objectiveValue.perSimulationScenario_rawRotational.resize( - totalNumberOfSimulationScenarios); - results.objectiveValue.perSimulationScenario_total.resize(totalNumberOfSimulationScenarios); - for (int simulationScenarioIndex : global.simulationScenarioIndices) { + global.simulationScenarioIndices.size()); + results.objectiveValue.perSimulationScenario_total.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]); - results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex] = computeError( + results.objectiveValue.perSimulationScenario_total[i] = computeError( global.fullPatternResults[simulationScenarioIndex], reducedModelResults, global.reducedToFullInterfaceViMap, @@ -721,15 +752,13 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( global.fullPatternResults[simulationScenarioIndex].displacements, reducedModelResults.displacements, global.reducedToFullInterfaceViMap); - results.objectiveValue.perSimulationScenario_rawTranslational[simulationScenarioIndex] - = rawTranslationalError; + results.objectiveValue.perSimulationScenario_rawTranslational[i] = rawTranslationalError; //Raw rotational const double rawRotationalError = computeRawRotationalError( global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, reducedModelResults.rotationalDisplacementQuaternion, global.reducedToFullInterfaceViMap); - results.objectiveValue.perSimulationScenario_rawRotational[simulationScenarioIndex] - = rawRotationalError; + results.objectiveValue.perSimulationScenario_rawRotational[i] = rawRotationalError; //Normalized translational const double normalizedTranslationalError = computeDisplacementError( @@ -737,8 +766,7 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( reducedModelResults.displacements, global.reducedToFullInterfaceViMap, global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); - results.objectiveValue.perSimulationScenario_translational[simulationScenarioIndex] - = normalizedTranslationalError; + results.objectiveValue.perSimulationScenario_translational[i] = normalizedTranslationalError; // const double test_normalizedTranslationError = computeDisplacementError( // global.fullPatternResults[simulationScenarioIndex].displacements, // reducedModelResults.displacements, @@ -750,8 +778,7 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( reducedModelResults.rotationalDisplacementQuaternion, global.reducedToFullInterfaceViMap, global.rotationalDisplacementNormalizationValues[simulationScenarioIndex]); - results.objectiveValue.perSimulationScenario_rotational[simulationScenarioIndex] - = normalizedRotationalError; + results.objectiveValue.perSimulationScenario_rotational[i] = normalizedRotationalError; // const double test_normalizedRotationalError = computeRotationalError( // global.fullPatternResults[simulationScenarioIndex].rotationalDisplacementQuaternion, // reducedModelResults.rotationalDisplacementQuaternion, @@ -766,16 +793,15 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( std::cout << "translation normalization value:" << global.translationalDisplacementNormalizationValues[simulationScenarioIndex] << std::endl; - std::cout << "Translational error:" << normalizedTranslationalError << 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 << "Objective value:" - << results.objectiveValue.perSimulationScenario_total[simulationScenarioIndex] + std::cout << "Total Error value:" << results.objectiveValue.perSimulationScenario_total[i] << std::endl; results.objectiveValue.totalRaw += rawTranslationalError + rawRotationalError; std::cout << std::endl; @@ -786,6 +812,9 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::getResults( std::cout << "Finished optimizing." << endl; std::cout << "Total optimal objective value:" << results.objectiveValue.total << std::endl; assert(global.minY == optimizationResult_dlib.y); + if (global.minY != optimizationResult_dlib.y) { + std::cerr << "ERROR in objective value" << std::endl; + } } for (int simulationScenarioIndex : global.simulationScenarioIndices) { @@ -830,19 +859,239 @@ ReducedPatternOptimization::Results ReducedModelOptimizer::runOptimization(const 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 +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}; + } +} + +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 + / 5; + job.nodalExternalForces[viPair.second] + = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) + * forceMagnitude / 5; + 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 fullPatternMaxSimulationForceObjective(const double &forceMagnitude) +{ + const double &desiredDisplacementValue = global.desiredMaxDisplacementValue; + 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.shouldUseTranslationalKineticEnergyThreshold = true; + // settings.totalTranslationalKineticEnergyThreshold = 1e-7; + SimulationResults results = simulator.executeSimulation(std::make_shared(job), + settings); + const double error = std::abs( + // results.displacements[global.fullPatternInterfaceViPairs[1].first].getTranslation().norm() + results.displacements[global.interfaceViForComputingScenarioError].getTranslation().norm() + - desiredDisplacementValue); + + return error; +} + +double ReducedModelOptimizer::getFullPatternMaxSimulationForce(const BaseSimulationScenario &scenario) +{ + double forceMagnitude = 1; + global.desiredMaxDisplacementValue = 0.1 + * vcg::Distance(global.baseTriangle.cP(0), + (global.baseTriangle.cP(1) + + global.baseTriangle.cP(2)) + / 2); + const double optimizationEpsilon = global.desiredMaxDisplacementValue * 0.5; + switch (scenario) { + case Axial: + global.constructScenarioFunction = &ReducedModelOptimizer::constructAxialSimulationScenario; + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; + dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective, + forceMagnitude, + 1e-8, + 1e8, + optimizationEpsilon); + break; + case Shear: + global.constructScenarioFunction = &ReducedModelOptimizer::constructShearSimulationScenario; + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; + dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective, + forceMagnitude, + 1e-8, + 1e8, + optimizationEpsilon); + break; + case Bending: + global.constructScenarioFunction = &ReducedModelOptimizer::constructBendingSimulationScenario; + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; + dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective, + forceMagnitude, + 1e-8, + 1e8, + optimizationEpsilon, + 200); + break; + case Dome: + global.desiredMaxDisplacementValue *= 2; + global.constructScenarioFunction = &ReducedModelOptimizer::constructDomeSimulationScenario; + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[1].first; + dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective, + forceMagnitude, + 1e-8, + 1e8, + optimizationEpsilon, + 200); + break; + case Saddle: + global.desiredMaxDisplacementValue *= 2; + global.constructScenarioFunction = &ReducedModelOptimizer::constructSaddleSimulationScenario; + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; + dlib::find_min_single_variable(&fullPatternMaxSimulationForceObjective, + forceMagnitude, + 1e-8, + 1e8, + optimizationEpsilon, + 150); + break; + } + + return forceMagnitude; +} + +//TODO: Make more compact std::vector> ReducedModelOptimizer::createFullPatternSimulationScenarios( const std::shared_ptr &pMesh) { std::vector> scenarios; scenarios.resize(totalNumberOfSimulationScenarios); - std::unordered_map> fixedVertices; - std::unordered_map nodalForces; + + SimulationJob job; + job.pMesh = pMesh; //// Axial - const double maxForceMagnitude_axial = 500; - const double minForceMagnitude_axial = -500; + const double maxForceMagnitude_axial = getFullPatternMaxSimulationForce( + BaseSimulationScenario::Axial); + const double minForceMagnitude_axial = -maxForceMagnitude_axial; const int numberOfSimulationScenarios_axial = simulationScenariosResolution[BaseSimulationScenario::Axial]; const double forceMagnitudeStep_axial = numberOfSimulationScenarios_axial == 1 @@ -853,38 +1102,27 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.begin() + BaseSimulationScenario::Axial, 0); - for (int axialSimulationScenarioIndex = 0; axialSimulationScenarioIndex < numberOfSimulationScenarios_axial; axialSimulationScenarioIndex++) { + job.nodalExternalForces.clear(); + job.constrainedVertices.clear(); + job.nodalForcedDisplacements.clear(); + job.label = baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_" + + std::to_string(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}; - } - } + constructAxialSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job); scenarios[baseSimulationScenarioIndexOffset_axial + axialSimulationScenarioIndex] - = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[BaseSimulationScenario::Axial] + "_" - + std::to_string(axialSimulationScenarioIndex), - fixedVertices, - nodalForces, - {})); + = std::make_shared(job); } //// Shear - const double maxForceMagnitude_shear = 50; - const double minForceMagnitude_shear = -50; + const double maxForceMagnitude_shear = getFullPatternMaxSimulationForce( + BaseSimulationScenario::Shear); + const double minForceMagnitude_shear = -maxForceMagnitude_shear; const int numberOfSimulationScenarios_shear = simulationScenariosResolution[BaseSimulationScenario::Shear]; const double forceMagnitudeStep_shear = numberOfSimulationScenarios_shear == 1 @@ -898,35 +1136,21 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int shearSimulationScenarioIndex = 0; shearSimulationScenarioIndex < numberOfSimulationScenarios_shear; shearSimulationScenarioIndex++) { - fixedVertices.clear(); - nodalForces.clear(); + job.constrainedVertices.clear(); + job.nodalExternalForces.clear(); + job.label = baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_" + + std::to_string(shearSimulationScenarioIndex); 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}; - } - } + constructShearSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job); scenarios[baseSimulationScenarioIndexOffset_shear + shearSimulationScenarioIndex] - = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[BaseSimulationScenario::Shear] + "_" - + std::to_string(shearSimulationScenarioIndex), - fixedVertices, - nodalForces, - {})); + = std::make_shared(job); } //// Bending - const double maxForceMagnitude_bending = 0.005; - const double minForceMagnitude_bending = -0.005; + const double maxForceMagnitude_bending = getFullPatternMaxSimulationForce( + BaseSimulationScenario::Bending); + const double minForceMagnitude_bending = -maxForceMagnitude_bending; const int numberOfSimulationScenarios_bending = simulationScenariosResolution[BaseSimulationScenario::Bending]; const double forceMagnitudeStep_bending = numberOfSimulationScenarios_bending == 1 @@ -941,27 +1165,23 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int bendingSimulationScenarioIndex = 0; bendingSimulationScenarioIndex < numberOfSimulationScenarios_bending; bendingSimulationScenarioIndex++) { - fixedVertices.clear(); - nodalForces.clear(); + job.nodalExternalForces.clear(); + job.constrainedVertices.clear(); + job.label = baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_" + + std::to_string(bendingSimulationScenarioIndex); 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}; - } + constructBendingSimulationScenario(forceMagnitude, + m_fullPatternOppositeInterfaceViPairs, + job); scenarios[baseSimulationScenarioIndexOffset_bending + bendingSimulationScenarioIndex] - = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[BaseSimulationScenario::Bending] + "_" - + std::to_string(bendingSimulationScenarioIndex), - fixedVertices, - nodalForces, - {})); + = std::make_shared(job); } //// Dome - const double maxForceMagnitude_dome = 0.025; - const double minForceMagnitude_dome = -0.025; + const double maxForceMagnitude_dome = getFullPatternMaxSimulationForce( + BaseSimulationScenario::Dome); + const double minForceMagnitude_dome = -maxForceMagnitude_dome; const int numberOfSimulationScenarios_dome = simulationScenariosResolution[BaseSimulationScenario::Dome]; const double forceMagnitudeStep_dome = numberOfSimulationScenarios_dome == 1 @@ -975,58 +1195,22 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int domeSimulationScenarioIndex = 0; domeSimulationScenarioIndex < numberOfSimulationScenarios_dome; domeSimulationScenarioIndex++) { - fixedVertices.clear(); - nodalForces.clear(); - std::unordered_map nodalForcedDisplacements; + job.constrainedVertices.clear(); + job.nodalExternalForces.clear(); + job.nodalForcedDisplacements.clear(); + job.label = baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_" + + std::to_string(domeSimulationScenarioIndex); 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}; - } - } + constructDomeSimulationScenario(forceMagnitude, m_fullPatternOppositeInterfaceViPairs, job); scenarios[baseSimulationScenarioIndexOffset_dome + domeSimulationScenarioIndex] - = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[BaseSimulationScenario::Dome] + "_" - + std::to_string(domeSimulationScenarioIndex), - fixedVertices, - nodalForces, - nodalForcedDisplacements)); + = std::make_shared(job); } //// Saddle - const double maxForceMagnitude_saddle = 0.005; - const double minForceMagnitude_saddle = -0.005; + const double maxForceMagnitude_saddle = getFullPatternMaxSimulationForce( + BaseSimulationScenario::Saddle); + const double minForceMagnitude_saddle = -maxForceMagnitude_saddle; const int numberOfSimulationScenarios_saddle = simulationScenariosResolution[BaseSimulationScenario::Saddle]; const double forceMagnitudeStep_saddle = numberOfSimulationScenarios_saddle == 1 @@ -1041,37 +1225,18 @@ ReducedModelOptimizer::createFullPatternSimulationScenarios( for (int saddleSimulationScenarioIndex = 0; saddleSimulationScenarioIndex < numberOfSimulationScenarios_saddle; saddleSimulationScenarioIndex++) { - fixedVertices.clear(); - nodalForces.clear(); + job.constrainedVertices.clear(); + job.nodalExternalForces.clear(); + job.nodalForcedDisplacements.clear(); + job.label = baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_" + + std::to_string(saddleSimulationScenarioIndex); 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}) * forceMagnitude - / 2; - nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude - / 2; - } - } + constructSaddleSimulationScenario(forceMagnitude, + m_fullPatternOppositeInterfaceViPairs, + job); scenarios[baseSimulationScenarioIndexOffset_saddle + saddleSimulationScenarioIndex] - = std::make_shared( - SimulationJob(pMesh, - baseSimulationScenarioNames[BaseSimulationScenario::Saddle] + "_" - + std::to_string(saddleSimulationScenarioIndex), - fixedVertices, - nodalForces, - {})); + = std::make_shared(job); } return scenarios; @@ -1181,32 +1346,38 @@ Results ReducedModelOptimizer::optimize( global.numOfSimulationCrashes = 0; global.numberOfFunctionCalls = 0; global.optimizationSettings = optimizationSettings; + global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh; global.fullPatternSimulationJobs = createFullPatternSimulationScenarios( m_pFullPatternSimulationMesh); // polyscope::removeAllStructures(); DRMSimulationModel::Settings simulationSettings; simulationSettings.shouldDraw = false; - // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( - // ReducedPatternOptimization::Colors::fullInitial); + const bool drawFullPatternSimulationResults = false; + if (drawFullPatternSimulationResults) { + global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( + 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(ReducedPatternOptimization::Colors::fullDeformed, - // true, - // true); - // fullPatternResults_linear.labelPrefix += "_linear"; - // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, - // true, - // true); - // polyscope::show(); - // fullPatternResults.unregister(); - // fullPatternResults_linear.unregister(); + if (drawFullPatternSimulationResults) { + // 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); + polyscope::show(); + fullPatternResults.unregister(); + // fullPatternResults_linear.unregister(); + } global.fullPatternResults[simulationScenarioIndex] = fullPatternResults; SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; @@ -1217,12 +1388,14 @@ Results ReducedModelOptimizer::optimize( = std::make_shared(reducedPatternSimulationJob); // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; } - // global.fullPatternSimulationJobs[0]->pMesh->unregister(); - + if (drawFullPatternSimulationResults) { + global.fullPatternSimulationJobs[0]->pMesh->unregister(); + } // if (global.optimizationSettings.normalizationStrategy // != Settings::NormalizationStrategy::NonNormalized) { computeObjectiveValueNormalizationFactors(); // } + PolyscopeInterface::deinitPolyscope(); Results optResults = runOptimization(optimizationSettings); return optResults; diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index bcac219..9b1f60d 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -30,12 +30,13 @@ class ReducedModelOptimizer std::unordered_map> slotToNode; public: - constexpr static std::array simulationScenariosResolution = {2, 2, 6, 6, 6}; + constexpr static std::array simulationScenariosResolution = {4, 4, 4, 4, 4}; inline static int totalNumberOfSimulationScenarios = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.end(), 0); inline static int fanSize{6}; + inline static double initialHexagonSize{0.3}; inline static VectorType patternPlaneNormal{0, 0, 1}; ReducedPatternOptimization::Results optimize( const ReducedPatternOptimization::Settings &xRanges, @@ -132,6 +133,35 @@ public: &reducedToFullInterfaceViMap, const double &normalizationFactor_translationalDisplacement, const double &normalizationFactor_rotationalDisplacement); + 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); private: static void computeDesiredReducedModelDisplacements( @@ -153,6 +183,9 @@ private: static ReducedPatternOptimization::Results getResults( const dlib::function_evaluation &optimizationResult_dlib, const ReducedPatternOptimization::Settings &settings); + std::vector getFullPatternMaxSimulationForces(); + double getFullPatternMaxSimulationForce( + const ReducedPatternOptimization::BaseSimulationScenario &scenario); }; void updateMesh(long n, const double *x); #endif // REDUCEDMODELOPTIMIZER_HPP