[MPM] Compute sum reactions for a given model part

applications/MPMApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -23,7 +23,8 @@
 #include "spaces/ublas_space.h"
 #include "linear_solvers/linear_solver.h"
 #include "custom_utilities/material_point_search_utility.h"
-#include "custom_utilities/material_point_generator_utility.cpp"
+#include "custom_utilities/material_point_generator_utility.h"
+#include "custom_utilities/reaction_utilities.h"
 
 
 namespace Kratos{
@@ -71,6 +72,16 @@ namespace Python{
         m.def("GenerateMaterialPointElement", GenerateMaterialPointElementAccordingToDimension);
         m.def("GenerateMaterialPointCondition", GenerateMaterialPointConditionAccordingToDimension);
         m.def("GenerateLagrangeNodes", MaterialPointGeneratorUtility::GenerateLagrangeNodes);
+
+        namespace py = pybind11;
+
+        // Calculate reaction forces
+        py::class_<ReactionUtilities> (m,"ReactionUtilities")
+            .def(py::init<>())
+            .def_static("CalculateGridConformingReaction", &ReactionUtilities::CalculateGridConformingReaction)
+            .def_static("CalculateNonConformingReaction", &ReactionUtilities::CalculateNonConformingReaction)
+            ;
+
     }
 
 }  // namespace Python.

applications/MPMApplication/custom_utilities/material_point_generator_utility.cpp

@@ -1039,19 +1039,19 @@ namespace Kratos::MaterialPointGeneratorUtility
     //
 
 
-    template void GenerateMaterialPointElement<2>(ModelPart& rBackgroundGridModelPart,
+    template KRATOS_API(MPM_APPLICATION) void GenerateMaterialPointElement<2>(ModelPart& rBackgroundGridModelPart,
                                         ModelPart& rInitialModelPart,
                                         ModelPart& rMPMModelPart,
                                         bool IsMixedFormulation);
-    template void GenerateMaterialPointElement<3>(ModelPart& rBackgroundGridModelPart,
+    template KRATOS_API(MPM_APPLICATION) void GenerateMaterialPointElement<3>(ModelPart& rBackgroundGridModelPart,
                                         ModelPart& rInitialModelPart,
                                         ModelPart& rMPMModelPart,
                                         bool IsMixedFormulation);
 
-    template void GenerateMaterialPointCondition<2>(ModelPart& rBackgroundGridModelPart,
+    template KRATOS_API(MPM_APPLICATION) void GenerateMaterialPointCondition<2>(ModelPart& rBackgroundGridModelPart,
                                             ModelPart& rInitialModelPart,
                                             ModelPart& rMPMModelPart);
-    template void GenerateMaterialPointCondition<3>(ModelPart& rBackgroundGridModelPart,
+    template KRATOS_API(MPM_APPLICATION) void GenerateMaterialPointCondition<3>(ModelPart& rBackgroundGridModelPart,
                                             ModelPart& rInitialModelPart,
                                             ModelPart& rMPMModelPart);
 

applications/MPMApplication/custom_utilities/material_point_generator_utility.h

@@ -65,7 +65,7 @@ namespace Kratos::MaterialPointGeneratorUtility
      * @details Generating material points using a designated shape functions
      */
     template<SizeType TDimension>
-    void GenerateMaterialPointElement(  ModelPart& rBackgroundGridModelPart,
+    void KRATOS_API(MPM_APPLICATION) GenerateMaterialPointElement(  ModelPart& rBackgroundGridModelPart,
                                         ModelPart& rInitialModelPart,
                                         ModelPart& rMPMModelPart,
                                         bool IsMixedFormulation=false);
@@ -74,15 +74,15 @@ namespace Kratos::MaterialPointGeneratorUtility
      * @details Generating material point condition using a designated shape functions
      */
     template<SizeType TDimension>
-    void GenerateMaterialPointCondition(ModelPart& rBackgroundGridModelPart,
+    void KRATOS_API(MPM_APPLICATION) GenerateMaterialPointCondition(ModelPart& rBackgroundGridModelPart,
                                             ModelPart& rInitialModelPart,
                                             ModelPart& rMPMModelPart);
 
     /**
      * @brief Function to add dofs to elements for Lagrange multiplier.
      * @details non-conforming Perturbed Lagrangian or Lagrange multiplier conditions
      */
-    void GenerateLagrangeNodes(ModelPart& rBackgroundGridModelPart);
+    void KRATOS_API(MPM_APPLICATION) GenerateLagrangeNodes(ModelPart& rBackgroundGridModelPart);
     /**
      * @brief Function to Initiate material point condition.
      * @details Generating material point condition using a designated shape functions

applications/MPMApplication/custom_utilities/reaction_utilities.cpp

@@ -0,0 +1,38 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   Multi-Physics
+//
+//  License:         BSD License
+//                   Kratos default license: kratos/license.txt
+//
+//  Main authors:    Nicolo Crescenzio
+
+// System includes
+
+// External includes
+
+// Project includes
+#include "utilities/variable_utils.h"
+
+// Application includes
+#include "mpm_application_variables.h"
+#include "custom_utilities/reaction_utilities.h"
+
+namespace Kratos
+{
+    array_1d<double,3> ReactionUtilities::CalculateGridConformingReaction(ModelPart& rModelPart) {
+        return VariableUtils().SumHistoricalVariable<array_1d<double,3>>(REACTION, rModelPart, 0);
+    }
+
+    array_1d<double,3> ReactionUtilities::CalculateNonConformingReaction(ModelPart& rModelPart) {
+        const ProcessInfo& r_current_process_info = rModelPart.GetProcessInfo();
+        std::vector<array_1d<double,3>> mpc_reaction{ ZeroVector(3) };
+        return block_for_each<SumReduction<array_1d<double,3>>>(rModelPart.Conditions(), mpc_reaction,
+            [&r_current_process_info](auto& condition, auto& rReaction) {
+                condition.CalculateOnIntegrationPoints(MPC_CONTACT_FORCE, rReaction, r_current_process_info);
+                return rReaction[0];
+            });
+    }
+}

applications/MPMApplication/custom_utilities/reaction_utilities.h

@@ -0,0 +1,50 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   Multi-Physics
+//
+//  License:         BSD License
+//                   Kratos default license: kratos/license.txt
+//
+//  Main author:     Nicolo Crescenzio
+
+#pragma once
+
+// System includes
+
+// External includes
+
+// Project includes
+#include "includes/define.h"
+#include "includes/model_part.h"
+
+namespace Kratos
+{
+
+// Auxiliary utility to compute the reaction over a line or surface
+// The line or surface can be a part of the boundary of the background grid
+// (conforming boundary conditions) or can be defined through material
+// point conditions (non-conforming boundary conditions)
+class KRATOS_API(MPM_APPLICATION) ReactionUtilities
+{
+public:
+
+    /**
+    * @brief Computes the sum of the REACTION variable for all nodes in a given model part
+    * @param rModelPart reference to the model part where the force is to be computed
+    * @return An array containing the force value
+    */
+    static array_1d<double, 3> CalculateGridConformingReaction(ModelPart &rModelPart);
+
+    /**
+    * @brief Computes the sum of the MPC_CONTACT_FORCE variable for all material point
+    * conditions in a given model part
+    * @param rModelPart reference to the model part where the force is to be computed
+    * @return An array containing the force value
+    */
+    static array_1d<double, 3> CalculateNonConformingReaction(ModelPart &rModelPart);
+
+}; // Class ReactionUtilities
+
+} // namespace Kratos

applications/MPMApplication/python_scripts/mpm_grid_conforming_reaction_output_process.py

@@ -0,0 +1,29 @@
+import KratosMultiphysics
+import KratosMultiphysics.MPMApplication as KratosMPM
+from KratosMultiphysics.MPMApplication.mpm_reaction_output_process import MPMReactionOutputProcess
+
+def Factory(settings: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> KratosMultiphysics.OutputProcess:
+    if not isinstance(model, KratosMultiphysics.Model):
+        raise Exception("expected input shall be a Model object")
+    if not isinstance(settings, KratosMultiphysics.Parameters):
+        raise Exception("expected input shall be a Parameters object, encapsulating a json string")
+    return MPMGridConformingReactionOutputProcess(model, settings["Parameters"])
+
+class MPMGridConformingReactionOutputProcess(MPMReactionOutputProcess):
+    """
+    The specific implementation for the output of reaction on grid conforming boundaries
+    """
+    def __init__(self, model: KratosMultiphysics.Model, settings: KratosMultiphysics.Parameters) -> None:
+        super().__init__(model, settings)
+        if not self.model_part_name.startswith('Background_Grid'):
+            err_msg  = f"'{self.__class__.__name__}' prints the sum of the reaction of a given SubModelPart of the background grid.\n"
+            err_msg += f"Model part name must start with 'Background_Grid'. Input model part name: '{self.model_part_name}'"
+            raise Exception(err_msg)
+
+    def _GetFileHeader(self):
+        header  = f'# Grid conforming reaction for model part {self.model_part_name}\n'
+        header +=  '# Time Fx Fy Fz\n'
+        return header
+
+    def _GetReaction(self):
+        return KratosMPM.ReactionUtilities.CalculateGridConformingReaction(self.model_part)

applications/MPMApplication/python_scripts/mpm_non_conforming_reaction_output_process.py

@@ -0,0 +1,28 @@
+import KratosMultiphysics
+import KratosMultiphysics.MPMApplication as KratosMPM
+from KratosMultiphysics.MPMApplication.mpm_reaction_output_process import MPMReactionOutputProcess
+
+def Factory(settings: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> KratosMultiphysics.OutputProcess:
+    if not isinstance(model, KratosMultiphysics.Model):
+        raise Exception("expected input shall be a Model object")
+    if not isinstance(settings, KratosMultiphysics.Parameters):
+        raise Exception("expected input shall be a Parameters object, encapsulating a json string")
+    return MPMNonConformingReactionOutputProcess(model, settings["Parameters"])
+
+class MPMNonConformingReactionOutputProcess(MPMReactionOutputProcess):
+    """
+    The specific implementation for the output of reaction on non-conforming boundaries
+    defined through material point conditions
+    """
+    def __init__(self, model: KratosMultiphysics.Model, settings: KratosMultiphysics.Parameters) -> None:
+        super().__init__(model, settings)
+        if self.model_part_name.startswith('Background_Grid'):
+            self.model_part_name = self.model_part_name.replace('Background_Grid','MPM_Material')
+
+    def _GetFileHeader(self):
+        header  = f'# Non-conforming reaction for model part {self.model_part_name}\n'
+        header +=  '# Time Fx Fy Fz\n'
+        return header
+
+    def _GetReaction(self):
+        return KratosMPM.ReactionUtilities.CalculateNonConformingReaction(self.model_part)

applications/MPMApplication/python_scripts/mpm_reaction_output_process.py

@@ -0,0 +1,83 @@
+import KratosMultiphysics
+from KratosMultiphysics.time_based_ascii_file_writer_utility import TimeBasedAsciiFileWriterUtility
+
+def Factory(settings: KratosMultiphysics.Parameters, model: KratosMultiphysics.Model) -> KratosMultiphysics.OutputProcess:
+    if not isinstance(model, KratosMultiphysics.Model):
+        raise Exception("expected input shall be a Model object")
+    if not isinstance(settings, KratosMultiphysics.Parameters):
+        raise Exception("expected input shall be a Parameters object, encapsulating a json string")
+    return MPMReactionOutputProcess(model, settings["Parameters"])
+
+class MPMReactionOutputProcess(KratosMultiphysics.OutputProcess):
+    """
+    Auxiliary base class to output total reaction of over a line or surface.
+    A derived class needs to be implemented to be able to use this functionality
+    (in both conforming and non-conforming boundary conditions),
+    as calling the base class alone is not enough.
+    """
+    def __init__(self, model: KratosMultiphysics.Model, settings: KratosMultiphysics.Parameters) -> None:
+        super().__init__()
+
+        self.model = model
+        self.params = settings
+        self.params.ValidateAndAssignDefaults(self.GetDefaultParameters())
+
+        self.output_file = None
+        self.controller = KratosMultiphysics.OutputController(self.model, self.params)
+        self.format = self.params["print_format"].GetString()
+
+        self.model_part_name = self.params["model_part_name"].GetString()
+        if not self.model_part_name:
+            raise Exception('No "model_part_name" was specified!')
+
+    @staticmethod
+    def GetDefaultParameters():
+        return KratosMultiphysics.Parameters("""
+            {
+                "model_part_name"           : "",
+                "output_interval"           : 1.0,
+                "output_control_type"       : "step",
+                "print_format"              : ".8f",
+                "output_file_settings"      : {}
+            }
+            """)
+
+    def ExecuteBeforeSolutionLoop(self) -> None:
+        self.model_part = self.model[self.model_part_name]
+        file_handler_params = KratosMultiphysics.Parameters(self.params["output_file_settings"])
+
+        # If output file name is not specified, use default value: <model_part_name>_reaction.dat
+        if not file_handler_params.Has("file_name"):
+            output_file_name = f"{self.model_part_name}_reaction"
+            wrn_msg = f"File name not specified, using the default name: {output_file_name}"
+            KratosMultiphysics.Logger.PrintWarning(self.__class__.__name__,wrn_msg)
+            file_handler_params.AddString("file_name",output_file_name)
+
+        self.output_file = TimeBasedAsciiFileWriterUtility(self.model_part, file_handler_params, self._GetFileHeader()).file
+        self.PrintOutput()
+
+    def PrintOutput(self) -> None:
+        if self.output_file:
+            time = self.model_part.ProcessInfo[KratosMultiphysics.TIME]
+            reaction_force = self._GetReaction() # Compute reaction force
+            out = f"{str(time)} {reaction_force[0]:{self.format}} {reaction_force[1]:{self.format}} {reaction_force[2]:{self.format}}\n"
+            self.output_file.write(out)
+            self.output_file.flush()
+            self.controller.Update()
+
+    def IsOutputStep(self) -> bool:
+        return self.controller.Evaluate()
+
+    def ExecuteFinalize(self) -> None:
+        if self.output_file:
+            self.output_file.close()
+
+    def _GetFileHeader(self):
+        err_msg  = '{self.__class__.__name__}: _GetFileHeader called in base class\n'
+        err_msg += 'this needs to be implemented and called from derived class'
+        raise Exception(err_msg)
+
+    def _GetReaction(self):
+        err_msg  = '{self.__class__.__name__}: _GetReaction called in base class\n'
+        err_msg += 'this needs to be implemented and called from derived class'
+        raise Exception(err_msg)