[GeoMechanicsApplication] Fix prescribed acceleration, velocity and d…

…t_water_pressure (#12208)

This commit fixes the issue that the time integration schemes would overwrite values for prescribed acceleration/velocity/dt_water_pressure, even if these variables were 'Fixed'

applications/GeoMechanicsApplication/custom_strategies/schemes/backward_euler_scheme.hpp

@@ -12,6 +12,7 @@
 
 #pragma once
 
+#include "custom_utilities/node_utilities.h"
 #include "geo_mechanics_application_variables.h"
 #include "geomechanics_time_integration_scheme.hpp"
 
@@ -52,13 +53,13 @@ class BackwardEulerScheme : public GeoMechanicsTimeIntegrationScheme<TSparseSpac
         KRATOS_TRY
 
         block_for_each(rModelPart.Nodes(), [this](Node& rNode) {
-            // For the Backward Euler schemes the first derivatives should be
-            // updated before calculating the second derivatives
             UpdateVectorTimeDerivatives(rNode);
 
             for (const auto& r_first_order_scalar_variable : this->GetFirstOrderScalarVariables()) {
-                SetDerivative(r_first_order_scalar_variable.first_time_derivative,
-                              r_first_order_scalar_variable.instance, rNode);
+                if (rNode.IsFixed(r_first_order_scalar_variable.first_time_derivative)) continue;
+
+                rNode.FastGetSolutionStepValue(r_first_order_scalar_variable.first_time_derivative) =
+                    CalculateDerivative(r_first_order_scalar_variable.instance, rNode);
             }
         });
 
@@ -70,23 +71,28 @@ class BackwardEulerScheme : public GeoMechanicsTimeIntegrationScheme<TSparseSpac
         for (const auto& r_second_order_vector_variable : this->GetSecondOrderVectorVariables()) {
             if (!rNode.SolutionStepsDataHas(r_second_order_vector_variable.instance)) continue;
 
-            SetDerivative(r_second_order_vector_variable.first_time_derivative,
-                          r_second_order_vector_variable.instance, rNode);
+            const auto updated_first_time_derivative =
+                CalculateDerivative(r_second_order_vector_variable.instance, rNode);
+
+            NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(
+                rNode, r_second_order_vector_variable.first_time_derivative, updated_first_time_derivative);
 
             // Make sure that setting the second_time_derivative is done
             // after setting the first_time_derivative.
-            SetDerivative(r_second_order_vector_variable.second_time_derivative,
-                          r_second_order_vector_variable.first_time_derivative, rNode);
+            const auto updated_second_time_derivative =
+                CalculateDerivative(r_second_order_vector_variable.first_time_derivative, rNode);
+
+            NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(
+                rNode, r_second_order_vector_variable.second_time_derivative, updated_second_time_derivative);
         }
     }
 
     template <class T>
-    void SetDerivative(const Variable<T>& derivative_variable, const Variable<T>& instance_variable, Node& rNode) const
+    T CalculateDerivative(const Variable<T>& rInstanceVariable, Node& rNode) const
     {
-        rNode.FastGetSolutionStepValue(derivative_variable) =
-            (rNode.FastGetSolutionStepValue(instance_variable, 0) -
-             rNode.FastGetSolutionStepValue(instance_variable, 1)) /
-            this->GetDeltaTime();
+        return (rNode.FastGetSolutionStepValue(rInstanceVariable, 0) -
+                rNode.FastGetSolutionStepValue(rInstanceVariable, 1)) /
+               this->GetDeltaTime();
     }
 };
 

applications/GeoMechanicsApplication/custom_strategies/schemes/generalized_newmark_scheme.hpp

@@ -11,6 +11,8 @@
 //
 #pragma once
 
+#include "custom_utilities/node_utilities.h"
+#include "custom_utilities/variables_utilities.hpp"
 #include "geomechanics_time_integration_scheme.hpp"
 #include "includes/model_part.h"
 #include <optional>
@@ -131,6 +133,8 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
 
     void UpdateScalarTimeDerivative(Node& rNode, const Variable<double>& variable, const Variable<double>& dt_variable) const
     {
+        if (rNode.IsFixed(dt_variable)) return;
+
         const auto delta_variable =
             rNode.FastGetSolutionStepValue(variable, 0) - rNode.FastGetSolutionStepValue(variable, 1);
         const auto previous_dt_variable = rNode.FastGetSolutionStepValue(dt_variable, 1);
@@ -145,12 +149,15 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
         for (const auto& r_second_order_vector_variable : this->GetSecondOrderVectorVariables()) {
             if (!rNode.SolutionStepsDataHas(r_second_order_vector_variable.instance)) continue;
 
-            noalias(rNode.FastGetSolutionStepValue(r_second_order_vector_variable.first_time_derivative, 0)) =
+            const auto updated_first_derivative =
                 rNode.FastGetSolutionStepValue(r_second_order_vector_variable.first_time_derivative, 1) +
                 (1.0 - GetGamma()) * this->GetDeltaTime() *
                     rNode.FastGetSolutionStepValue(r_second_order_vector_variable.second_time_derivative, 1) +
                 GetGamma() * this->GetDeltaTime() *
                     rNode.FastGetSolutionStepValue(r_second_order_vector_variable.second_time_derivative, 0);
+
+            NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(
+                rNode, r_second_order_vector_variable.first_time_derivative, updated_first_derivative);
         }
     }
 
@@ -159,14 +166,17 @@ class GeneralizedNewmarkScheme : public GeoMechanicsTimeIntegrationScheme<TSpars
         for (const auto& r_second_order_vector_variable : this->GetSecondOrderVectorVariables()) {
             if (!rNode.SolutionStepsDataHas(r_second_order_vector_variable.instance)) continue;
 
-            noalias(rNode.FastGetSolutionStepValue(r_second_order_vector_variable.second_time_derivative, 0)) =
+            const auto updated_second_time_derivative =
                 ((rNode.FastGetSolutionStepValue(r_second_order_vector_variable.instance, 0) -
                   rNode.FastGetSolutionStepValue(r_second_order_vector_variable.instance, 1)) -
                  this->GetDeltaTime() * rNode.FastGetSolutionStepValue(
                                             r_second_order_vector_variable.first_time_derivative, 1) -
                  (0.5 - GetBeta()) * this->GetDeltaTime() * this->GetDeltaTime() *
                      rNode.FastGetSolutionStepValue(r_second_order_vector_variable.second_time_derivative, 1)) /
                 (GetBeta() * this->GetDeltaTime() * this->GetDeltaTime());
+
+            NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(
+                rNode, r_second_order_vector_variable.second_time_derivative, updated_second_time_derivative);
         }
     }
 

applications/GeoMechanicsApplication/custom_strategies/schemes/geomechanics_time_integration_scheme.hpp

@@ -11,6 +11,7 @@
 //
 #pragma once
 
+#include "custom_utilities/variables_utilities.hpp"
 #include "geo_mechanics_application_variables.h"
 #include "solving_strategies/schemes/scheme.h"
 
@@ -315,12 +316,6 @@ class GeoMechanicsTimeIntegrationScheme : public Scheme<TSparseSpace, TDenseSpac
     }
 
 protected:
-    const Variable<double>& GetComponentFromVectorVariable(const Variable<array_1d<double, 3>>& rSource,
-                                                           const std::string& rComponent) const
-    {
-        return KratosComponents<Variable<double>>::Get(rSource.Name() + "_" + rComponent);
-    }
-
     virtual inline void SetTimeFactors(ModelPart& rModelPart)
     {
         mDeltaTime = rModelPart.GetProcessInfo()[DELTA_TIME];
@@ -363,8 +358,8 @@ class GeoMechanicsTimeIntegrationScheme : public Scheme<TSparseSpace, TDenseSpac
                 // We don't check for "Z", since it is optional (in case of a 2D problem)
                 std::vector<std::string> components{"X", "Y"};
                 for (const auto& component : components) {
-                    const auto& variable_component =
-                        GetComponentFromVectorVariable(r_second_order_vector_variable.instance, component);
+                    const auto& variable_component = VariablesUtilities::GetComponentFromVectorVariable(
+                        r_second_order_vector_variable.instance.Name(), component);
                     this->CheckDof(r_node, variable_component);
                 }
             }

applications/GeoMechanicsApplication/custom_strategies/schemes/newmark_dynamic_U_Pw_scheme.hpp

@@ -17,6 +17,7 @@
 
 // Application includes
 #include "custom_strategies/schemes/newmark_quasistatic_U_Pw_scheme.hpp"
+#include "custom_utilities/variables_utilities.hpp"
 #include "geo_mechanics_application_variables.h"
 
 namespace Kratos
@@ -80,15 +81,15 @@ class NewmarkDynamicUPwScheme : public GeneralizedNewmarkScheme<TSparseSpace, TD
         const std::vector<std::string> components = {"X", "Y", "Z"};
 
         for (const auto& component : components) {
-            const auto& instance_component =
-                this->GetComponentFromVectorVariable(rSecondOrderVariables.instance, component);
+            const auto& instance_component = VariablesUtilities::GetComponentFromVectorVariable(
+                rSecondOrderVariables.instance.Name(), component);
 
             if (!rNode.HasDofFor(instance_component)) continue;
 
-            const auto& first_time_derivative_component = this->GetComponentFromVectorVariable(
-                rSecondOrderVariables.first_time_derivative, component);
-            const auto& second_time_derivative_component = this->GetComponentFromVectorVariable(
-                rSecondOrderVariables.second_time_derivative, component);
+            const auto& first_time_derivative_component = VariablesUtilities::GetComponentFromVectorVariable(
+                rSecondOrderVariables.first_time_derivative.Name(), component);
+            const auto& second_time_derivative_component = VariablesUtilities::GetComponentFromVectorVariable(
+                rSecondOrderVariables.second_time_derivative.Name(), component);
 
             const double previous_variable = rNode.FastGetSolutionStepValue(instance_component, 1);
             const double current_first_time_derivative =

applications/GeoMechanicsApplication/custom_utilities/node_utilities.cpp

@@ -0,0 +1,41 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Richard Faasse
+//
+#include "node_utilities.h"
+#include "includes/node.h"
+#include "variables_utilities.hpp"
+
+#include <boost/range/combine.hpp>
+#include <string>
+#include <vector>
+
+namespace Kratos
+{
+
+void NodeUtilities::AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode,
+                                                                    const Variable<array_1d<double, 3>>& rDestinationVariable,
+                                                                    const array_1d<double, 3>& rNewValues)
+{
+    const std::vector<std::string> components = {"X", "Y", "Z"};
+    for (const auto& zipped : boost::combine(rNewValues, components)) {
+        double      new_value = 0.0;
+        std::string component;
+        boost::tie(new_value, component) = zipped;
+
+        if (const auto& component_variable =
+                VariablesUtilities::GetComponentFromVectorVariable(rDestinationVariable.Name(), component);
+            !rNode.IsFixed(component_variable)) {
+            rNode.FastGetSolutionStepValue(component_variable, 0) = new_value;
+        }
+    }
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/node_utilities.h

@@ -0,0 +1,31 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Richard Faasse
+//
+
+#pragma once
+
+#include "containers/array_1d.h"
+#include "containers/variable.h"
+
+namespace Kratos
+{
+
+class Node;
+
+class KRATOS_API(GEO_MECHANICS_APPLICATION) NodeUtilities
+{
+public:
+    static void AssignUpdatedVectorVariableToNonFixedComponents(Node& rNode,
+                                                                const Variable<array_1d<double, 3>>& rDestinationVariable,
+                                                                const array_1d<double, 3>& rNewValues);
+};
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/variables_utilities.cpp

@@ -0,0 +1,32 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//
+//  Main authors:    Richard Faasse
+//
+
+#pragma once
+
+#include "variables_utilities.hpp"
+#include "includes/kratos_components.h"
+
+namespace Kratos
+{
+
+const Variable<double>& VariablesUtilities::GetComponentFromVectorVariable(const std::string& rSourceVariableName,
+                                                                           const std::string& rComponent)
+{
+    const auto component_name = rSourceVariableName + "_" + rComponent;
+    KRATOS_ERROR_IF_NOT(KratosComponents<Variable<double>>::Has(component_name))
+        << "The component \"" << component_name << "\" is not registered!" << std::endl;
+
+    return KratosComponents<Variable<double>>::Get(component_name);
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/variables_utilities.hpp

@@ -20,19 +20,19 @@
 namespace Kratos
 {
 
-class VariablesUtilities
+class KRATOS_API(GEO_MECHANICS_APPLICATION) VariablesUtilities
 {
 public:
     template <typename GeometryType, typename OutputIt>
-    static OutputIt GetNodalValues(const GeometryType& rGeometry,
-                                   const Variable<double>& rNodalVariable,
-                                   OutputIt FirstOut)
+    static OutputIt GetNodalValues(const GeometryType& rGeometry, const Variable<double>& rNodalVariable, OutputIt FirstOut)
     {
-        return std::transform(rGeometry.begin(), rGeometry.end(), FirstOut,
-                              [&rNodalVariable](const auto& node) {
-                                  return node.FastGetSolutionStepValue(rNodalVariable);
-                              });
+        return std::transform(rGeometry.begin(), rGeometry.end(), FirstOut, [&rNodalVariable](const auto& node) {
+            return node.FastGetSolutionStepValue(rNodalVariable);
+        });
     }
+
+    static const Variable<double>& GetComponentFromVectorVariable(const std::string& rSourceVariableName,
+                                                                  const std::string& rComponent);
 };
 
-}
+} // namespace Kratos

applications/GeoMechanicsApplication/tests/cpp_tests/test_backward_euler_UPw_scheme.cpp

@@ -115,4 +115,69 @@ KRATOS_TEST_CASE_IN_SUITE(BackwardEulerUPwSchemePredict_UpdatesVariablesDerivati
         expected_dt_water_pressure);
 }
 
+KRATOS_TEST_CASE_IN_SUITE(BackwardEulerUPwSchemeUpdate_DoesNotUpdateFixedScalarVariable, KratosGeoMechanicsFastSuite)
+{
+    BackwardEulerUPwSchemeTester tester;
+
+    tester.mScheme.Initialize(tester.GetModelPart()); // This is needed to set the time factors
+    ModelPart::DofsArrayType dof_set;
+    CompressedMatrix         A;
+    Vector                   Dx;
+    Vector                   b;
+    tester.GetModelPart().Nodes()[0].Fix(DT_WATER_PRESSURE);
+
+    tester.mScheme.Update(tester.GetModelPart(), dof_set, A, Dx, b);
+
+    // should be the same as the original
+    const auto expected_dt_water_pressure = 0.0;
+
+    const auto actual_dt_water_pressure =
+        tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(DT_WATER_PRESSURE, 0);
+    KRATOS_EXPECT_DOUBLE_EQ(actual_dt_water_pressure, expected_dt_water_pressure);
+}
+
+KRATOS_TEST_CASE_IN_SUITE(BackwardEulerUPwSchemeUpdate_DoesNotUpdateFixedSecondDerivativeVectorVariable,
+                          KratosGeoMechanicsFastSuite)
+{
+    BackwardEulerUPwSchemeTester tester;
+
+    tester.mScheme.Initialize(tester.GetModelPart()); // This is needed to set the time factors
+    ModelPart::DofsArrayType dof_set;
+    CompressedMatrix         A;
+    Vector                   Dx;
+    Vector                   b;
+    tester.GetModelPart().Nodes()[0].Fix(ACCELERATION_X);
+    tester.GetModelPart().Nodes()[0].Fix(ACCELERATION_Z);
+
+    tester.mScheme.Update(tester.GetModelPart(), dof_set, A, Dx, b);
+
+    // first and last term should be the same as original, while the middle value is updated
+    const auto expected_acceleration = Kratos::array_1d<double, 3>{0.0, -1.0, 0.0};
+
+    const auto actual_acceleration =
+        tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(ACCELERATION, 0);
+    KRATOS_EXPECT_VECTOR_NEAR(actual_acceleration, expected_acceleration, 1e-6)
+}
+
+KRATOS_TEST_CASE_IN_SUITE(BackwardEulerUPwSchemeUpdate_DoesNotUpdateFixedFirstDerivativeVectorVariable,
+                          KratosGeoMechanicsFastSuite)
+{
+    BackwardEulerUPwSchemeTester tester;
+
+    tester.mScheme.Initialize(tester.GetModelPart()); // This is needed to set the time factors
+    ModelPart::DofsArrayType dof_set;
+    CompressedMatrix         A;
+    Vector                   Dx;
+    Vector                   b;
+    tester.GetModelPart().Nodes()[0].Fix(VELOCITY_Y);
+
+    tester.mScheme.Update(tester.GetModelPart(), dof_set, A, Dx, b);
+
+    // first and last term should be updated, while the middle value is fixed
+    const auto expected_velocity = Kratos::array_1d<double, 3>{-1.75, 0.0, -2.25};
+
+    const auto actual_velocity = tester.GetModelPart().Nodes()[0].FastGetSolutionStepValue(VELOCITY, 0);
+    KRATOS_EXPECT_VECTOR_NEAR(actual_velocity, expected_velocity, 1e-6)
+}
+
 } // namespace Kratos::Testing