diff --git a/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp b/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp index a985ff46abcf..2e1e5183b4b4 100644 --- a/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp +++ b/applications/GeoMechanicsApplication/custom_python/add_custom_strategies_to_python.cpp @@ -110,18 +110,18 @@ void AddCustomStrategiesToPython(pybind11::module& m) py::class_( m, "GeoMechanicsNewtonRaphsonStrategy") - .def(py::init()); py::class_( m, "GeoMechanicsNewtonRaphsonErosionProcessStrategy") - .def(py::init()); py::class_( m, "GeoMechanicsRammArcLengthStrategy") - .def(py::init()) .def("UpdateLoads", &GeoMechanicsRammArcLengthStrategyType::UpdateLoads); diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp index 2431333c34f5..aca6fca4a05d 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_erosion_process_strategy.hpp @@ -54,7 +54,6 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy GeoMechanicsNewtonRaphsonErosionProcessStrategy(ModelPart& model_part, typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -63,16 +62,7 @@ class GeoMechanicsNewtonRaphsonErosionProcessStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : GeoMechanicsNewtonRaphsonStrategy( - model_part, - pScheme, - pNewLinearSolver, - pNewConvergenceCriteria, - pNewBuilderAndSolver, - rParameters, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { rank = model_part.GetCommunicator().MyPID(); mPipingIterations = rParameters["max_piping_iterations"].GetInt(); diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp index f8360def0513..c9017fe3965d 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_newton_raphson_strategy.hpp @@ -52,16 +52,14 @@ class GeoMechanicsNewtonRaphsonStrategy * @brief Default constructor * @param rModelPart The model part of the problem * @param pScheme The integration scheme - * @param pNewLinearSolver The linear solver employed * @param pNewConvergenceCriteria The convergence criteria employed * @param MaxIterations The maximum number of iterations * @param CalculateReactions The flag for the reaction calculation * @param ReformDofSetAtEachStep The flag that allows to compute the modification of the DOF * @param MoveMeshFlag The flag that allows to move the mesh */ - GeoMechanicsNewtonRaphsonStrategy(ModelPart& model_part, - typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, + GeoMechanicsNewtonRaphsonStrategy(ModelPart& model_part, + typename TSchemeType::Pointer pScheme, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -70,15 +68,7 @@ class GeoMechanicsNewtonRaphsonStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : ResidualBasedNewtonRaphsonStrategy( - model_part, - pScheme, - /*pNewLinearSolver,*/ - pNewConvergenceCriteria, - pNewBuilderAndSolver, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { // only include validation with c++11 since raw_literals do not exist in c++03 Parameters default_parameters(R"( @@ -100,62 +90,6 @@ class GeoMechanicsNewtonRaphsonStrategy // Validate against defaults -- this also ensures no type mismatch rParameters.ValidateAndAssignDefaults(default_parameters); - - // Set Load SubModelParts and Variable names - if (rParameters["loads_sub_model_part_list"].size() > 0) { - mSubModelPartList.resize(rParameters["loads_sub_model_part_list"].size()); - mVariableNames.resize(rParameters["loads_variable_list"].size()); - - if (mSubModelPartList.size() != mVariableNames.size()) - KRATOS_ERROR << "For each SubModelPart there must be a corresponding nodal Variable" - << std::endl; - - for (unsigned int i = 0; i < mVariableNames.size(); ++i) { - mSubModelPartList[i] = - &(model_part.GetSubModelPart(rParameters["loads_sub_model_part_list"][i].GetString())); - mVariableNames[i] = rParameters["loads_variable_list"][i].GetString(); - } - } - } - -protected: - /// Member Variables - std::vector mSubModelPartList; /// List of every SubModelPart associated to an external load - std::vector mVariableNames; /// Name of the nodal variable associated to every SubModelPart - - int Check() override - { - KRATOS_TRY - - return MotherType::Check(); - - KRATOS_CATCH("") - } - - double CalculateReferenceDofsNorm(DofsArrayType& rDofSet) - { - double ReferenceDofsNorm = 0.0; - - int NumThreads = ParallelUtilities::GetNumThreads(); - OpenMPUtils::PartitionVector DofSetPartition; - OpenMPUtils::DivideInPartitions(rDofSet.size(), NumThreads, DofSetPartition); - -#pragma omp parallel reduction(+ : ReferenceDofsNorm) - { - int k = OpenMPUtils::ThisThread(); - - typename DofsArrayType::iterator DofsBegin = rDofSet.begin() + DofSetPartition[k]; - typename DofsArrayType::iterator DofsEnd = rDofSet.begin() + DofSetPartition[k + 1]; - - for (typename DofsArrayType::iterator itDof = DofsBegin; itDof != DofsEnd; ++itDof) { - if (itDof->IsFree()) { - const double& temp = itDof->GetSolutionStepValue(); - ReferenceDofsNorm += temp * temp; - } - } - } - - return sqrt(ReferenceDofsNorm); } }; // Class GeoMechanicsNewtonRaphsonStrategy diff --git a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp index 2ecd12520b26..e13287c0bcaf 100644 --- a/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp +++ b/applications/GeoMechanicsApplication/custom_strategies/strategies/geo_mechanics_ramm_arc_length_strategy.hpp @@ -19,6 +19,9 @@ // Application includes #include "geo_mechanics_application_variables.h" +#include +#include + namespace Kratos { @@ -50,12 +53,9 @@ class GeoMechanicsRammArcLengthStrategy using GrandMotherType::mpDx; // Delta x of iteration i using GrandMotherType::mpScheme; using GrandMotherType::mReformDofSetAtEachStep; - using MotherType::mSubModelPartList; - using MotherType::mVariableNames; - GeoMechanicsRammArcLengthStrategy(ModelPart& model_part, - typename TSchemeType::Pointer pScheme, - typename TLinearSolver::Pointer pNewLinearSolver, + GeoMechanicsRammArcLengthStrategy(ModelPart& model_part, + typename TSchemeType::Pointer pScheme, typename TConvergenceCriteriaType::Pointer pNewConvergenceCriteria, typename TBuilderAndSolverType::Pointer pNewBuilderAndSolver, Parameters& rParameters, @@ -64,17 +64,24 @@ class GeoMechanicsRammArcLengthStrategy bool ReformDofSetAtEachStep = false, bool MoveMeshFlag = false) : GeoMechanicsNewtonRaphsonStrategy( - model_part, - pScheme, - pNewLinearSolver, - pNewConvergenceCriteria, - pNewBuilderAndSolver, - rParameters, - MaxIterations, - CalculateReactions, - ReformDofSetAtEachStep, - MoveMeshFlag) + model_part, pScheme, pNewConvergenceCriteria, pNewBuilderAndSolver, rParameters, MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag) { + // Set Load SubModelParts and Variable names + if (rParameters["loads_sub_model_part_list"].size() > 0) { + mSubModelPartList.resize(rParameters["loads_sub_model_part_list"].size()); + mVariableNames.resize(rParameters["loads_variable_list"].size()); + + if (mSubModelPartList.size() != mVariableNames.size()) + KRATOS_ERROR << "For each SubModelPart there must be a corresponding nodal Variable" + << std::endl; + + for (unsigned int i = 0; i < mVariableNames.size(); ++i) { + mSubModelPartList[i] = + &(model_part.GetSubModelPart(rParameters["loads_sub_model_part_list"][i].GetString())); + mVariableNames[i] = rParameters["loads_variable_list"][i].GetString(); + } + } + mDesiredIterations = rParameters["desired_iterations"].GetInt(); mMaxRadiusFactor = rParameters["max_radius_factor"].GetDouble(); mMinRadiusFactor = rParameters["min_radius_factor"].GetDouble(); @@ -306,7 +313,7 @@ class GeoMechanicsRammArcLengthStrategy else if (mRadius < mMinRadiusFactor * mRadius_0) mRadius = mMinRadiusFactor * mRadius_0; // Update Norm of x - mNormxEquilibrium = this->CalculateReferenceDofsNorm(rDofSet); + mNormxEquilibrium = CalculateReferenceDofsNorm(rDofSet); } else { std::cout << "************ NO CONVERGENCE: restoring equilibrium path ************" << std::endl; @@ -534,6 +541,23 @@ class GeoMechanicsRammArcLengthStrategy // Save the applied Lambda factor mLambda_old = mLambda; } + +private: + static double CalculateReferenceDofsNorm(DofsArrayType& rDofSet) + { + auto is_free_dof = [](const auto& rDof) { return rDof.IsFree(); }; + auto free_dofs = rDofSet | boost::adaptors::filtered(is_free_dof); + + auto free_dof_values = Vector{boost::size(free_dofs)}; + auto get_dof_value = [](const auto& rDof) { return rDof.GetSolutionStepValue(); }; + std::transform(std::begin(free_dofs), std::end(free_dofs), free_dof_values.begin(), get_dof_value); + + return norm_2(free_dof_values); + } + + std::vector mSubModelPartList; // List of every SubModelPart associated to an external load + std::vector mVariableNames; // Name of the nodal variable associated to every SubModelPart + }; // Class GeoMechanicsRammArcLengthStrategy } // namespace Kratos \ No newline at end of file diff --git a/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp b/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp index 03296529d279..6e21ce30c7c4 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/solving_strategy_factory.hpp @@ -69,7 +69,7 @@ class SolvingStrategyFactory ParametersUtilities::CopyOptionalParameters(rSolverSettings, strategy_entries); auto result = std::make_unique>( - rModelPart, scheme, solver, criteria, builder_and_solver, strategy_parameters, + rModelPart, scheme, criteria, builder_and_solver, strategy_parameters, max_iterations, calculate_reactions, reform_dof_set_at_each_step, move_mesh_flag); result->SetEchoLevel(echo_level); return result; diff --git a/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp b/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp index dbc7c4350839..e0abc0639d35 100644 --- a/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp +++ b/applications/GeoMechanicsApplication/custom_workflows/dgeoflow.cpp @@ -138,8 +138,8 @@ KratosExecute::GeoMechanicsNewtonRaphsonErosionProcessStrategyType::Pointer Krat auto pSolvingStrategy = Kratos::make_unique>( - rModelPart, p_scheme, p_solver, p_criteria, p_builder_and_solver, p_parameters, - MaxIterations, CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag); + rModelPart, p_scheme, p_criteria, p_builder_and_solver, p_parameters, MaxIterations, + CalculateReactions, ReformDofSetAtEachStep, MoveMeshFlag); pSolvingStrategy->Check(); return pSolvingStrategy; diff --git a/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py b/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py index c5580cb9f167..4a0d43e24856 100644 --- a/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py +++ b/applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py @@ -489,7 +489,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("loads_variable_list",self.settings["loads_variable_list"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsNewtonRaphsonStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, @@ -504,7 +503,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("max_piping_iterations", self.settings["max_piping_iterations"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsNewtonRaphsonErosionProcessStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params, @@ -545,7 +543,6 @@ def _ConstructSolver(self, builder_and_solver, strategy_type): self.strategy_params.AddValue("loads_variable_list",self.settings["loads_variable_list"]) solving_strategy = GeoMechanicsApplication.GeoMechanicsRammArcLengthStrategy(self.computing_model_part, self.scheme, - self.linear_solver, self.convergence_criterion, builder_and_solver, self.strategy_params,