diff --git a/reducedmodeloptimizer.cpp b/reducedmodeloptimizer.cpp index 7a3ac04..252c4c0 100644 --- a/reducedmodeloptimizer.cpp +++ b/reducedmodeloptimizer.cpp @@ -4,14 +4,17 @@ #include "simulationhistoryplotter.hpp" #include "trianglepatterngeometry.hpp" #include "trianglepattterntopology.hpp" -#include +#ifdef USE_ENSMALLEN +#include "ensmallen.hpp" +#endif #include #include #include #include #include -#include +//#include //#define USE_SCENARIO_WEIGHTS +#include using namespace ReducedModelOptimization; @@ -233,6 +236,7 @@ double ReducedModelOptimizer::objective(const std::vector &x, // for (const int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { const std::shared_ptr &reducedJob = optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]; + //#ifdef POLYSCOPE_DEFINED // std::cout << reducedJob->getLabel() << ":" << std::endl; //#endif @@ -243,7 +247,19 @@ double ReducedModelOptimizer::objective(const std::vector &x, = std::numeric_limits::max(); optimizationState.numOfSimulationCrashes++; #ifdef POLYSCOPE_DEFINED - std::cout << "Failed simulation" << std::endl; + std::cout << "Failed simulation with x:" << std::endl; + std::cout.precision(17); + for (int i = 0; i < x.size(); i++) { + std::cout << x[i] << " "; + } + std::cout << std::endl; + // pReducedPatternSimulationMesh->registerForDrawing(); + // polyscope::show(); + // pReducedPatternSimulationMesh->unregister(); + // std::filesystem::create_directories("failedJob"); + // reducedJob->save("failedJob"); + // Results debugRes; + #endif } else { // const double simulationScenarioError_PE = std::abs( @@ -380,7 +396,7 @@ void ReducedModelOptimizer::createSimulationMeshes(PatternGeometry &fullModel, void ReducedModelOptimizer::computeMaps( const std::unordered_map> &slotToNode, PatternGeometry &fullPattern, - PatternGeometry &reducedPattern, + ReducedModel &reducedPattern, std::unordered_map &reducedToFullInterfaceViMap, std::unordered_map @@ -427,11 +443,11 @@ void ReducedModelOptimizer::computeMaps( const size_t reducedModelBaseTriangleInterfaceVi = pu_reducedModel.remap[baseTriangleInterfaceVi]; const size_t reducedModelInterfaceVertexOffset = reducedPattern.VN() /*- 1*/ /*- reducedModelBaseTriangleInterfaceVi*/; - Results::applyOptimizationResults_innerHexagon(initialHexagonSize, - 0, - baseTriangle, - reducedPattern); - reducedPattern.createFan({0}); //TODO: should be an input parameter from main + Results::applyOptimizationResults_reducedModel_nonFanned(initialHexagonSize, + 0, + baseTriangle, + reducedPattern); + reducedPattern.createFan(); //TODO: should be an input parameter from main for (size_t fanIndex = 0; fanIndex < fanSize; fanIndex++) { reducedToFullInterfaceViMap[reducedModelInterfaceVertexOffset * fanIndex + reducedModelBaseTriangleInterfaceVi] @@ -503,12 +519,11 @@ void ReducedModelOptimizer::computeMaps( #endif } -void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, - PatternGeometry &reducedPattern) +void ReducedModelOptimizer::computeMaps(PatternGeometry &fullPattern, ReducedModel &reducedModel) { ReducedModelOptimizer::computeMaps(slotToNode, fullPattern, - reducedPattern, + reducedModel, optimizationState.reducedToFullInterfaceViMap, m_fullToReducedInterfaceViMap, m_fullPatternOppositeInterfaceViPairs); @@ -529,38 +544,42 @@ ReducedModelOptimizer::ReducedModelOptimizer() scenarioIsSymmetrical[BaseSimulationScenario::Shear] = false; constructBaseScenarioFunctions[BaseSimulationScenario::Bending] = &constructBendingSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::Bending] = false; + scenarioIsSymmetrical[BaseSimulationScenario::Bending] = true; constructBaseScenarioFunctions[BaseSimulationScenario::Dome] = &constructDomeSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::Dome] = false; + scenarioIsSymmetrical[BaseSimulationScenario::Dome] = true; constructBaseScenarioFunctions[BaseSimulationScenario::Saddle] = &constructSaddleSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = false; + scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true; constructBaseScenarioFunctions[BaseSimulationScenario::S] = &constructSSimulationScenario; - scenarioIsSymmetrical[BaseSimulationScenario::S] = false; + scenarioIsSymmetrical[BaseSimulationScenario::S] = true; } void ReducedModelOptimizer::initializePatterns( PatternGeometry &fullPattern, - PatternGeometry &reducedPattern, + ReducedModel &reducedModel, const std::array &optimizationParameters) { - assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN()); + assert(fullPattern.VN() == reducedModel.VN() && fullPattern.EN() >= reducedModel.EN()); -#if POLYSCOPE_DEFINED +#ifdef POLYSCOPE_DEFINED polyscope::removeAllStructures(); #endif // Create copies of the input models PatternGeometry copyFullPattern; - PatternGeometry copyReducedPattern; + ReducedModel copyReducedModel; copyFullPattern.copy(fullPattern); - copyReducedPattern.copy(reducedPattern); + copyReducedModel.copy(reducedModel); +#ifdef POLYSCOPE_DEFINED copyFullPattern.updateEigenEdgeAndVertices(); - copyReducedPattern.updateEigenEdgeAndVertices(); - baseTriangle = copyReducedPattern.getBaseTriangle(); +#endif + copyReducedModel.updateEigenEdgeAndVertices(); + baseTriangle = copyReducedModel.getBaseTriangle(); - computeMaps(copyFullPattern, copyReducedPattern); + computeMaps(copyFullPattern, copyReducedModel); optimizationState.fullPatternOppositeInterfaceViPairs = m_fullPatternOppositeInterfaceViPairs; - createSimulationMeshes(copyFullPattern, copyReducedPattern); + g_baseScenarioMaxDisplacementHelperData.fullPatternOppositeInterfaceViPairs + = m_fullPatternOppositeInterfaceViPairs; + createSimulationMeshes(copyFullPattern, copyReducedModel); initializeUpdateReducedPatternFunctions(); initializeOptimizationParameters(m_pFullPatternSimulationMesh, optimizationParameters); } @@ -570,9 +589,9 @@ void ReducedModelOptimizer::optimize(ConstPatternGeometry &fullPattern, ReducedModelOptimization::Results &optimizationResults) { //scale pattern and reduced model - fullPattern.scale(optimizationSettings.targetBaseTriangleSize, interfaceNodeIndex); + fullPattern.scale(optimizationSettings.targetBaseTriangleSize); ReducedModel reducedModel; - reducedModel.scale(optimizationSettings.targetBaseTriangleSize, interfaceNodeIndex); + reducedModel.scale(optimizationSettings.targetBaseTriangleSize); auto start = std::chrono::system_clock::now(); optimizationResults.label = fullPattern.getLabel(); optimizationResults.baseTriangleFullPattern.copy(fullPattern); @@ -589,18 +608,23 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() { const vcg::Triangle3 &baseTriangle = this->baseTriangle; optimizationState.functions_updateReducedPatternParameter[R] = - [=](const double &newR, std::shared_ptr &pReducedPatternSimulationMesh) { + [=](const double &newR, std::shared_ptr &pReducedModelSimulationMesh) { + // std::shared_ptr pReducedModel_edgeMesh = std::dynamic_pointer_cast( + // pReducedModelSimulationMesh); + // ReducedModel::updateFannedGeometry_R(newR, pReducedModel_edgeMesh); const CoordType barycentricCoordinates_hexagonBaseTriangleVertex(1 - newR, newR / 2, newR / 2); + // const vcg::Triangle3 &baseTriangle = getBaseTriangle(); const CoordType hexagonBaseTriangleVertexPosition = baseTriangle.cP(0) * barycentricCoordinates_hexagonBaseTriangleVertex[0] + baseTriangle.cP(1) * barycentricCoordinates_hexagonBaseTriangleVertex[1] + baseTriangle.cP(2) * barycentricCoordinates_hexagonBaseTriangleVertex[2]; + // constexpr int fanSize = 6; for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { - pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + pReducedModelSimulationMesh->vert[2 * rotationCounter].P() = vcg::RotationMatrix(PatternGeometry::DefaultNormal, vcg::math::ToRad(60.0 * rotationCounter)) * hexagonBaseTriangleVertexPosition; @@ -608,17 +632,21 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions() }; optimizationState.functions_updateReducedPatternParameter[Theta] = - [](const double &newTheta, std::shared_ptr &pReducedPatternSimulationMesh) { - const CoordType baseTriangleHexagonVertexPosition - = pReducedPatternSimulationMesh->vert[0].cP(); + [](const double &newTheta_degrees, + std::shared_ptr &pReducedModelSimulationMesh) { + // std::shared_ptr pReducedModel_edgeMesh + // = std::dynamic_pointer_cast(pReducedModelSimulationMesh); + // ReducedModel::updateFannedGeometry_theta(newTheta_degrees, pReducedModel_edgeMesh); + const CoordType baseTriangleHexagonVertexPosition = pReducedModelSimulationMesh->vert[0] + .cP(); const CoordType thetaRotatedHexagonBaseTriangleVertexPosition - = vcg::RotationMatrix(PatternGeometry::DefaultNormal, vcg::math::ToRad(newTheta)) + = vcg::RotationMatrix(PatternGeometry::DefaultNormal, + vcg::math::ToRad(newTheta_degrees)) * baseTriangleHexagonVertexPosition; - for (int rotationCounter = 0; rotationCounter < ReducedModelOptimizer::fanSize; rotationCounter++) { - pReducedPatternSimulationMesh->vert[2 * rotationCounter].P() + pReducedModelSimulationMesh->vert[2 * rotationCounter].P() = vcg::RotationMatrix(PatternGeometry::DefaultNormal, vcg::math::ToRad(60.0 * rotationCounter)) * thetaRotatedHexagonBaseTriangleVertexPosition; @@ -814,9 +842,9 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv } #ifdef POLYSCOPE_DEFINED std::cout << std::endl; + m_pReducedPatternSimulationMesh->updateEigenEdgeAndVertices(); #endif m_pReducedPatternSimulationMesh->reset(); - results.fullPatternYoungsModulus = youngsModulus; // Compute obj value per simulation scenario and the raw objective value // updateMeshFunction(optimalX); @@ -844,8 +872,9 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv optimizationState.simulationScenarioIndices.size()); for (int i = 0; i < optimizationState.simulationScenarioIndices.size(); i++) { const int simulationScenarioIndex = optimizationState.simulationScenarioIndices[i]; - SimulationResults reducedModelResults = simulator.executeSimulation( - optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]); + std::shared_ptr &pReducedJob + = optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex]; + SimulationResults reducedModelResults = simulator.executeSimulation(pReducedJob); #ifdef POLYSCOPE_DEFINED #ifdef USE_SCENARIO_WEIGHTS @@ -943,11 +972,12 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv << std::endl; std::cout << std::endl; - // reducedModelResults.registerForDrawing(Colors::reducedDeformed); - // optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing(Colors::fullDeformed); - // polyscope::show(); - // reducedModelResults.unregister(); - // optimizationState.fullPatternResults[simulationScenarioIndex].unregister(); + reducedModelResults.registerForDrawing(Colors::reducedDeformed); + optimizationState.fullPatternResults[simulationScenarioIndex].registerForDrawing( + Colors::fullDeformed); + polyscope::show(); + reducedModelResults.unregister(); + optimizationState.fullPatternResults[simulationScenarioIndex].unregister(); #endif } @@ -968,6 +998,7 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv // results.draw(); } +#ifdef DLIB_DEFINED std::array ReducedModelOptimizer::computeFullPatternMaxSimulationForces( const std::vector &desiredBaseSimulationScenario) const @@ -1028,6 +1059,7 @@ ReducedModelOptimizer::getFullPatternMaxSimulationForces( return fullPatternSimulationScenarioMaxMagnitudes; } +#endif void ReducedModelOptimizer::runOptimization(const Settings &settings, ReducedModelOptimization::Results &results) @@ -1087,8 +1119,8 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, = optimizationState.parametersInitialValue[parameterIndex]; return x[xIndex] * parameterInitialValue; }(); - // std::cout << "Optimization parameter:" << parameterIndex << std::endl; - // std::cout << "New value:" << parameterNewValue << std::endl; + // std::cout << "Optimization parameter:" << parameterIndex << std::endl; + // std::cout << "New value:" << parameterNewValue << std::endl; optimizationState .functions_updateReducedPatternParameter[parameterIndex](parameterNewValue, pMesh); @@ -1141,7 +1173,7 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, #else ens::SA optimizer; #endif - // ens::LBestPSO optimizer; + // ens::LBestPSO optimizer; parameterGroup_optimalResult.y = optimizer.Optimize(optimizationFunction, x); // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { // if (x[xIndex] > optimizationState.xMax[xIndex]) { @@ -1155,7 +1187,6 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, std::cout << "optimal x:" << "\n" << x << std::endl; - std::cout << "Minima:" << minima << std::endl; std::cout << "optimal objective:" << optimizationFunction.Evaluate(x) << std::endl; #endif parameterGroup_optimalResult.x.clear(); @@ -1382,36 +1413,47 @@ void ReducedModelOptimizer::constructDomeSimulationScenario( - job.pMesh->vert[viPair.second].cP()); VectorType momentAxis = vcg::RotationMatrix(VectorType(0, 0, 1), vcg::math::ToRad(90.0)) * interfaceVector.Normalize(); - if (viPairIt == oppositeInterfaceViPairs.begin()) { - job.nodalForcedDisplacements[viPair.first] - = Eigen::Vector3d(-interfaceVector[0], - -interfaceVector[1], - 0) - * 0.001 - * std::abs( - forceMagnitude); //NOTE:Should the forced displacement change relatively to the magnitude? - // * std::abs(forceMagnitude / maxForceMagnitude_dome); - job.nodalForcedDisplacements[viPair.second] = Eigen::Vector3d(interfaceVector[0], - interfaceVector[1], - 0) - * 0.001 * std::abs(forceMagnitude); - // * std::abs(forceMagnitude / maxForceMagnitude_dome); - // CoordType v = (pMesh->vert[viPair.first].cP() - pMesh->vert[viPair.second].cP()) - // ^ CoordType(0, 0, -1).Normalize(); - // nodalForces[viPair.first] = Vector6d({0, 0, 0, v[0], v[1], 0}) * forceMagnitude - // * 0.0001; - // nodalForces[viPair.second] = Vector6d({0, 0, 0, -v[0], -v[1], 0}) * forceMagnitude - // * 0.0001; - } else { - job.nodalExternalForces[viPair.first] - = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude - / 3; - job.nodalExternalForces[viPair.second] - = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) - * forceMagnitude / 3; - job.constrainedVertices[viPair.first] = std::unordered_set{2}; - job.constrainedVertices[viPair.second] = std::unordered_set{2}; - } + // if (viPairIt == oppositeInterfaceViPairs.begin()) { + // // job.nodalForcedDisplacements[viPair.first] = Eigen::Vector3d(-interfaceVector[0], + // // -interfaceVector[1], + // // 0) + // // * 0.01 + // /** 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.01 /** 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; + // job.constrainedVertices[viPair.first] = {0, 1, 2}; + // job.constrainedVertices[viPair.second] = {2}; + // job.constrainedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + // job.constrainedVertices[viPair.second] = std::unordered_set{0, 2}; + // } else { + // job.nodalExternalForces[viPair.first] + // = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude + // / 3; + // job.nodalExternalForces[viPair.second] + // = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) + // * forceMagnitude / 3; + // job.constrainedVertices[viPair.first] = std::unordered_set{2}; + // job.constrainedVertices[viPair.second] = std::unordered_set{2}; + // } + job.nodalExternalForces[viPair.first] + = Vector6d({0, 0, 0, momentAxis[0], momentAxis[1], momentAxis[2]}) * forceMagnitude / 3; + job.nodalExternalForces[viPair.second] + = Vector6d({0, 0, 0, -momentAxis[0], -momentAxis[1], -momentAxis[2]}) * forceMagnitude + / 3; + job.constrainedVertices[viPair.first] = std::unordered_set{0, 1, 2}; + job.constrainedVertices[viPair.second] = std::unordered_set{2}; } } @@ -1564,7 +1606,7 @@ double fullPatternMaxSimulationForceTranslationalObjective(const double &forceMa settings.totalExternalForcesNormPercentageTermination = 1e-2; settings.totalTranslationalKineticEnergyThreshold = 1e-8; settings.viscousDampingFactor = 5e-3; - settings.useKineticDamping = true; + settings.useTranslationalKineticEnergyForKineticDamping = true; // settings.averageResidualForcesCriterionThreshold = 1e-5; // settings.useAverage = true; // settings.totalTranslationalKineticEnergyThreshold = 1e-10; @@ -1616,6 +1658,7 @@ struct ForceMagnitudeOptimization }; #endif +#ifdef DLIB_DEFINED double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( const BaseSimulationScenario &scenario) const { @@ -1733,6 +1776,7 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( return forceMagnitude; } +#endif void ReducedModelOptimizer::computeScenarioWeights( const std::vector &baseSimulationScenarios, @@ -1818,7 +1862,7 @@ std::vector> ReducedModelOptimizer::createFullPat const double forceMagnitudeStep = numberOfSimulationScenarios == 1 ? maxForceMagnitude : (maxForceMagnitude - minForceMagnitude) - / (numberOfSimulationScenarios); + / (numberOfSimulationScenarios - 1); const int baseSimulationScenarioIndexOffset = std::accumulate(simulationScenariosResolution.begin(), simulationScenariosResolution.begin() + baseScenario, @@ -1949,11 +1993,15 @@ void ReducedModelOptimizer::optimize( continue; } } + std::set set_desiredScenarios(desiredBaseSimulationScenarios.begin(), + desiredBaseSimulationScenarios.end()); + std::set set_zeroScenarios(zeroMagnitudeScenarios.begin(), + zeroMagnitudeScenarios.end()); std::vector usedBaseSimulationScenarios; - std::set_difference(desiredBaseSimulationScenarios.begin(), - desiredBaseSimulationScenarios.end(), - zeroMagnitudeScenarios.begin(), - zeroMagnitudeScenarios.end(), + std::set_difference(set_desiredScenarios.begin(), + set_desiredScenarios.end(), + set_zeroScenarios.begin(), + set_zeroScenarios.end(), std::back_inserter(usedBaseSimulationScenarios)); assert(!optimizationSettings.optimizationStrategy.empty()); @@ -1979,11 +2027,16 @@ void ReducedModelOptimizer::optimize( 0)); } - // constexpr bool recomputeForceMagnitudes = true; + // constexpr bool recomputeForceMagnitudes = false; // std::array fullPatternSimulationScenarioMaxMagnitudes // = getFullPatternMaxSimulationForces(usedBaseSimulationScenarios, // optimizationSettings.intermediateResultsDirectoryPath, // recomputeForceMagnitudes); + // for (int baseScenario = Axial; baseScenario != NumberOfBaseSimulationScenarios; baseScenario++) { + // fullPatternSimulationScenarioMaxMagnitudes[baseScenario] + // = std::min(fullPatternSimulationScenarioMaxMagnitudes[baseScenario], + // optimizationSettings.baseScenarioMaxMagnitudes[baseScenario]); + // } // optimizationSettings.baseScenarioMaxMagnitudes = fullPatternSimulationScenarioMaxMagnitudes; std::array fullPatternSimulationScenarioMaxMagnitudes = optimizationSettings.baseScenarioMaxMagnitudes; @@ -1997,47 +2050,60 @@ void ReducedModelOptimizer::optimize( results.baseTriangle = baseTriangle; - DRMSimulationModel::Settings simulationSettings; - simulationSettings.totalExternalForcesNormPercentageTermination = 1e-5; - simulationSettings.maxDRMIterations = 200000; - simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-15; -// simulationSettings.viscousDampingFactor = 5e-3; -// simulationSettings.debugModeStep = 100000; -// simulationSettings.shouldDraw = true; -// simulationSettings.beVerbose = true; -// simulationSettings.useKineticDamping = true; -// simulationSettings.save(std::filesystem::path(std::string(__FILE__)) -// .parent_path() -// .parent_path() -// .append("DefaultSettings") -// .append("DRMSettings") -// .append("defaultDRMSimulationSettings.json")); + DRMSimulationModel::Settings drmSettings; + drmSettings.totalExternalForcesNormPercentageTermination = 1e-5; + drmSettings.maxDRMIterations = 200000; + // drmSettings.totalTranslationalKineticEnergyThreshold = 1e-10; + drmSettings.totalTranslationalKineticEnergyThreshold = 1e-9; + // drmSettings.gradualForcedDisplacementSteps = 20; + // drmSettings.linearGuessForceScaleFactor = 1; + // drmSettings.viscousDampingFactor = 5e-3; + // drmSettings.useTotalRotationalKineticEnergyForKineticDamping = true; + // drmSettings.shouldDraw = true; -// simulationSettings.averageResidualForcesCriterionThreshold = 1e-5; -// simulationSettings.viscousDampingFactor = 1e-3; -// simulationSettings.beVerbose = true; -// simulationSettings.shouldDraw = true; -// simulationSettings.isDebugMode = true; -// simulationSettings.debugModeStep = 100000; -#ifdef POLYSCOPE_DEFINED + // drmSettings.totalExternalForcesNormPercentageTermination = 1e-2; + drmSettings.totalTranslationalKineticEnergyThreshold = 1e-8; + // drmSettings.viscousDampingFactor = 5e-3; + // simulationSettings.viscousDampingFactor = 5e-3; + // simulationSettings.linearGuessForceScaleFactor = 1; + // simulationSettings.gamma = 0.3; + // simulationSettings.xi = 0.9999; + // drmSettings.debugModeStep = 100; + // drmSettings.shouldDraw = true; + // drmSettings.shouldCreatePlots = true; + // drmSettings.beVerbose = true; + // simulationSettings.useKineticDamping = true; + // simulationSettings.save(std::filesystem::path(std::string(__FILE__)) + // .parent_path() + // .parent_path() + // .append("DefaultSettings") + // .append("DRMSettings") + // .append("defaultDRMSimulationSettings.json")); + + // simulationSettings.averageResidualForcesCriterionThreshold = 1e-5; + // simulationSettings.viscousDampingFactor = 1e-3; + // simulationSettings.beVerbose = true; + // simulationSettings.shouldDraw = true; + // simulationSettings.isDebugMode = true; + // simulationSettings.debugModeStep = 100000; + //#ifdef POLYSCOPE_DEFINED constexpr bool drawFullPatternSimulationResults = false; - if (drawFullPatternSimulationResults) { - optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( - ReducedModelOptimization::Colors::fullInitial); - } -#endif + auto t1 = std::chrono::high_resolution_clock::now(); + //#endif results.wasSuccessful = true; optimizationState.fullPatternResults.resize(totalNumberOfSimulationScenarios); optimizationState.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios); std::for_each( + //#ifndef POLYSCOPE_DEFINED std::execution::par_unseq, + //#endif optimizationState.simulationScenarioIndices.begin(), optimizationState.simulationScenarioIndices.end(), [&](const int &simulationScenarioIndex) { - // for (int simulationScenarioIndex : optimizationState.simulationScenarioIndices) { const std::shared_ptr &pFullPatternSimulationJob = optimizationState.fullPatternSimulationJobs[simulationScenarioIndex]; + // std::cout << "Executing " << pFullPatternSimulationJob->getLabel() << std::endl; std::filesystem::path jobResultsDirectoryPath( std::filesystem::path(optimizationSettings.intermediateResultsDirectoryPath) @@ -2050,7 +2116,7 @@ void ReducedModelOptimizer::optimize( //#else // constexpr bool recomputeFullPatternResults = true; //#endif - SimulationResults fullPatternResults; + SimulationResults patternDrmResults; // if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) { // fullPatternResults // .load(std::filesystem::path(jobResultsDirectoryPath).append("Results"), @@ -2078,76 +2144,101 @@ void ReducedModelOptimizer::optimize( //#else DRMSimulationModel simulator; - fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, - simulationSettings); -// fullPatternResults.save(jobResultsDirectoryPath); -// } -//#endif -// if (!fullPatternResults.converged) { -// DRMSimulationModel::Settings simulationSettings_secondRound; -// simulationSettings_secondRound.viscousDampingFactor = 2e-3; -// simulationSettings_secondRound.useKineticDamping = true; -// simulationSettings.maxDRMIterations = 200000; -// fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, -// simulationSettings_secondRound); -#ifdef POLYSCOPE_DEFINED -// std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel() -// << " used viscous damping." << std::endl; -#endif + patternDrmResults = simulator.executeSimulation(pFullPatternSimulationJob, drmSettings); + // fullPatternResults.save(jobResultsDirectoryPath); + // } + //#endif + // if (!fullPatternResults.converged) { + // DRMSimulationModel::Settings simulationSettings_secondRound; + // simulationSettings_secondRound.viscousDampingFactor = 2e-3; + // simulationSettings_secondRound.useKineticDamping = true; + // simulationSettings.maxDRMIterations = 200000; + // fullPatternResults = simulator.executeSimulation(pFullPatternSimulationJob, + // simulationSettings_secondRound); - if (!fullPatternResults.converged) { - results.wasSuccessful = false; - std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() - << " did not converge." << std::endl; + if (!patternDrmResults.converged) { #ifdef POLYSCOPE_DEFINED - // DRMSimulationModel::Settings debugSimulationSettings; - // debugSimulationSettings.debugModeStep = 50; - // // // debugSimulationSettings.maxDRMIterations = 100000; - // debugSimulationSettings.shouldDraw = true; - // // debugSimulationSettings.drawingStep = debugSimulationSettings.debugModeStep; - // debugSimulationSettings.shouldCreatePlots = true; - // // // debugSimulationSettings.Dtini = 0.06; - // debugSimulationSettings.beVerbose = true; - // debugSimulationSettings.averageResidualForcesCriterionThreshold = 1e-5; - // debugSimulationSettings.maxDRMIterations = 100000; - // debugSimulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8; - // debugSimulationSettings.viscousDampingFactor = 1e-2; - // // // debugSimulationSettings.totalExternalForcesNormPercentageTermination = 1e-3; - // // // debugSimulationSettings.shouldUseTranslationalKineticEnergyThreshold = true; - // auto debugResults = simulator.executeSimulation(pFullPatternSimulationJob, - // debugSimulationSettings); - // debugResults.setLabelPrefix("debugResults"); - // debugResults.registerForDrawing(); - // polyscope::show(); - // debugResults.unregister(); std::filesystem::path outputPath( std::filesystem::path("../nonConvergingJobs") .append(m_pFullPatternSimulationMesh->getLabel()) .append("final_" + pFullPatternSimulationJob->getLabel())); std::filesystem::create_directories(outputPath); pFullPatternSimulationJob->save(outputPath); - simulationSettings.save(outputPath); - std::terminate(); - return; +// drmSettings_secondTry.save(outputPath); #endif - // } - } + //#ifdef POLYSCOPE_DEFINED + // SimulationResults fullPatternResults_linear; + // if (drawFullPatternSimulationResults) { + // LinearSimulationModel linearSimulator; + // fullPatternResults_linear = linearSimulator.executeSimulation( + // pFullPatternSimulationJob); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear.registerForDrawing(std::array{0, 0, 255}, + // true); + // std::cout << pFullPatternSimulationJob->getLabel() + // << " did not converge on the first try. Trying again.." << std::endl; + // } + //#endif + DRMSimulationModel::Settings drmSettings_secondTry = drmSettings; + // drmSettings_secondTry.shouldDraw = true; + // drmSettings_secondTry.debugModeStep = 20000; + // drmSettings_secondTry.shouldCreatePlots = true; + // drmSettings_secondTry.beVerbose = true; + drmSettings_secondTry.linearGuessForceScaleFactor = 1; + drmSettings_secondTry.viscousDampingFactor = 4e-8; + *drmSettings_secondTry.maxDRMIterations *= 5; + drmSettings_secondTry.useTotalRotationalKineticEnergyForKineticDamping = true; + // drmSettings.totalTranslationalKineticEnergyThreshold = std::nullopt; + // drmSettings.totalExternalForcesNormPercentageTermination = 1e-3; + drmSettings.totalTranslationalKineticEnergyThreshold = 1e-10; + // drmSettings_secondTry.gamma = 2; + // drmSettings_secondTry.xi = 0.98; + SimulationResults drmResults_secondTry + = simulator.executeSimulation(pFullPatternSimulationJob, drmSettings_secondTry); + + if (drmResults_secondTry.converged) { #ifdef POLYSCOPE_DEFINED - if (drawFullPatternSimulationResults) { - // SimulationResults fullPatternResults_linear = linearSimulator.executeSimulation( - // pFullPatternSimulationJob); - fullPatternResults.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, - true); - // fullPatternResults_linear.labelPrefix += "_linear"; - // fullPatternResults_linear.registerForDrawing(ReducedModelOptimization::Colors::fullDeformed, - // true, - // true); - polyscope::show(); - fullPatternResults.unregister(); - // fullPatternResults_linear.unregister(); - } + std::cout << "Simulation job " << pFullPatternSimulationJob->getLabel() + << " converged after a second try." << std::endl; #endif - optimizationState.fullPatternResults[simulationScenarioIndex] = fullPatternResults; + patternDrmResults = drmResults_secondTry; + } else { + results.wasSuccessful = false; + std::cerr << "Simulation job " << pFullPatternSimulationJob->getLabel() + << " of pattern " << pFullPatternSimulationJob->pMesh->getLabel() + << " did not converge." << std::endl; +#ifdef POLYSCOPE_DEFINED + std::filesystem::path outputPath( + std::filesystem::path("../nonConvergingJobs") + .append(m_pFullPatternSimulationMesh->getLabel()) + .append("final_" + pFullPatternSimulationJob->getLabel())); + std::filesystem::create_directories(outputPath); + pFullPatternSimulationJob->save(outputPath); + drmSettings_secondTry.save(outputPath); + + if (drawFullPatternSimulationResults) { + optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( + ReducedModelOptimization::Colors::fullInitial); + pFullPatternSimulationJob->registerForDrawing( + optimizationState.fullPatternSimulationJobs[0]->pMesh->getLabel(), true); + patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); + // SimulationResults fullPatternResults_linear + // = linearSimulator.executeSimulation(pFullPatternSimulationJob); + // patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear + // .registerForDrawing(std::array{0, 0, 255}, true); + polyscope::show(); + patternDrmResults.unregister(); + // fullPatternResults_linear.unregister(); + optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); + } +#endif + return; + } + } + + optimizationState.fullPatternResults[simulationScenarioIndex] = patternDrmResults; SimulationJob reducedPatternSimulationJob; reducedPatternSimulationJob.pMesh = m_pReducedPatternSimulationMesh; computeReducedModelSimulationJob(*pFullPatternSimulationJob, @@ -2155,13 +2246,42 @@ void ReducedModelOptimizer::optimize( reducedPatternSimulationJob); optimizationState.reducedPatternSimulationJobs[simulationScenarioIndex] = std::make_shared(reducedPatternSimulationJob); - // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; + // std::cout << "Ran sim scenario:" << simulationScenarioIndex << std::endl; +#ifdef POLYSCOPE_DEFINED + if (drawFullPatternSimulationResults) { + optimizationState.fullPatternSimulationJobs[0]->pMesh->registerForDrawing( + ReducedModelOptimization::Colors::fullInitial); + pFullPatternSimulationJob->registerForDrawing(optimizationState + .fullPatternSimulationJobs[0] + ->pMesh->getLabel(), + true); + patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); + // SimulationResults fullPatternResults_linear + // = linearSimulator.executeSimulation(pFullPatternSimulationJob); + // patternDrmResults.registerForDrawing(std::array{0, 255, 0}, true); + // fullPatternResults_linear.labelPrefix += "_linear"; + // fullPatternResults_linear + // .registerForDrawing(std::array{0, 0, 255}, true); + polyscope::show(); + patternDrmResults.unregister(); + // fullPatternResults_linear.unregister(); + optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); + } +#endif }); + auto t2 = std::chrono::high_resolution_clock::now(); + auto s_int = std::chrono::duration_cast(t2 - t1); + std::cout << s_int.count() << std::endl; + results.wasSuccessful = false; + return; #ifdef POLYSCOPE_DEFINED + std::cout << "Computed ground of truth pattern results in:" << s_int.count() << " seconds." + << std::endl; if (drawFullPatternSimulationResults) { optimizationState.fullPatternSimulationJobs[0]->pMesh->unregister(); } + optimizationState.reducedPatternSimulationJobs[0]->pMesh->updateEigenEdgeAndVertices(); #endif // if (optimizationState.optimizationSettings.normalizationStrategy // == Settings::NormalizationStrategy::Epsilon) { @@ -2171,5 +2291,4 @@ void ReducedModelOptimizer::optimize( std::cout << "Running reduced model optimization" << std::endl; #endif runOptimization(optimizationSettings, results); - results.notes = optimizationNotes; } diff --git a/reducedmodeloptimizer.hpp b/reducedmodeloptimizer.hpp index 26644ca..0ba6e50 100644 --- a/reducedmodeloptimizer.hpp +++ b/reducedmodeloptimizer.hpp @@ -117,7 +117,7 @@ public: inline static std::array parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"}; constexpr static std::array - simulationScenariosResolution = {11, 11, 21, 21, 21, 21}; + simulationScenariosResolution = {12, 12, 22, 22, 22, 22}; constexpr static std::array baseScenarioWeights = {1, 1, 2, 3, 2}; inline static int totalNumberOfSimulationScenarios @@ -150,7 +150,7 @@ public: std::shared_ptr &pReducedPatternSimulationMesh); void computeMaps(const std::unordered_map> &slotToNode, PatternGeometry &fullPattern, - PatternGeometry &reducedPattern, + ReducedModel &reducedPattern, std::unordered_map &reducedToFullInterfaceViMap, std::unordered_map @@ -263,7 +263,7 @@ public: ReducedModelOptimization::Results optimize( ConstPatternGeometry &fullPattern, - const ReducedModelOptimization::Settings &optimizatifullyonSettings); + const ReducedModelOptimization::Settings &optimizationSettings); void optimize(ConstPatternGeometry &fullPattern, ReducedModelOptimization::Settings &optimizationSettings, @@ -285,7 +285,7 @@ private: ReducedModelOptimization::S})); void initializePatterns(PatternGeometry &fullPattern, - PatternGeometry &reducedPattern, + ReducedModel &reducedPattern, const std::array &optimizationParameters); @@ -295,7 +295,7 @@ private: Eigen::MatrixX3d &optimalDisplacementsOfReducedModel); void runOptimization(const ReducedModelOptimization::Settings &settings, ReducedModelOptimization::Results &results); - void computeMaps(PatternGeometry &fullModel, PatternGeometry &reducedPattern); + void computeMaps(PatternGeometry &fullModel, ReducedModel &reducedModel); void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel); void initializeOptimizationParameters( const std::shared_ptr &mesh,