diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ddb879..2f94c9a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,8 @@ cmake_minimum_required(VERSION 2.8) project(ReducedModelOptimization) +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + message(STATUS "The compiler ${CMAKE_CXX_COMPILER}") #Add the project cmake scripts to the module path list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) diff --git a/src/main.cpp b/src/main.cpp index b20fca1..5f10e44 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -66,8 +66,7 @@ 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); + ReducedModelOptimization::Results optimizationResults = optimizer.optimize(settings_optimization); // Export results const bool input_resultDirectoryDefined = argc >= 5; diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 48114e4..ede5cae 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -10,7 +10,6 @@ using namespace ReducedModelOptimization; struct GlobalOptimizationVariables { - std::vector g_optimalReducedModelDisplacements; // std::vector> fullPatternDisplacements; std::vector fullPatternResults; std::vector translationalDisplacementNormalizationValues; @@ -22,8 +21,7 @@ struct GlobalOptimizationVariables { matplot::line_handle gPlotHandle; std::vector objectiveValueHistory; Eigen::VectorXd initialParameters; - std::vector - simulationScenarioIndices; + std::vector simulationScenarioIndices; std::vector g_innerHexagonVectors{6, VectorType(0, 0, 0)}; double innerHexagonInitialRotationAngle{30}; double innerHexagonInitialPos{0}; @@ -36,8 +34,6 @@ struct GlobalOptimizationVariables { vcg::Triangle3 baseTriangle; } global; -std::vector reducedPatternMaximumDisplacementSimulationJobs; - double ReducedModelOptimizer::computeDisplacementError( const std::vector &fullPatternDisplacements, const std::vector &reducedPatternDisplacements, @@ -48,11 +44,7 @@ double ReducedModelOptimizer::computeDisplacementError( const double rawError = computeRawDisplacementError(fullPatternDisplacements, reducedPatternDisplacements, reducedToFullInterfaceViMap); - if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { - return rawError / normalizationFactor; - } - return rawError; + return rawError / normalizationFactor; } double ReducedModelOptimizer::computeRawDisplacementError( @@ -75,7 +67,6 @@ double ReducedModelOptimizer::computeRawDisplacementError( fullPatternDisplacements[fullModelVi][2]); Eigen::Vector3d errorVector = fullVertexDisplacement - reducedVertexDisplacement; - // error += errorVector.squaredNorm(); error += errorVector.norm(); } @@ -181,7 +172,7 @@ double ReducedModelOptimizer::objective(long n, const double *x) { reducedModelResults.internalPotentialEnergy - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); } else { - simulationScenarioError += computeError( + simulationScenarioError = computeError( global.fullPatternResults[simulationScenarioIndex], reducedModelResults, global.reducedToFullInterfaceViMap, @@ -287,14 +278,15 @@ void ReducedModelOptimizer::createSimulationMeshes( void ReducedModelOptimizer::computeMaps( const std::unordered_map> &slotToNode, - PatternGeometry &fullPattern, PatternGeometry &reducedPattern, + PatternGeometry &fullPattern, + PatternGeometry &reducedPattern, std::unordered_map &reducedToFullInterfaceViMap, std::unordered_map &fullToReducedInterfaceViMap, - std::unordered_map - &fullPatternOppositeInterfaceViMap) { - + std::vector> + &fullPatternOppositeInterfaceViMap) +{ // Compute the offset between the interface nodes const size_t interfaceSlotIndex = 4; // bottom edge assert(slotToNode.find(interfaceSlotIndex) != slotToNode.end() && @@ -349,12 +341,14 @@ void ReducedModelOptimizer::computeMaps( // Create opposite vertex map fullPatternOppositeInterfaceViMap.clear(); - for (int fanIndex = fanSize / 2 - 1; fanIndex >= 0; fanIndex--) { - const size_t vi0 = fullModelBaseTriangleInterfaceVi + - fanIndex * fullPatternInterfaceVertexOffset; + fullPatternOppositeInterfaceViMap.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[vi0] = vi1; + fullPatternOppositeInterfaceViMap.emplace_back(std::make_pair(vi0, vi1)); } const bool debugMapping = false; @@ -422,11 +416,11 @@ void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, reducedPattern, global.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap, - m_fullPatternOppositeInterfaceViMap); + m_fullPatternOppositeInterfaceViPairs); } -ReducedModelOptimizer::ReducedModelOptimizer( - const std::vector &numberOfNodesPerSlot) { +ReducedModelOptimizer::ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot) +{ FlatPatternTopology::constructNodeToSlotMap(numberOfNodesPerSlot, nodeToSlot); FlatPatternTopology::constructSlotToNodeMap(nodeToSlot, slotToNode); } @@ -561,103 +555,102 @@ void ReducedModelOptimizer::computeReducedModelSimulationJob( simulationJobOfReducedModel.nodalForcedDisplacements = reducedNodalForcedDisplacements; } -#if POLYSCOPE_DEFINED -void ReducedModelOptimizer::visualizeResults( - const std::vector> - &fullPatternSimulationJobs, - const std::vector> - &reducedPatternSimulationJobs, - const std::vector &simulationScenarios, - const std::unordered_map - &reducedToFullInterfaceViMap) { - DRMSimulationModel simulator; - std::shared_ptr pFullPatternSimulationMesh = - fullPatternSimulationJobs[0]->pMesh; - pFullPatternSimulationMesh->registerForDrawing(); - pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel() - + "_undeformed.ply"); - double totalError = 0; - for (const int simulationScenarioIndex : simulationScenarios) { - const std::shared_ptr &pFullPatternSimulationJob = - fullPatternSimulationJobs[simulationScenarioIndex]; - pFullPatternSimulationJob->registerForDrawing( - pFullPatternSimulationMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // fullModelResults.saveDeformedModel(); - const std::shared_ptr &pReducedPatternSimulationJob = - reducedPatternSimulationJobs[simulationScenarioIndex]; - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - double normalizationFactor = 1; - if (global.optimizationSettings.normalizationStrategy != - ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) { - normalizationFactor - = global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; - } - reducedModelResults.saveDeformedModel(); - fullModelResults.saveDeformedModel(); - double error = computeDisplacementError(reducedModelResults.displacements, - fullModelResults.displacements, - reducedToFullInterfaceViMap, - normalizationFactor); - std::cout << "Error of simulation scenario " - << ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex] << " is " - << error << std::endl; - totalError += error; - reducedModelResults.registerForDrawing(); - // firstOptimizationRoundResults[simulationScenarioIndex].registerForDrawing(); - // registerWorldAxes(); - const std::string screenshotFilename = - "/home/iason/Coding/Projects/Approximating shapes with flat " - "patterns/RodModelOptimizationForPatterns/build/OptimizationResults/" - "Images/" + - pFullPatternSimulationMesh->getLabel() + "_" + - ReducedModelOptimization::simulationScenarioStrings[simulationScenarioIndex]; - polyscope::show(); - polyscope::screenshot(screenshotFilename, false); - fullModelResults.unregister(); - reducedModelResults.unregister(); - // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); - } - std::cout << "Total error:" << totalError << std::endl; -} +//#if POLYSCOPE_DEFINED +//void ReducedModelOptimizer::visualizeResults( +// const std::vector> &fullPatternSimulationJobs, +// const std::vector> &reducedPatternSimulationJobs, +// const std::vector &simulationScenarios, +// const std::unordered_map +// &reducedToFullInterfaceViMap) +//{ +// DRMSimulationModel simulator; +// std::shared_ptr pFullPatternSimulationMesh = +// fullPatternSimulationJobs[0]->pMesh; +// pFullPatternSimulationMesh->registerForDrawing(); +// pFullPatternSimulationMesh->save(pFullPatternSimulationMesh->getLabel() + "_undeformed.ply"); +// reducedPatternSimulationJobs[0]->pMesh->save(reducedPatternSimulationJobs[0]->pMesh->getLabel() +// + "_undeformed.ply"); +// double totalError = 0; +// for (const int simulationScenarioIndex : simulationScenarios) { +// const std::shared_ptr &pFullPatternSimulationJob = +// fullPatternSimulationJobs[simulationScenarioIndex]; +// pFullPatternSimulationJob->registerForDrawing( +// pFullPatternSimulationMesh->getLabel()); +// SimulationResults fullModelResults = +// simulator.executeSimulation(pFullPatternSimulationJob); +// fullModelResults.registerForDrawing(); +// // fullModelResults.saveDeformedModel(); +// const std::shared_ptr &pReducedPatternSimulationJob = +// reducedPatternSimulationJobs[simulationScenarioIndex]; +// SimulationResults reducedModelResults = +// simulator.executeSimulation(pReducedPatternSimulationJob); +// double normalizationFactor = 1; +// if (global.optimizationSettings.normalizationStrategy != +// ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) { +// normalizationFactor +// = global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; +// } +// reducedModelResults.saveDeformedModel(); +// fullModelResults.saveDeformedModel(); +// double error = computeDisplacementError(reducedModelResults.displacements, +// fullModelResults.displacements, +// reducedToFullInterfaceViMap, +// normalizationFactor); +// std::cout << "Error of simulation scenario " +// getLabel() + "_" +// + baseSimulationScenarioNames[simulationScenarioIndex]; +// polyscope::show(); +// polyscope::screenshot(screenshotFilename, false); +// fullModelResults.unregister(); +// reducedModelResults.unregister(); +// // firstOptimizationRoundResults[simulationScenarioIndex].unregister(); +// } +// std::cout << "Total error:" << totalError << std::endl; +//} -void ReducedModelOptimizer::registerResultsForDrawing( - const std::shared_ptr &pFullPatternSimulationJob, - const std::shared_ptr &pReducedPatternSimulationJob, - const std::unordered_map - &reducedToFullInterfaceViMap) { - DRMSimulationModel simulator; - std::shared_ptr pFullPatternSimulationMesh = - pFullPatternSimulationJob->pMesh; - pFullPatternSimulationMesh->registerForDrawing(); - // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() - // + - // "_undeformed.ply"); - // reducedPatternSimulationJobs[0]->pMesh->savePly( - // reducedPatternSimulationJobs[0]->pMesh->getLabel() + - // "_undeformed.ply"); - pFullPatternSimulationJob->registerForDrawing( - pFullPatternSimulationMesh->getLabel()); - SimulationResults fullModelResults = - simulator.executeSimulation(pFullPatternSimulationJob); - fullModelResults.registerForDrawing(); - // fullModelResults.saveDeformedModel(); - SimulationResults reducedModelResults = - simulator.executeSimulation(pReducedPatternSimulationJob); - // reducedModelResults.saveDeformedModel(); - // fullModelResults.saveDeformedModel(); - double error = computeRawDisplacementError( - reducedModelResults.displacements, fullModelResults.displacements, - reducedToFullInterfaceViMap/*, - global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); - std::cout << "Error is " << error << std::endl; - reducedModelResults.registerForDrawing(); -} -#endif // POLYSCOPE_DEFINED +//void ReducedModelOptimizer::registerResultsForDrawing( +// const std::shared_ptr &pFullPatternSimulationJob, +// const std::shared_ptr &pReducedPatternSimulationJob, +// const std::unordered_map +// &reducedToFullInterfaceViMap) { +// DRMSimulationModel simulator; +// std::shared_ptr pFullPatternSimulationMesh = +// pFullPatternSimulationJob->pMesh; +// pFullPatternSimulationMesh->registerForDrawing(); +// // pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() +// // + +// // "_undeformed.ply"); +// // reducedPatternSimulationJobs[0]->pMesh->savePly( +// // reducedPatternSimulationJobs[0]->pMesh->getLabel() + +// // "_undeformed.ply"); +// pFullPatternSimulationJob->registerForDrawing( +// pFullPatternSimulationMesh->getLabel()); +// SimulationResults fullModelResults = +// simulator.executeSimulation(pFullPatternSimulationJob); +// fullModelResults.registerForDrawing(); +// // fullModelResults.saveDeformedModel(); +// SimulationResults reducedModelResults = +// simulator.executeSimulation(pReducedPatternSimulationJob); +// // reducedModelResults.saveDeformedModel(); +// // fullModelResults.saveDeformedModel(); +// double error = computeRawDisplacementError( +// reducedModelResults.displacements, fullModelResults.displacements, +// reducedToFullInterfaceViMap/*, +// global.reducedPatternMaximumDisplacementNormSum[simulationScenarioIndex]*/); +// std::cout << "Error is " << error << std::endl; +// reducedModelResults.registerForDrawing(); +//} +//#endif // POLYSCOPE_DEFINED void ReducedModelOptimizer::computeDesiredReducedModelDisplacements( const SimulationResults &fullModelResults, @@ -709,31 +702,31 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { results.optimalXNameValuePairs[xVariableIndex] = std::make_pair(settings.xRanges[xVariableIndex].label, global.minX[xVariableIndex]); } - results.objectiveValue = global.minY; -#ifdef POLYSCOPE_DEFINED - std::cout<<"Total optimal objective value:"< optimalX(results.optimalXNameValuePairs.size()); - for(int xVariableIndex=0;xVariableIndex> ReducedModelOptimizer::createScenarios( const std::shared_ptr &pMesh) { std::vector> scenarios; - scenarios.resize(SimulationScenario::NumberOfSimulationScenarios); + scenarios.resize(totalNumberOfSimulationScenarios); std::unordered_map> fixedVertices; std::unordered_map nodalForces; - const double forceMagnitude = 10; //// Axial - SimulationScenario scenarioName = SimulationScenario::Axial; - // NewMethod - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { - if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { + 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}) * - forceMagnitude * 10; + 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}; } } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + scenarios[simulationScenarioIndex] = std::make_shared( + SimulationJob(pMesh, + "Axial_" + std::to_string(simulationScenarioIndex), + fixedVertices, + nodalForces, + {})); //// Shear - scenarioName = SimulationScenario::Shear; + simulationScenarioIndex = BaseSimulationScenario::Shear; fixedVertices.clear(); nodalForces.clear(); // NewMethod - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { - if (viPairIt != m_fullPatternOppositeInterfaceViMap.begin()) { + 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 * 4; + 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}; } } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + scenarios[simulationScenarioIndex] = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[simulationScenarioIndex], + fixedVertices, + nodalForces, + {})); //// Bending - scenarioName = SimulationScenario::Bending; + simulationScenarioIndex = BaseSimulationScenario::Bending; fixedVertices.clear(); nodalForces.clear(); - for (const auto &viPair : m_fullPatternOppositeInterfaceViMap) { - nodalForces[viPair.first] = Vector6d({0, 0, forceMagnitude, 0, 0, 0}) * 0.0005; + 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}; } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + scenarios[simulationScenarioIndex] = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[simulationScenarioIndex], + fixedVertices, + nodalForces, + {})); //// Dome - scenarioName = SimulationScenario::Dome; + simulationScenarioIndex = BaseSimulationScenario::Dome; fixedVertices.clear(); nodalForces.clear(); std::unordered_map nodalForcedDisplacements; - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + for (auto viPairIt = m_fullPatternOppositeInterfaceViPairs.begin(); + viPairIt != m_fullPatternOppositeInterfaceViPairs.end(); + viPairIt++) { const auto viPair = *viPairIt; - if (viPairIt == m_fullPatternOppositeInterfaceViMap.begin()) { + if (viPairIt == m_fullPatternOppositeInterfaceViPairs.begin()) { CoordType interfaceV = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()); nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceV[0], @@ -865,47 +903,45 @@ ReducedModelOptimizer::createScenarios( 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}) * forceMagnitude - * 0.0005; - nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude - * 0.0005; + 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; } } - scenarios[scenarioName] = std::make_shared( + scenarios[simulationScenarioIndex] = std::make_shared( SimulationJob(pMesh, - simulationScenarioStrings[scenarioName], + baseSimulationScenarioNames[simulationScenarioIndex], fixedVertices, nodalForces, nodalForcedDisplacements)); //// Saddle - scenarioName = SimulationScenario::Saddle; + simulationScenarioIndex = BaseSimulationScenario::Saddle; fixedVertices.clear(); nodalForces.clear(); - for (auto viPairIt = m_fullPatternOppositeInterfaceViMap.begin(); - viPairIt != m_fullPatternOppositeInterfaceViMap.end(); viPairIt++) { + 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_fullPatternOppositeInterfaceViMap.begin()) { - nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.0002 - * forceMagnitude; - nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.0002 - * forceMagnitude; + 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}; - nodalForces[viPair.first] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * 0.0001 - * forceMagnitude; - nodalForces[viPair.second] = Vector6d({0, 0, 0, v[0], v[1], 0}) * 0.0001 - * forceMagnitude; + 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; } } - scenarios[scenarioName] = std::make_shared( - SimulationJob(pMesh, simulationScenarioStrings[scenarioName], - fixedVertices, nodalForces, {})); + scenarios[simulationScenarioIndex] = std::make_shared( + SimulationJob(pMesh, + baseSimulationScenarioNames[simulationScenarioIndex], + fixedVertices, + nodalForces, + {})); return scenarios; } @@ -913,11 +949,11 @@ ReducedModelOptimizer::createScenarios( void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { // Compute the sum of the displacement norms - std::vector fullPatternTranslationalDisplacementNormSum(NumberOfSimulationScenarios); - std::vector fullPatternAngularDistance(NumberOfSimulationScenarios); + std::vector fullPatternTranslationalDisplacementNormSum( + totalNumberOfSimulationScenarios); + std::vector fullPatternAngularDistance(totalNumberOfSimulationScenarios); for (int simulationScenarioIndex : global.simulationScenarioIndices) { double translationalDisplacementNormSum = 0; - double angularDistanceSum = 0; for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { const int fullPatternVi = interfaceViPair.second; //If the full pattern vertex is translationally constrained dont take it into account @@ -936,6 +972,10 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { .fullPatternResults[simulationScenarioIndex] .displacements[fullPatternVi]; translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); + } + double angularDistanceSum = 0; + for (auto interfaceViPair : global.reducedToFullInterfaceViMap) { + const int fullPatternVi = interfaceViPair.second; //If the full pattern vertex is rotationally constrained dont take it into account if (global.fullPatternSimulationJobs[simulationScenarioIndex] ->constrainedVertices.contains(fullPatternVi)) { @@ -951,6 +991,7 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { .rotationalDisplacementQuaternion[fullPatternVi] .angularDistance(Eigen::Quaterniond::Identity()); } + fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] = translationalDisplacementNormSum; fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; @@ -964,72 +1005,91 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() { global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], epsilon_translationalDisplacement); - const double epsilon_rotationalDisplacement = vcg::math::ToRad(2.0); + // const double epsilon_rotationalDisplacement = vcg::math::ToRad(10.0); global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] - = std::max(fullPatternAngularDistance[simulationScenarioIndex], - epsilon_rotationalDisplacement); + = /*std::max(*/ fullPatternAngularDistance[simulationScenarioIndex] /*, + epsilon_rotationalDisplacement)*/ + ; } else { global.translationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] = 1; } } - } +} Results ReducedModelOptimizer::optimize( const Settings &optimizationSettings, - const std::vector &simulationScenarios) { - - global.simulationScenarioIndices = simulationScenarios; - if (global.simulationScenarioIndices.empty()) { - global.simulationScenarioIndices = {SimulationScenario::Axial, - SimulationScenario::Shear, - SimulationScenario::Bending, - SimulationScenario::Dome, - SimulationScenario::Saddle}; + const std::vector &desiredBaseSimulationScenarioIndices) +{ + 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)); + } + if (desiredBaseSimulationScenarioIndices.empty()) { + global.simulationScenarioIndices.resize(totalNumberOfSimulationScenarios); + std::iota(global.simulationScenarioIndices.begin(), + global.simulationScenarioIndices.end(), + 0); } - global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios); - global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios); - global.fullPatternResults.resize(NumberOfSimulationScenarios); - global.translationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios); - global.rotationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios); + 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.fullPatternSimulationJobs = - createScenarios(m_pFullPatternSimulationMesh); - reducedPatternMaximumDisplacementSimulationJobs.resize( - NumberOfSimulationScenarios); + global.fullPatternSimulationJobs = createScenarios(m_pFullPatternSimulationMesh); // polyscope::removeAllStructures(); DRMSimulationModel::Settings simulationSettings; simulationSettings.shouldDraw = false; // global.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( // ReducedModelOptimization::Colors::fullInitial); + // LinearSimulationModel linearSimulator; for (int simulationScenarioIndex : global.simulationScenarioIndices) { - const std::shared_ptr &pFullPatternSimulationJob = - global.fullPatternSimulationJobs[simulationScenarioIndex]; + const std::shared_ptr &pFullPatternSimulationJob + = global.fullPatternSimulationJobs[simulationScenarioIndex]; SimulationResults fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, simulationSettings); - // fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed); + // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); + // fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + // true, + // true); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, + // true, + // true); // polyscope::show(); - // fullPatternResults.unregister(); + // fullPatternResults.unregister(); + // fullPatternResults_linear.unregister(); global.fullPatternResults[simulationScenarioIndex] = fullPatternResults; SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, m_fullToReducedInterfaceViMap, reducedPatternSimulationJob); - global.reducedPatternSimulationJobs[simulationScenarioIndex] = - std::make_shared(reducedPatternSimulationJob); + global.reducedPatternSimulationJobs[simulationScenarioIndex] + = std::make_shared(reducedPatternSimulationJob); } - // global.fullPatternSimulationJobs[0]->pMesh->unregister(); + // global.fullPatternSimulationJobs[0]->pMesh->unregister(); - if (global.optimizationSettings.normalizationStrategy != - Settings::NormalizationStrategy::NonNormalized) { - computeObjectiveValueNormalizationFactors(); - } + // if (global.optimizationSettings.normalizationStrategy + // != Settings::NormalizationStrategy::NonNormalized) { + computeObjectiveValueNormalizationFactors(); + // } Results optResults = runOptimization(optimizationSettings); for (int simulationScenarioIndex : global.simulationScenarioIndices) { optResults.fullPatternSimulationJobs.push_back( diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index ca4561d..1da7286 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -22,127 +22,132 @@ class ReducedModelOptimizer std::shared_ptr m_pFullPatternSimulationMesh; std::unordered_map m_fullToReducedInterfaceViMap; - std::unordered_map - m_fullPatternOppositeInterfaceViMap; + std::vector> + m_fullPatternOppositeInterfaceViPairs; std::unordered_map nodeToSlot; std::unordered_map> slotToNode; -#ifdef POLYSCOPE_DEFINED -#endif // POLYSCOPE_DEFINED public: - 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()); - double operator()(const Eigen::VectorXd &x, Eigen::VectorXd &) const; + constexpr static std::array simulationScenariosResolution = {1, 1, 1, 1, 1}; + 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()); + 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); + ReducedModelOptimizer(const std::vector &numberOfNodesPerSlot); + static void computeReducedModelSimulationJob( + const SimulationJob &simulationJobOfFullModel, + const std::unordered_map &fullToReducedMap, + SimulationJob &simulationJobOfReducedModel); - SimulationJob - getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); + SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob); - void initializePatterns(PatternGeometry &fullPattern, - PatternGeometry &reducedPatterm, - const int &optimizationParameters); + void initializePatterns(PatternGeometry &fullPattern, + PatternGeometry &reducedPatterm, + const int &optimizationParameters); - static void runSimulation(const std::string &filename, - std::vector &x); + static void runSimulation(const std::string &filename, std::vector &x); - static double objective(double x2, double A, double J, double I2, double I3, double x3, - double innerHexagonRotationAngle); - static double objective(double b, double r, double E); + static double objective(double x2, + double A, + double J, + double I2, + double I3, + double x3, + double innerHexagonRotationAngle); + static double objective(double b, double r, double E); - static std::vector> - createScenarios(const std::shared_ptr &pMesh, - const std::unordered_map - &fullPatternOppositeInterfaceViMap); + static std::vector> createScenarios( + 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::unordered_map - &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 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 computeRawDisplacementError( - const std::vector &fullPatternDisplacements, - const std::vector &reducedPatternDisplacements, - const std::unordered_map - &reducedToFullInterfaceViMap); + static double computeRawDisplacementError( + 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 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 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); + 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); - private: - static void computeDesiredReducedModelDisplacements( - const SimulationResults &fullModelResults, - const std::unordered_map &displacementsReducedToFullMap, - Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); - static ReducedModelOptimization::Results runOptimization( - const ReducedModelOptimization::Settings &settings); - std::vector> createScenarios( - const std::shared_ptr &pMesh); - void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); - void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); - static void initializeOptimizationParameters(const std::shared_ptr &mesh, - const int &optimizationParamters); +private: + static void computeDesiredReducedModelDisplacements( + const SimulationResults &fullModelResults, + const std::unordered_map &displacementsReducedToFullMap, + Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); + static ReducedModelOptimization::Results runOptimization( + const ReducedModelOptimization::Settings &settings); + std::vector> createScenarios( + const std::shared_ptr &pMesh); + void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); + void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); + static void initializeOptimizationParameters(const std::shared_ptr &mesh, + const int &optimizationParamters); - static double objective(long n, const double *x); - DRMSimulationModel simulator; - void computeObjectiveValueNormalizationFactors(); + static double objective(long n, const double *x); + DRMSimulationModel simulator; + void computeObjectiveValueNormalizationFactors(); }; void updateMesh(long n, const double *x); #endif // REDUCEDMODELOPTIMIZER_HPP