From 107a2ce66c6abf8599505fe82e621096cf5757ca Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Tue, 30 Mar 2021 12:35:00 +0300 Subject: [PATCH] Objective value takes into account the rotational displacemenets between the full and the reduced pattern. --- src/main.cpp | 4 +- src/reducedmodeloptimizer.cpp | 104 ++++++++++++++++++++++++---------- 2 files changed, 75 insertions(+), 33 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index fe4cdbb..f1ddb01 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -56,7 +56,7 @@ int main(int argc, char *argv[]) { settings_optimization.normalizationStrategy = ReducedModelOptimization::Settings::NormalizationStrategy::Epsilon; settings_optimization.normalizationParameter = 0.0003; - settings_optimization.solutionAccuracy = 0.01; + settings_optimization.solutionAccuracy = 0.001; // Optimize pair const std::string pairName = @@ -97,7 +97,7 @@ int main(int argc, char *argv[]) { // .string(), // false); csv_results << "Name"; - optimizationResults.writeHeaderTo(settings_optimization, csv_results); + optimizationResults.writeHeaderTo(csv_results); settings_optimization.writeHeaderTo(csv_results); csv_results << endrow; csv_results << pairName; diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index fb0f687..7ab6843 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -13,7 +13,8 @@ struct GlobalOptimizationVariables { std::vector g_optimalReducedModelDisplacements; // std::vector> fullPatternDisplacements; std::vector fullPatternResults; - std::vector objectiveNormalizationValues; + std::vector translationalDisplacementNormalizationValues; + std::vector rotationalDisplacementNormalizationValues; std::vector> fullPatternSimulationJobs; std::vector> reducedPatternSimulationJobs; std::unordered_map @@ -67,12 +68,6 @@ double ReducedModelOptimizer::computeRawDisplacementError( reducedPatternDisplacements[reducedModelVi][0], reducedPatternDisplacements[reducedModelVi][1], reducedPatternDisplacements[reducedModelVi][2]); - if (!std::isfinite(reducedVertexDisplacement[0]) || - !std::isfinite(reducedVertexDisplacement[1]) || - !std::isfinite(reducedVertexDisplacement[2])) { - std::cout << "Displacements are not finite" << std::endl; - std::terminate(); - } const VertexIndex fullModelVi = reducedFullViPair.second; Eigen::Vector3d fullVertexDisplacement( fullPatternDisplacements[fullModelVi][0], @@ -134,11 +129,24 @@ double ReducedModelOptimizer::objective(long n, const double *x) { reducedModelResults.internalPotentialEnergy - global.fullPatternResults[simulationScenarioIndex].internalPotentialEnergy); } else { - simulationScenarioError = computeDisplacementError( + const double translationalError = computeDisplacementError( reducedModelResults.displacements, global.fullPatternResults[simulationScenarioIndex].displacements, global.reducedToFullInterfaceViMap, - global.objectiveNormalizationValues[simulationScenarioIndex]); + global.translationalDisplacementNormalizationValues[simulationScenarioIndex]); + double rotationalError = 0; + for (const auto &reducedToFullInterfaceViPair : global.reducedToFullInterfaceViMap) { + const double vertexRotationalError + = global.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion[reducedToFullInterfaceViPair.second] + .angularDistance(reducedModelResults.rotationalDisplacementQuaternion + [reducedToFullInterfaceViPair.first]); + rotationalError += vertexRotationalError; + } + simulationScenarioError += translationalError + + rotationalError + / global.rotationalDisplacementNormalizationValues + [simulationScenarioIndex]; } //#ifdef POLYSCOPE_DEFINED // std::cout << "sim error:" << simulationScenarioError << std::endl; @@ -549,10 +557,9 @@ void ReducedModelOptimizer::visualizeResults( std::shared_ptr pFullPatternSimulationMesh = fullPatternSimulationJobs[0]->pMesh; pFullPatternSimulationMesh->registerForDrawing(); - pFullPatternSimulationMesh->savePly(pFullPatternSimulationMesh->getLabel() + - "_undeformed.ply"); - reducedPatternSimulationJobs[0]->pMesh->savePly( - reducedPatternSimulationJobs[0]->pMesh->getLabel() + "_undeformed.ply"); + 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 = @@ -570,8 +577,8 @@ void ReducedModelOptimizer::visualizeResults( double normalizationFactor = 1; if (global.optimizationSettings.normalizationStrategy != ReducedModelOptimization::Settings::NormalizationStrategy::NonNormalized) { - normalizationFactor = - global.objectiveNormalizationValues[simulationScenarioIndex]; + normalizationFactor + = global.translationalDisplacementNormalizationValues[simulationScenarioIndex]; } reducedModelResults.saveDeformedModel(); fullModelResults.saveDeformedModel(); @@ -680,10 +687,10 @@ ReducedModelOptimizer::runOptimization(const Settings &settings) { std::chrono::duration_cast(end - start); ReducedModelOptimization::Results results; results.numberOfSimulationCrashes = global.numOfSimulationCrashes; - results.optimalXNameValuePairs.reserve(settings.xRanges.size()); + results.optimalXNameValuePairs.resize(settings.xRanges.size()); for(int xVariableIndex=0;xVariableIndex fullPatternDisplacementNormSum( - NumberOfSimulationScenarios); + std::vector fullPatternTranslationalDisplacementNormSum(NumberOfSimulationScenarios); + std::vector fullPatternAngularDistance(NumberOfSimulationScenarios); for (int simulationScenarioIndex : global.simulationScenarioIndices) { - double displacementNormSum = 0; + 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 + if (global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.contains(fullPatternVi)) { + const std::unordered_set constrainedDof + = global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.at(fullPatternVi); + if (constrainedDof.contains(0) && constrainedDof.contains(1) + && constrainedDof.contains(2)) { + continue; + } + } + const Vector6d &vertexDisplacement = global .fullPatternResults[simulationScenarioIndex] - .displacements[interfaceViPair.second]; - displacementNormSum += vertexDisplacement.getTranslation().norm(); + .displacements[fullPatternVi]; + translationalDisplacementNormSum += vertexDisplacement.getTranslation().norm(); + //If the full pattern vertex is rotationally constrained dont take it into account + if (global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.contains(fullPatternVi)) { + const std::unordered_set constrainedDof + = global.fullPatternSimulationJobs[simulationScenarioIndex] + ->constrainedVertices.at(fullPatternVi); + if (constrainedDof.contains(3) && constrainedDof.contains(5) + && constrainedDof.contains(4)) { + continue; + } + } + angularDistanceSum += global.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion[fullPatternVi] + .angularDistance(Eigen::Quaterniond::Identity()); } - fullPatternDisplacementNormSum[simulationScenarioIndex] = - displacementNormSum; + fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] + = translationalDisplacementNormSum; + fullPatternAngularDistance[simulationScenarioIndex] = angularDistanceSum; } for (int simulationScenarioIndex : global.simulationScenarioIndices) { if (global.optimizationSettings.normalizationStrategy == Settings::NormalizationStrategy::Epsilon) { - const double epsilon = - global.optimizationSettings.normalizationParameter; - global.objectiveNormalizationValues[simulationScenarioIndex] = std::max( - fullPatternDisplacementNormSum[simulationScenarioIndex], epsilon); + const double epsilon_translationalDisplacement = global.optimizationSettings + .normalizationParameter; + global.translationalDisplacementNormalizationValues[simulationScenarioIndex] + = std::max(fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex], + epsilon_translationalDisplacement); + const double epsilon_rotationalDisplacement = vcg::math::ToRad(2.0); + global.rotationalDisplacementNormalizationValues[simulationScenarioIndex] + = std::max(fullPatternAngularDistance[simulationScenarioIndex], + epsilon_rotationalDisplacement); } } } @@ -919,7 +960,8 @@ Results ReducedModelOptimizer::optimize( global.g_optimalReducedModelDisplacements.resize(NumberOfSimulationScenarios); global.reducedPatternSimulationJobs.resize(NumberOfSimulationScenarios); global.fullPatternResults.resize(NumberOfSimulationScenarios); - global.objectiveNormalizationValues.resize(NumberOfSimulationScenarios); + global.translationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios); + global.rotationalDisplacementNormalizationValues.resize(NumberOfSimulationScenarios); global.minY = std::numeric_limits::max(); global.numOfSimulationCrashes = 0; global.numberOfFunctionCalls = 0;