diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..8708fe2 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "deps/MySources"] + path = deps/MySources + url = https://gitea-s2i2s.isti.cnr.it/manolas/MySources diff --git a/deps/MySources b/deps/MySources new file mode 160000 index 0000000..f79e0ac --- /dev/null +++ b/deps/MySources @@ -0,0 +1 @@ +Subproject commit f79e0acb4f3361939cc495956b7afb52a318d02b diff --git a/src/reducedmodeloptimizer.cpp b/src/reducedmodeloptimizer.cpp index 3b10f31..314023f 100644 --- a/src/reducedmodeloptimizer.cpp +++ b/src/reducedmodeloptimizer.cpp @@ -108,9 +108,9 @@ double ReducedModelOptimizer::computeRawRotationalError( { double rawRotationalError = 0; for (const auto &reducedToFullInterfaceViPair : reducedToFullInterfaceViMap) { - const double vertexRotationalError - = rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance( - rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first]); + const double vertexRotationalError = std::abs( + rotatedQuaternion_fullPattern[reducedToFullInterfaceViPair.second].angularDistance( + rotatedQuaternion_reducedPattern[reducedToFullInterfaceViPair.first])); rawRotationalError += vertexRotationalError; } @@ -140,27 +140,27 @@ double ReducedModelOptimizer::computeError( const double &scenarioWeight, const Settings::ObjectiveWeights &objectiveWeights) { - const double translationalError - = computeDisplacementError(simulationResults_fullPattern.displacements, - simulationResults_reducedPattern.displacements, - reducedToFullInterfaceViMap, - normalizationFactor_translationalDisplacement); // const double translationalError - // = computeRawTranslationalError(simulationResults_fullPattern.displacements, - // simulationResults_reducedPattern.displacements, - // reducedToFullInterfaceViMap); + // = computeDisplacementError(simulationResults_fullPattern.displacements, + // simulationResults_reducedPattern.displacements, + // reducedToFullInterfaceViMap, + // normalizationFactor_translationalDisplacement); + const double translationalError + = computeRawTranslationalError(simulationResults_fullPattern.displacements, + simulationResults_reducedPattern.displacements, + reducedToFullInterfaceViMap); // std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl; // std::cout << "trans error:" << translationalError << std::endl; - const double rotationalError - = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, - simulationResults_reducedPattern.rotationalDisplacementQuaternion, - reducedToFullInterfaceViMap, - normalizationFactor_rotationalDisplacement); // const double rotationalError - // = computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, - // simulationResults_reducedPattern.rotationalDisplacementQuaternion, - // reducedToFullInterfaceViMap); + // = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, + // simulationResults_reducedPattern.rotationalDisplacementQuaternion, + // reducedToFullInterfaceViMap, + // normalizationFactor_rotationalDisplacement); + const double rotationalError + = computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion, + simulationResults_reducedPattern.rotationalDisplacementQuaternion, + reducedToFullInterfaceViMap); // std::cout << "rot error:" << rotationalError<< std::endl; const double simulationError = objectiveWeights.translational * translationalError + objectiveWeights.rotational * rotationalError; @@ -186,18 +186,16 @@ struct EnsmallenOptimizationObjective //#ifdef POLYSCOPE_DEFINED // std::cout << "Out of range" << std::endl; //#endif - // return std::numeric_limits::max(); - // return 1e7; - // return 1000; - x[xi] = global.xMax[xi]; + return std::numeric_limits::max(); + // return 1e10; + // x[xi] = global.xMax[xi]; } else if (x[xi] < global.xMin[xi]) { //#ifdef POLYSCOPE_DEFINED // std::cout << "Out of range" << std::endl; //#endif - // return std::numeric_limits::max(); - // return 1e7; - // return 1000; - x[xi] = global.xMin[xi]; + return std::numeric_limits::max(); + // return 1e10; + // x[xi] = global.xMin[xi]; } } return ReducedModelOptimizer::objective(x); @@ -314,9 +312,9 @@ double ReducedModelOptimizer::objective(const std::vector &x) // double totalError = std::accumulate(simulationErrorsPerScenario.begin(), // simulationErrorsPerScenario.end(), // 0); - double totalError = std::reduce(std::execution::par_unseq, - simulationErrorsPerScenario.begin(), - simulationErrorsPerScenario.end()); + const double totalError = std::reduce(std::execution::par_unseq, + simulationErrorsPerScenario.begin(), + simulationErrorsPerScenario.end()); // std::cout << totalError << std::endl; // global.objectiveValueHistory.push_back(totalError); // global.plotColors.push_back(10); @@ -1095,16 +1093,16 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, // that is able to handle arbitrary functions. EnsmallenOptimizationObjective optimizationFunction; //Set min max values - // ens::SA optimizer; - // ens::CNE optimizer; - // ens::DE optimizer; - // ens::SPSA optimizer; - // arma::mat xMin_arma(global.xMin); - // arma::mat xMax_arma(global.xMax); + // ens::SA optimizer; + // ens::CNE optimizer; + // ens::DE optimizer; + // ens::SPSA optimizer; + arma::mat xMin_arma(global.xMin); + arma::mat xMax_arma(global.xMax); // ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000); - ens::LBestPSO optimizer(64, - 1, - 1, + ens::LBestPSO optimizer(200, + xMin_arma, + xMax_arma, settings.numberOfFunctionCalls, 500, settings.solverAccuracy, @@ -1112,13 +1110,13 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings, 2.35); // ens::LBestPSO optimizer; const double minima = optimizer.Optimize(optimizationFunction, x); - for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { - if (x[xIndex] > global.xMax[xIndex]) { - x[xIndex] = global.xMax[xIndex]; - } else if (x[xIndex] < global.xMin[xIndex]) { - x[xIndex] = global.xMin[xIndex]; - } - } + // for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) { + // if (x[xIndex] > global.xMax[xIndex]) { + // x[xIndex] = global.xMax[xIndex]; + // } else if (x[xIndex] < global.xMin[xIndex]) { + // x[xIndex] = global.xMin[xIndex]; + // } + // } #ifdef POLYSCOPE_DEFINED std::cout << "optimal x:" @@ -1327,14 +1325,14 @@ void ReducedModelOptimizer::constructDomeSimulationScenario( = Eigen::Vector3d(-interfaceVector[0], -interfaceVector[1], 0) - * 0.005 + * 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.005 * std::abs(forceMagnitude); + * 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(); @@ -1565,19 +1563,22 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce( global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight; break; case Shear: - global.desiredMaxDisplacementValue = 0.03 * baseTriangleHeight; + global.desiredMaxDisplacementValue = 0.04 * baseTriangleHeight; break; case Bending: + global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; + global.desiredMaxDisplacementValue = 0.05 * baseTriangleHeight; break; case Dome: - global.desiredMaxRotationAngle = vcg::math::ToRad(35.0); + global.desiredMaxRotationAngle = vcg::math::ToRad(20.0); objectiveFunction = &fullPatternMaxSimulationForceRotationalObjective; forceMagnitudeEpsilon *= 1e-2; - objectiveEpsilon = vcg::math::ToRad(3.0); + objectiveEpsilon = vcg::math::ToRad(1.0); forceMagnitude = 0.6; break; case Saddle: global.interfaceViForComputingScenarioError = global.fullPatternInterfaceViPairs[0].first; + global.desiredMaxDisplacementValue = 0.05 * baseTriangleHeight; break; } @@ -1801,9 +1802,9 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors() continue; } } - angularDistanceSum += global.fullPatternResults[simulationScenarioIndex] - .rotationalDisplacementQuaternion[fullPatternVi] - .angularDistance(Eigen::Quaterniond::Identity()); + angularDistanceSum += std::abs(global.fullPatternResults[simulationScenarioIndex] + .rotationalDisplacementQuaternion[fullPatternVi] + .angularDistance(Eigen::Quaterniond::Identity())); } fullPatternTranslationalDisplacementNormSum[simulationScenarioIndex] diff --git a/src/reducedmodeloptimizer.hpp b/src/reducedmodeloptimizer.hpp index d370b86..7de30d0 100644 --- a/src/reducedmodeloptimizer.hpp +++ b/src/reducedmodeloptimizer.hpp @@ -66,7 +66,7 @@ public: // inline constexpr static ParameterLabels parameterLabels(); inline static std::array - parameterLabels = {"R", "A", "I2", "I3", "J", "Theta", "R"}; + parameterLabels = {"E", "A", "I2", "I3", "J", "Theta", "R"}; constexpr static std::array simulationScenariosResolution = {11, 11, 20, 20, 20}; constexpr static std::array