[GeoMechanicsApplication] Adopt structural truss element with reset displacement

applications/GeoMechanicsApplication/custom_processes/README.md

@@ -5,6 +5,7 @@ This folder contains the custom processes that are used in the GeoMechanicsAppli
 Documented processes:
 - $c-\phi$ reduction process
 - [GeoExtrapolateIntegrationPointValuesToNodesProcess](#extrapolation-of-integration-values-to-nodes)
+- [ResetDisplacementProcess](#reset-displacement-process)
 
 ## $c-\phi$ reduction process
 For the assessment of a safety factor to characterize slope stability, a Mohr-Coulomb material based $c-\phi$ reduction 
@@ -59,5 +60,15 @@ Where the `model_part_name` should contain the name of the model part where the
 
 When this process is added to the `ProjectParameters.json`, the variables specified in `list_of_variables` can be exported as nodal output (e.g. as `nodal_results` in the `GiDOutputProcess`). 
 
+## Reset displacement process
+The `ResetDisplacementProcess` can be used to change the reference point of the displacements to the displacement at the start of that stage.
+
+### Requirements
+For this process to work, the following requirements have to be met:
+1. The elements in the model part that the process is applied to should have an implementation for `CalculateOnIntegrationPoints` that calculates the PK2_STRESS_VECTOR as well as an overload of `CalculateOnIntegrationPoints` that returns a list of ConstitutiveLaw::Pointer objects for each integration point.
+2. The input type of the model can only be "rest" (short for restarted), to ensure that the state of the model is retained from the previous stage. The reason for this, is that the constitutive laws are used at the start of a state to calculate the initial stresses based on the history. If the model is not 'restarted', the constitutive laws will be cleared and the initial stresses can not be calculated correctly.
+3. The ConstitutiveLaw used in the elements this process is applied to should use the `InitialState` to apply the initial stresses to the calculated stresses.
+
+
 ## References
 <a id="1">[1]</a> Brinkgreve, R.B.J., Bakker, H.L., 1991. Non-linear finite element analysis of safety factors, Computer Methods and Advances in Geomechanics, Beer, Booker & Carterr (eds), Balkema, Rotterdam.

applications/GeoMechanicsApplication/custom_processes/reset_displacement_process.cpp

@@ -0,0 +1,75 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Richard Faasse
+//
+
+#include "reset_displacement_process.h"
+#include "includes/initial_state.h"
+#include "includes/kratos_parameters.h"
+#include "includes/model_part.h"
+#include "includes/ublas_interface.h"
+#include "includes/variables.h"
+
+#include <vector>
+
+namespace Kratos
+{
+
+ResetDisplacementProcess::ResetDisplacementProcess(ModelPart& rModelPart, const Parameters&)
+    : mrModelPart(rModelPart)
+{
+}
+
+void ResetDisplacementProcess::ExecuteInitialize()
+{
+    block_for_each(mrModelPart.Elements(), [this](Element& rElement) {
+        std::vector<Vector> stresses_on_integration_points;
+        rElement.CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stresses_on_integration_points,
+                                              mrModelPart.GetProcessInfo());
+        std::vector<ConstitutiveLaw::Pointer> constitutive_laws;
+        rElement.CalculateOnIntegrationPoints(CONSTITUTIVE_LAW, constitutive_laws, mrModelPart.GetProcessInfo());
+
+        CheckRetrievedElementData(constitutive_laws, stresses_on_integration_points, rElement.GetId());
+
+        for (auto i = std::size_t{0}; i < constitutive_laws.size(); ++i) {
+            auto p_initial_state = make_intrusive<InitialState>();
+            p_initial_state->SetInitialStressVector(stresses_on_integration_points[i]);
+            p_initial_state->SetInitialStrainVector(ZeroVector{constitutive_laws[i]->GetStrainSize()});
+            constitutive_laws[i]->SetInitialState(p_initial_state);
+        }
+    });
+}
+
+int ResetDisplacementProcess::Check()
+{
+    KRATOS_ERROR_IF_NOT(mrModelPart.GetProcessInfo()[IS_RESTARTED])
+        << "The IS_RESTARTED flag must be set to true in the ProcessInfo of the "
+           "model part. Please use the \"rest\" option for the model input type";
+
+    return 0;
+}
+
+void ResetDisplacementProcess::CheckRetrievedElementData(const std::vector<ConstitutiveLaw::Pointer>& rConstitutiveLaws,
+                                                         const std::vector<Vector>& rStressesOnIntegrationPoints,
+                                                         IndexType ElementId)
+{
+    KRATOS_ERROR_IF(rConstitutiveLaws.empty())
+        << "The constitutive laws on the integration points could not be retrieved for element "
+        << ElementId << std::endl;
+    KRATOS_ERROR_IF(rStressesOnIntegrationPoints.empty())
+        << "The stress vectors on the integration points could not be retrieved for element "
+        << ElementId << std::endl;
+    KRATOS_ERROR_IF(rStressesOnIntegrationPoints.size() != rConstitutiveLaws.size())
+        << "Number of retrieved stress vectors (" << rStressesOnIntegrationPoints.size()
+        << ") does not match the number of constitutive laws (" << rConstitutiveLaws.size()
+        << ") for element " << ElementId << std::endl;
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_processes/reset_displacement_process.h

@@ -0,0 +1,46 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Richard Faasse
+//
+
+#pragma once
+
+#include "includes/constitutive_law.h"
+#include "processes/process.h"
+
+namespace Kratos
+{
+
+class ModelPart;
+class Parameters;
+class Element;
+
+class KRATOS_API(GEO_MECHANICS_APPLICATION) ResetDisplacementProcess : public Process
+{
+public:
+    KRATOS_CLASS_POINTER_DEFINITION(ResetDisplacementProcess);
+    ResetDisplacementProcess(ModelPart& rModelPart, const Parameters&);
+    ~ResetDisplacementProcess() override = default;
+
+    ResetDisplacementProcess(const ResetDisplacementProcess&)            = delete;
+    ResetDisplacementProcess& operator=(const ResetDisplacementProcess&) = delete;
+
+    void ExecuteInitialize() override;
+    int  Check() override;
+
+private:
+    ModelPart& mrModelPart;
+
+    static void CheckRetrievedElementData(const std::vector<ConstitutiveLaw::Pointer>& rConstitutiveLaws,
+                                          const std::vector<Vector>& rStressesOnIntegrationPoints,
+                                          IndexType                  ElementId);
+};
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_python/add_custom_processes_to_python.cpp

@@ -46,6 +46,7 @@
 #include "custom_processes/deactivate_conditions_on_inactive_elements_process.hpp"
 #include "custom_processes/find_neighbour_elements_of_conditions_process.hpp"
 #include "custom_processes/geo_extrapolate_integration_point_values_to_nodes_process.h"
+#include "custom_processes/reset_displacement_process.h"
 #include "custom_processes/set_absorbing_boundary_parameters_process.hpp"
 #include "custom_processes/set_multiple_moving_loads.h"
 #include "custom_processes/set_parameter_field_process.hpp"
@@ -176,6 +177,10 @@ void AddCustomProcessesToPython(pybind11::module& m)
     py::class_<CalculateIncrementalDisplacementProcess, CalculateIncrementalDisplacementProcess::Pointer, Process>(
         m, "CalculateIncrementalDisplacementProcess")
         .def(py::init<ModelPart&, const Parameters&>());
+
+    py::class_<ResetDisplacementProcess, ResetDisplacementProcess::Pointer, Process>(
+        m, "ResetDisplacementProcess")
+        .def(py::init<ModelPart&, const Parameters&>());
 }
 
 } // Namespace Kratos::Python.

applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py

@@ -6,6 +6,8 @@
 import KratosMultiphysics.StructuralMechanicsApplication as StructuralMechanicsApplication
 import KratosMultiphysics.GeoMechanicsApplication as GeoMechanicsApplication
 
+import json
+
 def CreateSolver(model, custom_settings):
     return GeoMechanicalSolver(model, custom_settings)
 
@@ -205,6 +207,14 @@ def PrepareModelPart(self):
                                                   self.settings["nodal_smoothing"].GetBool())
 
         self.main_model_part.ProcessInfo.SetValue(KratosMultiphysics.STEP, 0)
+        self.computing_model_part_name = "porous_computational_model_part"
+
+        sub_model_part_names = [f"sub_{name.GetString()}" for name in self.settings["body_domain_sub_model_part_list"]]
+        self.body_domain_sub_sub_model_part_list = KratosMultiphysics.Parameters(json.dumps(sub_model_part_names))
+
+        sub_model_part_names = [f"sub_{name.GetString()}" for name in self.settings["loads_sub_model_part_list"]]
+        self.loads_sub_sub_model_part_list = KratosMultiphysics.Parameters(json.dumps(sub_model_part_names))
+
         if not self.main_model_part.ProcessInfo[KratosMultiphysics.IS_RESTARTED]:
             ## Executes the check and prepare model process (Create computing_model_part and set constitutive law)
             self._ExecuteCheckAndPrepare()
@@ -402,22 +412,6 @@ def _GetLinearSolver(self):
         return self.linear_solver
 
     def _ExecuteCheckAndPrepare(self):
-
-        self.computing_model_part_name = "porous_computational_model_part"
-
-        # Create list of sub sub model parts (it is a copy of the standard lists with a different name)
-        import json
-
-        self.body_domain_sub_sub_model_part_list = []
-        for i in range(self.settings["body_domain_sub_model_part_list"].size()):
-            self.body_domain_sub_sub_model_part_list.append("sub_"+self.settings["body_domain_sub_model_part_list"][i].GetString())
-        self.body_domain_sub_sub_model_part_list = KratosMultiphysics.Parameters(json.dumps(self.body_domain_sub_sub_model_part_list))
-
-        self.loads_sub_sub_model_part_list = []
-        for i in range(self.settings["loads_sub_model_part_list"].size()):
-            self.loads_sub_sub_model_part_list.append("sub_"+self.settings["loads_sub_model_part_list"][i].GetString())
-        self.loads_sub_sub_model_part_list = KratosMultiphysics.Parameters(json.dumps(self.loads_sub_sub_model_part_list))
-
         # Auxiliary parameters object for the CheckAndPepareModelProcess
         params = KratosMultiphysics.Parameters("{}")
         params.AddEmptyValue("computing_model_part_name").SetString(self.computing_model_part_name)

applications/GeoMechanicsApplication/python_scripts/reset_displacement_process.py

@@ -0,0 +1,10 @@
+import KratosMultiphysics as Core
+import KratosMultiphysics.GeoMechanicsApplication as Geo
+
+
+def Factory(settings, model):
+    if not isinstance(settings, Core.Parameters):
+        raise TypeError("expected input shall be a Parameters object, encapsulating a json string")
+
+    model_part = model[settings["Parameters"]["model_part_name"].GetString()]
+    return Geo.ResetDisplacementProcess(model_part, settings["Parameters"])