[GeoMechanicsApplication] Fix ordering issue for reset displacement

applications/GeoMechanicsApplication/custom_elements/compressibility_calculator.cpp

@@ -38,7 +38,7 @@ Vector CompressibilityCalculator::RHSContribution(const Matrix& rCompressibility
 
 Matrix CompressibilityCalculator::LHSContribution(const Matrix& rCompressibilityMatrix) const
 {
-    return mInputProvider.GetMatrixScalarFactor()* rCompressibilityMatrix;
+    return mInputProvider.GetMatrixScalarFactor() * rCompressibilityMatrix;
 }
 
 std::pair<Matrix, Vector> CompressibilityCalculator::LocalSystemContribution()

applications/GeoMechanicsApplication/custom_elements/permeability_calculator.cpp

@@ -44,11 +44,11 @@ Vector PermeabilityCalculator::RHSContribution(const Matrix& rPermeabilityMatrix
 Matrix PermeabilityCalculator::CalculatePermeabilityMatrix() const
 {
     RetentionLaw::Parameters retention_parameters(mInputProvider.GetElementProperties());
-    const auto&              r_properties = mInputProvider.GetElementProperties();
-    auto integration_coefficients = mInputProvider.GetIntegrationCoefficients();
-    const auto   shape_function_gradients = mInputProvider.GetShapeFunctionGradients();
-    const auto   local_dimension          = shape_function_gradients[0].size2();
-    const Matrix constitutive_matrix =
+    const auto&              r_properties             = mInputProvider.GetElementProperties();
+    auto                     integration_coefficients = mInputProvider.GetIntegrationCoefficients();
+    const auto               shape_function_gradients = mInputProvider.GetShapeFunctionGradients();
+    const auto               local_dimension          = shape_function_gradients[0].size2();
+    const Matrix             constitutive_matrix =
         GeoElementUtilities::FillPermeabilityMatrix(r_properties, local_dimension);
 
     const auto   number_of_nodes           = shape_function_gradients[0].size1();

applications/GeoMechanicsApplication/custom_elements/permeability_calculator.h

@@ -26,7 +26,7 @@ class PermeabilityCalculator : public ContributionCalculator
     class InputProvider
     {
     public:
-        InputProvider(std::function<const Properties&()>                         GetElementProperties,
+        InputProvider(std::function<const Properties&()> GetElementProperties,
                       std::function<const std::vector<RetentionLaw::Pointer>&()> GetRetentionLaws,
                       std::function<Vector()>                        GetIntegrationCoefficients,
                       std::function<Vector(const Variable<double>&)> GetNodalValuesOf,
@@ -39,9 +39,9 @@ class PermeabilityCalculator : public ContributionCalculator
         {
         }
 
-        [[nodiscard]] const Properties& GetElementProperties() const;
+        [[nodiscard]] const Properties&                         GetElementProperties() const;
         [[nodiscard]] const std::vector<RetentionLaw::Pointer>& GetRetentionLaws() const;
-        [[nodiscard]] Vector GetIntegrationCoefficients() const;
+        [[nodiscard]] Vector                                    GetIntegrationCoefficients() const;
         [[nodiscard]] Vector GetNodalValues(const Variable<double>& rVariable) const;
         [[nodiscard]] Geometry<Node>::ShapeFunctionsGradientsType GetShapeFunctionGradients() const;
 

applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.cpp

@@ -15,15 +15,17 @@
 
 // Project includes
 #include "custom_python/add_custom_utilities_to_python.h"
-#include "includes/kratos_parameters.h"
 
-#include "custom_utilities/condition_utilities.hpp"
-#include "custom_utilities/element_utilities.hpp"
-#include "custom_utilities/interface_element_utilities.h"
+#include "custom_utilities/node_utilities.h"
 
 namespace Kratos::Python
 {
 
-void AddCustomUtilitiesToPython(pybind11::module&) { namespace py = pybind11; }
+void AddCustomUtilitiesToPython(const pybind11::module& rModule)
+{
+    pybind11::class_<NodeUtilities>(rModule, "NodeUtilities")
+        .def("AssignUpdatedVectorVariableToNonFixedComponentsOfNodes",
+             &NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes);
+}
 
 } // Namespace Kratos::Python.

applications/GeoMechanicsApplication/custom_python/add_custom_utilities_to_python.h

@@ -25,7 +25,7 @@
 namespace Kratos::Python
 {
 
-void AddCustomUtilitiesToPython(pybind11::module&);
+void AddCustomUtilitiesToPython(const pybind11::module& rModule);
 
 } // namespace Kratos::Python.
 

applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp

@@ -22,7 +22,8 @@ namespace Kratos
 
 void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode,
                                                                     const Variable<array_1d<double, 3>>& rDestinationVariable,
-                                                                    const array_1d<double, 3>& rNewValues)
+                                                                    const array_1d<double, 3>& rNewValues,
+                                                                    IndexType SolutionStepIndex)
 {
     const std::vector<std::string> components = {"X", "Y", "Z"};
     for (const auto& zipped : boost::combine(rNewValues, components)) {
@@ -33,9 +34,20 @@ void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode,
         if (const auto& component_variable =
                 VariablesUtilities::GetComponentFromVectorVariable(rDestinationVariable.Name(), component);
             !rNode.IsFixed(component_variable)) {
-            rNode.FastGetSolutionStepValue(component_variable, 0) = new_value;
+            rNode.FastGetSolutionStepValue(component_variable, SolutionStepIndex) = new_value;
         }
     }
 }
 
+void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(
+    const ModelPart::NodesContainerType& rNodes,
+    const Variable<array_1d<double, 3>>& rDestinationVariable,
+    const array_1d<double, 3>&           rNewValues,
+    IndexType                            SolutionStepIndex)
+{
+    block_for_each(rNodes, [&rDestinationVariable, &rNewValues, SolutionStepIndex](Node& rNode) {
+        AssignUpdatedVectorVariableToNonFixedComponents(rNode, rDestinationVariable, rNewValues, SolutionStepIndex);
+    });
+}
+
 } // namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/node_utilities.h

@@ -14,6 +14,7 @@
 
 #include "containers/array_1d.h"
 #include "containers/variable.h"
+#include "includes/model_part.h"
 
 namespace Kratos
 {
@@ -23,9 +24,17 @@ class Node;
 class KRATOS_API(GEO_MECHANICS_APPLICATION) NodeUtilities
 {
 public:
+    KRATOS_CLASS_POINTER_DEFINITION(NodeUtilities);
+
     static void AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode,
                                                                 const Variable<array_1d<double, 3>>& rDestinationVariable,
-                                                                const array_1d<double, 3>& rNewValues);
+                                                                const array_1d<double, 3>& rNewValues,
+                                                                IndexType SolutionStepIndex = 0);
+
+    static void AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(const ModelPart::NodesContainerType& rNodes,
+                                                                       const Variable<array_1d<double, 3>>& rDestinationVariable,
+                                                                       const array_1d<double, 3>& rNewValues,
+                                                                       IndexType SolutionStepIndex = 0);
 };
 
 } // namespace Kratos

applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.cpp

@@ -191,14 +191,6 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path&            rWorki
         std::vector<std::shared_ptr<Process>> processes = GetProcesses(project_parameters);
         std::vector<std::weak_ptr<Process>> process_observables(processes.begin(), processes.end());
 
-        if (mpTimeLoopExecutor) {
-            mpTimeLoopExecutor->SetCancelDelegate(rShouldCancel);
-            mpTimeLoopExecutor->SetProgressDelegate(rProgressDelegate);
-            mpTimeLoopExecutor->SetProcessObservables(process_observables);
-            mpTimeLoopExecutor->SetTimeIncrementor(MakeTimeIncrementor(project_parameters));
-            mpTimeLoopExecutor->SetSolverStrategyWrapper(MakeStrategyWrapper(project_parameters, rWorkingDirectory));
-        }
-
         for (const auto& process : processes) {
             process->ExecuteInitialize();
         }
@@ -208,6 +200,12 @@ int KratosGeoSettlement::RunStage(const std::filesystem::path&            rWorki
         }
 
         if (mpTimeLoopExecutor) {
+            mpTimeLoopExecutor->SetCancelDelegate(rShouldCancel);
+            mpTimeLoopExecutor->SetProgressDelegate(rProgressDelegate);
+            mpTimeLoopExecutor->SetProcessObservables(process_observables);
+            mpTimeLoopExecutor->SetTimeIncrementor(MakeTimeIncrementor(project_parameters));
+            mpTimeLoopExecutor->SetSolverStrategyWrapper(MakeStrategyWrapper(project_parameters, rWorkingDirectory));
+
             // For now, pass a dummy state. THIS PROBABLY NEEDS TO BE REFINED AT SOME POINT!
             TimeStepEndState start_of_loop_state;
             start_of_loop_state.convergence_state = TimeStepEndState::ConvergenceState::converged;
@@ -341,10 +339,8 @@ std::shared_ptr<StrategyWrapper> KratosGeoSettlement::MakeStrategyWrapper(const
     GetMainModelPart().CloneTimeStep();
 
     if (rProjectParameters["solver_settings"]["reset_displacements"].GetBool()) {
-        constexpr auto source_index      = std::size_t{0};
-        constexpr auto destination_index = std::size_t{1};
-        RestoreValuesOfNodalVariable(DISPLACEMENT, source_index, destination_index);
-        RestoreValuesOfNodalVariable(ROTATION, source_index, destination_index);
+        ResetValuesOfNodalVariable(DISPLACEMENT);
+        ResetValuesOfNodalVariable(ROTATION);
 
         VariableUtils{}.UpdateCurrentToInitialConfiguration(GetComputationalModelPart().Nodes());
     }

applications/GeoMechanicsApplication/custom_workflows/dgeosettlement.h

@@ -19,6 +19,7 @@
 #include "includes/kernel.h"
 #include "includes/kratos_export_api.h"
 
+#include "custom_utilities/node_utilities.h"
 #include "custom_utilities/process_factory.hpp"
 #include "geo_mechanics_application.h"
 #include "linear_solvers_application.h"
@@ -73,17 +74,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoSettlement
                             const std::stringstream&                rKratosLogBuffer) const;
 
     template <typename TVariableType>
-    void RestoreValuesOfNodalVariable(const TVariableType& rVariable, Node::IndexType SourceIndex, Node::IndexType DestinationIndex)
+    void ResetValuesOfNodalVariable(const TVariableType& rVariable)
     {
         if (!GetComputationalModelPart().HasNodalSolutionStepVariable(rVariable)) return;
 
-        VariableUtils{}.SetHistoricalVariableToZero(rVariable, GetComputationalModelPart().Nodes());
-
-        block_for_each(GetComputationalModelPart().Nodes(),
-                       [&rVariable, SourceIndex, DestinationIndex](auto& node) {
-            node.GetSolutionStepValue(rVariable, DestinationIndex) =
-                node.GetSolutionStepValue(rVariable, SourceIndex);
-        });
+        NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(
+            GetComputationalModelPart().Nodes(), rVariable, rVariable.Zero(), 0);
+        NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(
+            GetComputationalModelPart().Nodes(), rVariable, rVariable.Zero(), 1);
     }
 
     template <typename ProcessType>

applications/GeoMechanicsApplication/python_scripts/geomechanics_analysis.py

@@ -57,14 +57,15 @@ def _CalculateTotalDisplacement(self,node, old_total_displacement):
 
     def ResetIfHasNodalSolutionStepVariable(self, variable):
         if self._GetSolver().main_model_part.HasNodalSolutionStepVariable(variable):
-            KratosMultiphysics.VariableUtils().SetHistoricalVariableToZero(variable, self._GetSolver().GetComputingModelPart().Nodes)
-            for node in self._GetSolver().GetComputingModelPart().Nodes:
-                new_value = node.GetSolutionStepValue(variable, 0)
-                node.SetSolutionStepValue(variable, 1, new_value)
-
-    def ModifyInitialGeometry(self):
-        # Overrides the base class. Necessary to let reset_displacements function correctly i.c.w. prescribed displacements/rotations.
-        # The reset needs to take place befor the Initialize of the processes, as these will set the Dirichlet condition.
+            zero_vector = Kratos.Array3([0.0, 0.0, 0.0])
+            KratosGeo.NodeUtilities.AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(
+                self._GetSolver().GetComputingModelPart().Nodes, variable, zero_vector, 0)
+            KratosGeo.NodeUtilities.AssignUpdatedVectorVariableToNonFixedComponentsOfNodes(
+                self._GetSolver().GetComputingModelPart().Nodes, variable, zero_vector, 1)
+
+    def Initialize(self):
+        super().Initialize()
+
         self._GetSolver().main_model_part.ProcessInfo[KratosGeo.RESET_DISPLACEMENTS] = self.reset_displacements
         if self.reset_displacements:
             self.ResetIfHasNodalSolutionStepVariable(KratosMultiphysics.DISPLACEMENT)

applications/GeoMechanicsApplication/tests/truss_with_reset_displacement/ProjectParameters_stage2.json

@@ -30,7 +30,7 @@
         "block_builder":                      true,
         "solution_type":                      "Quasi-Static",
         "scheme_type":                        "Newmark",
-        "reset_displacements":                false,
+        "reset_displacements":                true,
         "newmark_beta":                       0.25,
         "newmark_gamma":                      0.5,
         "newmark_theta":                      0.5,