[Core] Add PMultigridBuilderAndSolver

applications/ContactStructuralMechanicsApplication/python_scripts/contact_structural_mechanics_implicit_dynamic_solver.py

@@ -137,16 +137,20 @@ def _CreateLinearSolver(self):
     def _CreateBuilderAndSolver(self):
         if self.contact_settings["mortar_type"].GetString() != "":
             linear_solver = self._GetLinearSolver()
-            if self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool():
+            builder_and_solver_type: str = self.settings["builder_and_solver_settings"]["type"].GetString()
+
+            if builder_and_solver_type == "block":
                 builder_and_solver = CSMA.ContactResidualBasedBlockBuilderAndSolver(linear_solver)
-            else:
+            elif builder_and_solver_type == "elimination":
                     # We use the elimination builder and solver
                     if self.settings["multi_point_constraints_used"].GetBool():
                         if (self.GetComputingModelPart().NumberOfMasterSlaveConstraints() > 0):
                             self.GetComputingModelPart().Set(KM.TO_SPLIT) # We set the flag for some operations
                         builder_and_solver = CSMA.ContactResidualBasedEliminationBuilderAndSolverWithConstraints(linear_solver)
                     else:
                         builder_and_solver = CSMA.ContactResidualBasedEliminationBuilderAndSolver(linear_solver)
+            else:
+                raise ValueError(f"Unsupported builder and solver type: {builder_and_solver_type}")
         else:
             builder_and_solver = super()._CreateBuilderAndSolver()
 

applications/ContactStructuralMechanicsApplication/python_scripts/contact_structural_mechanics_static_solver.py

@@ -138,7 +138,9 @@ def _CreateLinearSolver(self):
     def _CreateBuilderAndSolver(self):
         if self.contact_settings["mortar_type"].GetString() != "":
             linear_solver = self._GetLinearSolver()
-            if self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool():
+            builder_and_solver_type: str = self.settings["builder_and_solver_settings"]["type"].GetString()
+
+            if builder_and_solver_type == "block":
                 builder_and_solver = CSMA.ContactResidualBasedBlockBuilderAndSolver(linear_solver)
             else:
                     # We use the elimination builder and solver

applications/ContactStructuralMechanicsApplication/tests/ALM_frictionless_contact_test_2D/hyper_simple_patch_test_with_elimination_parameters.json

@@ -28,7 +28,7 @@
             "time_step" : 1.1
         },
         "builder_and_solver_settings" : {
-            "use_block_builder" : false
+            "type" : "elimination"
         },
         "clear_storage"                      : true,
         "compute_reactions"                  : false,

applications/ContactStructuralMechanicsApplication/tests/ALM_frictionless_contact_test_2D/hyper_simple_patch_test_with_elimination_with_constraints_parameters.json

@@ -28,7 +28,7 @@
             "time_step" : 1.1
         },
         "builder_and_solver_settings" : {
-            "use_block_builder" : false
+            "type" : "elimination"
         },
         "clear_storage"                      : true,
         "compute_reactions"                  : false,

applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_custom_scipy_base_solver.py

@@ -60,7 +60,8 @@ def _CreateLinearSolver(self):
         return KratosMultiphysics.LinearSolver()
 
     def _CreateSolutionStrategy(self):
-        if self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool():
+        builder_and_solver_type: str = self.settings["builder_and_solver_settings"]["type"].GetString()
+        if builder_and_solver_type == "block" or builder_and_solver_type == "block_lagrange":
             warn_msg = "In case an eigenvalue problem is computed an elimantion builder shall be used to ensure boundary conditions are applied correctly!"
             KratosMultiphysics.Logger.PrintWarning("CustomScipyBaseSolver", warn_msg)
 

applications/StructuralMechanicsApplication/python_scripts/structural_mechanics_solver.py

@@ -43,39 +43,7 @@ class MechanicalSolver(PythonSolver):
     settings -- Kratos parameters containing solver settings.
     """
     def __init__(self, model, custom_settings):
-        old_unused_settings = [
-            "use_computing_model_part",
-            "computing_model_part_name",
-            "problem_domain_sub_model_part_list",
-            "processes_sub_model_part_list"
-        ]
-
-        for old_setting in old_unused_settings:
-            if custom_settings.Has(old_setting):
-                KratosMultiphysics.Logger.PrintWarning("::[MechanicalSolver]:: ", 'Settings contain no longer used setting, please remove it: "{}"'.format(old_setting))
-                custom_settings.RemoveValue(old_setting)
-
-
-        settings_have_use_block_builder = custom_settings.Has("block_builder")
-        if settings_have_use_block_builder:
-            kratos_utilities.IssueDeprecationWarning('MechanicalSolver', 'Using "block_builder", please move it to "builder_and_solver_settings" as "use_block_builder"')
-            if not custom_settings.Has("builder_and_solver_settings"):
-                custom_settings.AddEmptyValue("builder_and_solver_settings")
-
-            custom_settings["builder_and_solver_settings"].AddValue("use_block_builder", custom_settings["block_builder"])
-            custom_settings.RemoveValue("block_builder")
-
-        settings_have_line_search = custom_settings.Has("line_search")
-        if settings_have_line_search:
-            kratos_utilities.IssueDeprecationWarning('MechanicalSolver', 'Using "line_search", please move it to "solving_strategy_settings" as "type"')
-            if custom_settings["line_search"].GetBool():
-                if not custom_settings.Has("solving_strategy_settings"):
-                    custom_settings.AddEmptyValue("solving_strategy_settings")
-
-                custom_settings["solving_strategy_settings"].AddEmptyValue("type")
-                custom_settings["solving_strategy_settings"]["type"].SetString("line_search")
-            custom_settings.RemoveValue("line_search")
-
+        self.__ConvertDeprecatedSettings(custom_settings)
         self._validate_settings_in_baseclass=True # To be removed eventually
         super().__init__(model, custom_settings)
 
@@ -140,8 +108,7 @@ def GetDefaultParameters(cls):
                 "advanced_settings" : { }
             },
             "builder_and_solver_settings" : {
-                "use_block_builder" : true,
-                "use_lagrange_BS"   : false,
+                "type" : "block",
                 "advanced_settings" : { }
             },
             "clear_storage": false,
@@ -379,19 +346,21 @@ def _GetLinearSolver(self):
         return self._linear_solver
 
     def _GetBuilderAndSolver(self):
+        builder_and_solver_type: str = self.settings["builder_and_solver_settings"]["type"].GetString()
         if not hasattr(self, '_builder_and_solver'):
             self._builder_and_solver = self._CreateBuilderAndSolver()
-        elif not self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool(): # Block builder and solver are unified with MPC and without. In the case of the elimination this could be a problem
+        elif builder_and_solver_type == "elimination": # Block builder and solver are unified with MPC and without. In the case of the elimination this could be a problem
             if self.GetComputingModelPart().NumberOfMasterSlaveConstraints() > 0 and not self.mpc_block_builder_initialized:
                 self.settings["multi_point_constraints_used"].SetBool(True)
                 self._builder_and_solver = self._CreateBuilderAndSolver()
                 self.mpc_block_builder_initialized = True
         return self._builder_and_solver
 
     def _GetSolutionStrategy(self):
+        builder_and_solver_type: str = self.settings["builder_and_solver_settings"]["type"].GetString()
         if not hasattr(self, '_mechanical_solution_strategy'):
             self._mechanical_solution_strategy = self._CreateSolutionStrategy()
-        elif not self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool(): # Block builder and solver are unified with MPC and without. In the case of the elimination this could be a problem
+        elif builder_and_solver_type == "elimination": # Block builder and solver are unified with MPC and without. In the case of the elimination this could be a problem
             if self.GetComputingModelPart().NumberOfMasterSlaveConstraints() > 0 and not self.mpc_block_builder_initialized:
                 self._mechanical_solution_strategy = self._CreateSolutionStrategy()
         return self._mechanical_solution_strategy
@@ -487,20 +456,33 @@ def _CreateLinearSolver(self):
             KratosMultiphysics.Logger.PrintInfo('::[MechanicalSolver]:: No linear solver was specified, using fastest available solver')
             return linear_solver_factory.CreateFastestAvailableDirectLinearSolver()
 
-    def _CreateBuilderAndSolver(self):
-        linear_solver = self._GetLinearSolver()
-        if self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool():
+    def _CreateBuilderAndSolver(self) -> KratosMultiphysics.BuilderAndSolver:
+        if not self.settings["builder_and_solver_settings"].Has("type"):
+            return self.__DeprecatedBuilderAndSolverFactory()
+
+        linear_solver: KratosMultiphysics.LinearSolver = self._GetLinearSolver()
+        builder_and_solver_name: str = self.settings["builder_and_solver_settings"]["type"].GetString()
+
+        if builder_and_solver_name == "block":
             bs_params = self.settings["builder_and_solver_settings"]["advanced_settings"]
-            if not self.settings["builder_and_solver_settings"]["use_lagrange_BS"].GetBool():
-                builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(linear_solver, bs_params)
-            else:
-                builder_and_solver = KratosMultiphysics.ResidualBasedBlockBuilderAndSolverWithLagrangeMultiplier(linear_solver, bs_params)
-        else:
+            return KratosMultiphysics.ResidualBasedBlockBuilderAndSolver(linear_solver, bs_params)
+        elif builder_and_solver_name == "block_lagrange":
+            return KratosMultiphysics.ResidualBasedBlockBuilderAndSolverWithLagrangeMultiplier(linear_solver, bs_params)
+        elif builder_and_solver_name == "elimination":
             if self.settings["multi_point_constraints_used"].GetBool():
-                builder_and_solver = KratosMultiphysics.ResidualBasedEliminationBuilderAndSolverWithConstraints(linear_solver)
+                return KratosMultiphysics.ResidualBasedEliminationBuilderAndSolverWithConstraints(linear_solver)
             else:
-                builder_and_solver = KratosMultiphysics.ResidualBasedEliminationBuilderAndSolver(linear_solver)
-        return builder_and_solver
+                return KratosMultiphysics.ResidualBasedEliminationBuilderAndSolver(linear_solver)
+        elif builder_and_solver_name == "p_multigrid":
+            parameters = self.settings["builder_and_solver_settings"]["advanced_settings"]
+            return KratosMultiphysics.PMultigridBuilderAndSolver(linear_solver, parameters)
+        else:
+            message: str = f"Invalid type for builder and solver '{builder_and_solver_name}'. Options are:\n"
+            message += "'block'\n"
+            message += "'block_lagrange'\n"
+            message += "'elimination'"
+            message += "'p_multigrid'"
+            raise ValueError(message)
 
     def _CreateScheme(self):
         """Create the solution scheme for the structural problem.
@@ -584,3 +566,59 @@ def _create_arc_length_strategy(self):
                                                                 self._GetBuilderAndSolver(),
                                                                 settings)
         return solving_strategy
+
+    def __ConvertDeprecatedSettings(self, settings: KratosMultiphysics.Parameters) -> None:
+        old_unused_settings = [
+            "use_computing_model_part",
+            "computing_model_part_name",
+            "problem_domain_sub_model_part_list",
+            "processes_sub_model_part_list"
+        ]
+
+        for old_setting in old_unused_settings:
+            if settings.Has(old_setting):
+                KratosMultiphysics.Logger.PrintWarning("::[MechanicalSolver]:: ", 'Settings contain no longer used setting, please remove it: "{}"'.format(old_setting))
+                settings.RemoveValue(old_setting)
+
+        settings_have_use_block_builder = settings.Has("block_builder")
+        if settings_have_use_block_builder:
+            kratos_utilities.IssueDeprecationWarning('MechanicalSolver', 'Using "block_builder", please move it to "builder_and_solver_settings" as "use_block_builder"')
+            if not settings.Has("builder_and_solver_settings"):
+                settings.AddEmptyValue("builder_and_solver_settings")
+
+            settings["builder_and_solver_settings"].AddValue("use_block_builder", settings["block_builder"])
+            settings.RemoveValue("block_builder")
+
+        settings_have_line_search = settings.Has("line_search")
+        if settings_have_line_search:
+            kratos_utilities.IssueDeprecationWarning('MechanicalSolver', 'Using "line_search", please move it to "solving_strategy_settings" as "type"')
+            if settings["line_search"].GetBool():
+                if not settings.Has("solving_strategy_settings"):
+                    settings.AddEmptyValue("solving_strategy_settings")
+
+                settings["solving_strategy_settings"].AddEmptyValue("type")
+                settings["solving_strategy_settings"]["type"].SetString("line_search")
+            settings.RemoveValue("line_search")
+
+        if settings.Has("builder_and_solver_settings") and not settings["builder_and_solver_settings"].Has("type"):
+            kratos_utilities.IssueDeprecationWarning(
+                "MechanicalSolver",
+                "Using deprecated builder and solver settings. Provide 'type' and 'advanced_settings' in the new system.")
+            bs_settings = settings["builder_and_solver_settings"]
+            updated_bs_settings = KratosMultiphysics.Parameters("""{}""")
+            if bs_settings.Has("use_block_builder"):
+                if bs_settings["use_block_builder"].GetBool():
+                    if bs_settings.Has("use_lagrange_BS") and bs_settings["use_lagrange_BS"].GetBool():
+                        updated_bs_settings.AddString("type", "block_lagrange")
+                    else:
+                        updated_bs_settings.AddString("type", "block")
+                else:
+                    updated_bs_settings.AddString("type", "elimination")
+            else:
+                updated_bs_settings.AddString("type", "block")
+
+            if bs_settings.Has("advanced_settings"):
+                updated_bs_settings.AddValue("advanced_settings", bs_settings["advanced_settings"])
+
+            settings.RemoveValue("builder_and_solver_settings")
+            settings.AddValue("builder_and_solver_settings", updated_bs_settings)

applications/StructuralMechanicsApplication/python_scripts/trilinos_structural_mechanics_solver.py

@@ -92,14 +92,19 @@ def _CreateBuilderAndSolver(self):
             guess_row_size = 15
         else:
             guess_row_size = 45
-        if self.settings["builder_and_solver_settings"]["use_block_builder"].GetBool():
+
+        builder_and_solver_type: str = self.settings["builder_and_solver_settings"]["type"].GetString()
+        if builder_and_solver_type == "block":
             builder_and_solver = TrilinosApplication.TrilinosBlockBuilderAndSolver(epetra_communicator,
                                                                                    guess_row_size,
                                                                                    linear_solver)
-        else:
+        elif builder_and_solver_type == "elimination":
             builder_and_solver = TrilinosApplication.TrilinosEliminationBuilderAndSolver(epetra_communicator,
                                                                                          guess_row_size,
                                                                                          linear_solver)
+        else:
+            raise ValueError(f"Unsupported builder and solver type '{builder_and_solver_type}'. Options are 'block', 'elimination'.")
+
         return builder_and_solver
 
     def _create_linear_strategy(self):

applications/StructuralMechanicsApplication/tests/rigid_test/test_elimination_mpc_parameters.json

@@ -22,7 +22,9 @@
         "time_stepping"                      : {
             "time_step" : 1.1
         },
-        "block_builder"                      : false,
+        "builder_and_solver_settings" : {
+            "type" : "elimination"
+        },
         "line_search"                        : false,
         "convergence_criterion"              : "residual_criterion",
         "displacement_relative_tolerance"    : 0.0001,

kratos/containers/nodal_data.h

@@ -186,7 +186,7 @@ class KRATOS_API(KRATOS_CORE) NodalData final
     ///@{
 
     /// Copy constructor.
-    NodalData(NodalData const& rOther);
+    NodalData(NodalData const& rOther) = default;
 
     ///@}
 

kratos/includes/dof.h

@@ -184,10 +184,18 @@ class Dof{
     {
     }
 
+    Dof(Dof&&) noexcept = default;
+
+    Dof(const Dof&) noexcept = default;
+
     ///@}
     ///@name Operators
     ///@{
 
+    Dof& operator=(Dof&&) noexcept = default;
+
+    Dof& operator=(const Dof&) noexcept = default;
+
     template<class TVariableType>
     typename TVariableType::Type& operator()(const TVariableType& rThisVariable, IndexType SolutionStepIndex = 0)
     {

kratos/includes/kratos_application.h

@@ -27,6 +27,7 @@
 #include "input_output/logger.h"
 #include "utilities/quaternion.h"
 #include "constraints/linear_master_slave_constraint.h"
+#include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // LinearMultifreedomConstraint
 
 // Geometries definition
 #include "geometries/register_kratos_components_for_geometry.h"
@@ -489,6 +490,7 @@ class KRATOS_API(KRATOS_CORE) KratosApplication {
     // Master-Slave base constraint
     const MasterSlaveConstraint mMasterSlaveConstraint;
     const LinearMasterSlaveConstraint mLinearMasterSlaveConstraint;
+    const LinearMultifreedomConstraint mLinearMultifreedomConstraint;
 
     // Periodic Condition
     const PeriodicCondition mPeriodicCondition;

kratos/python/add_constraint_to_python.cpp

@@ -20,6 +20,7 @@
 #include "includes/master_slave_constraint.h"
 #include "constraints/linear_master_slave_constraint.h"
 #include "constraints/slip_constraint.h"
+#include "solving_strategies/builder_and_solvers/p_multigrid/linear_multifreedom_constraint.hpp" // LinearMultifreedomConstraint
 
 namespace Kratos::Python
 {
@@ -133,6 +134,22 @@ void AddConstraintToPython(pybind11::module &m)
         >())
     ;
 
+    pybind11::class_<LinearMultifreedomConstraint,
+                     LinearMultifreedomConstraint::Pointer,
+                     MasterSlaveConstraint>(m, "LinearMultifreedomConstraint")
+        .def(pybind11::init([](const IndexType Id,
+                               const LinearMultifreedomConstraint::DofPointerVectorType& rDofs,
+                               const std::vector<std::size_t>& rConstraintLabels,
+                               const LinearMultifreedomConstraint::MatrixType& rRelationMatrix,
+                               const LinearMultifreedomConstraint::VectorType& rConstraintGaps){
+                                    return std::make_shared<LinearMultifreedomConstraint>(
+                                        Id,
+                                        LinearMultifreedomConstraint::DofPointerVectorType(rDofs),
+                                        rConstraintLabels,
+                                        rRelationMatrix,
+                                        rConstraintGaps);
+                               }))
+        ;
 
     pybind11::class_<SlipConstraint, SlipConstraint::Pointer, LinearMasterSlaveConstraint, Flags>
             (m, "SlipConstraint")