From a6aeead2f8d2be1e97cf6201aacdbaf268820992 Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Tue, 3 Sep 2024 16:19:34 +0200 Subject: [PATCH] Add a linear elastic constitutive law for line interfaces (#12654) Added a linear elastic constitutive law that uses an incremental formulation for line interface elements. This new constitutive law will eventually replace the existing classes `LinearElastic2DInterfaceLaw` and `LinearElastic3DInterfaceLaw`. So far it aims at 2D models only. In addition to the `README` file, the unit tests of class `GeoIncrementalLinearElasticInterfaceLaw` document its intended behavior. Other changes in this PR include adding two new variables to the GeoMechanicsApplication: `INTERFACE_NORMAL_STIFFNESS` and `INTERFACE_SHEAR_STIFFNESS`. --- .../custom_constitutive/README.md | 49 +++ ...cremental_linear_elastic_interface_law.cpp | 131 ++++++++ ...incremental_linear_elastic_interface_law.h | 53 ++++ .../geo_mechanics_application.cpp | 3 + .../geo_mechanics_application_variables.cpp | 3 + .../geo_mechanics_application_variables.h | 4 + ...cremental_linear_elastic_interface_law.cpp | 281 ++++++++++++++++++ 7 files changed, 524 insertions(+) create mode 100644 applications/GeoMechanicsApplication/custom_constitutive/README.md create mode 100644 applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp create mode 100644 applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.h create mode 100644 applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_interface_law.cpp diff --git a/applications/GeoMechanicsApplication/custom_constitutive/README.md b/applications/GeoMechanicsApplication/custom_constitutive/README.md new file mode 100644 index 000000000000..39904c6762df --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_constitutive/README.md @@ -0,0 +1,49 @@ +# Constitutive laws + + +## Incremental linear elastic interface law + +This constitutive law for an interface element linearly relates increments of tractions $\Delta \tau$ to increments of relative displacement $\Delta \Delta u$. +Relative displacement for interface element is the differential motion between the two sides of the interface. As a +consequence the relative displacement has unit of length $[\mathrm{L}]$ and the stiffness has unit of force over cubic length $[\mathrm{F/L^3}]$. + +### Relative displacement and traction +In 2D plane strain $y$ is the opening/closing direction of the interface, while differential motion in the tangential direction +gives shear. Similar to the continuum, the normal behavior is placed first in the relative displacement and traction vectors. The shear behavior +follows after the normal motion or traction in these vectors. + +$$ \Delta u = \begin{bmatrix} \Delta u_y \\ \Delta u_x \end{bmatrix} $$ + +$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$ + +Where: +* $\Delta u_y$: Relative displacement in the $y$-direction (normal to the interface). +* $\Delta u_x$: Relative displacement in the $x$-direction (tangential to the interface). +* $\tau_{yy}$: Traction in the $y$-direction (normal to the interface). +* $\tau_{xy}$: Shear traction in the $x$-direction (tangential to the interface). + +### 2D Elastic constitutive tensor + +The elastic behavior of the interface is characterized by the 2D elastic constitutive tensor $C$, which relates the traction and relative displacement vectors. The constitutive tensor in 2D is expressed as: + +$$ C = \begin{bmatrix} C_{yy} & 0 \\ + 0 & C_{xy} \end{bmatrix}$$ + +Where: +* $C$: Represents the 2D elastic constitutive tensor. +* $C_{yy}$: Represents the `INTERFACE_NORMAL_STIFFNESS`, which characterizes the stiffness in the normal direction. +* $C_{xy}$: Represents the `INTERFACE_SHEAR_STIFFNESS`, which characterizes the stiffness in the tangential direction. + +Both stiffness values have dimension $\mathrm{F/L^3}$. + +### Incremental relation + +The relation between the traction and the relative displacement for a given time increment is given by the following incremental equation: + +$$ \tau_{t + \Delta t} = \tau_t + C \cdot \Delta \Delta u $$ + +Where: +* $\tau_{t + \Delta t}$: The traction at the updated time $t + \Delta t$. +* $\tau_t$: The traction at the current time $t$. +* $C$: The elastic constitutive tensor. +* $\Delta \Delta u$: The incremental relative displacement vector. diff --git a/applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp b/applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp new file mode 100644 index 000000000000..2cc6039382eb --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp @@ -0,0 +1,131 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Wijtze Pieter Kikstra +// Anne van de Graaf +// + +#include "incremental_linear_elastic_interface_law.h" +#include "custom_geometries/line_interface_geometry.h" +#include "geo_mechanics_application_constants.h" +#include "geo_mechanics_application_variables.h" + +namespace Kratos +{ + +ConstitutiveLaw::Pointer GeoIncrementalLinearElasticInterfaceLaw::Clone() const +{ + return std::make_shared(*this); +} + +ConstitutiveLaw::SizeType GeoIncrementalLinearElasticInterfaceLaw::WorkingSpaceDimension() +{ + return 2; +} + +ConstitutiveLaw::SizeType GeoIncrementalLinearElasticInterfaceLaw::GetStrainSize() const +{ + return VOIGT_SIZE_2D_INTERFACE; +} + +Vector& GeoIncrementalLinearElasticInterfaceLaw::GetValue(const Variable& rThisVariable, Vector& rValue) +{ + if (rThisVariable == STRAIN) { + rValue = mPreviousRelativeDisplacement; + } else if (rThisVariable == CAUCHY_STRESS_VECTOR) { + rValue = mPreviousTraction; + } + + return rValue; +} + +ConstitutiveLaw::StressMeasure GeoIncrementalLinearElasticInterfaceLaw::GetStressMeasure() +{ + return ConstitutiveLaw::StressMeasure_Cauchy; +} + +bool GeoIncrementalLinearElasticInterfaceLaw::IsIncremental() { return true; } + +void GeoIncrementalLinearElasticInterfaceLaw::InitializeMaterial(const Properties&, + const ConstitutiveLaw::GeometryType&, + const Vector&) +{ + mPreviousRelativeDisplacement = + HasInitialState() ? GetInitialState().GetInitialStrainVector() : ZeroVector{GetStrainSize()}; + mPreviousTraction = + HasInitialState() ? GetInitialState().GetInitialStressVector() : ZeroVector{GetStrainSize()}; +} + +void GeoIncrementalLinearElasticInterfaceLaw::CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) +{ + rValues.GetStressVector() = + mPreviousTraction + + prod(MakeConstitutiveMatrix(rValues.GetMaterialProperties()[INTERFACE_NORMAL_STIFFNESS], + rValues.GetMaterialProperties()[INTERFACE_SHEAR_STIFFNESS]), + rValues.GetStrainVector() - mPreviousRelativeDisplacement); +} + +bool GeoIncrementalLinearElasticInterfaceLaw::RequiresInitializeMaterialResponse() { return false; } + +void GeoIncrementalLinearElasticInterfaceLaw::FinalizeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues) +{ + mPreviousRelativeDisplacement = rValues.GetStrainVector(); + mPreviousTraction = rValues.GetStressVector(); +} + +int GeoIncrementalLinearElasticInterfaceLaw::Check(const Properties& rMaterialProperties, + const ConstitutiveLaw::GeometryType& rElementGeometry, + const ProcessInfo& rCurrentProcessInfo) const +{ + const auto result = BaseType::Check(rMaterialProperties, rElementGeometry, rCurrentProcessInfo); + + KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_NORMAL_STIFFNESS)) + << "No interface normal stiffness is defined" << std::endl; + + KRATOS_ERROR_IF_NOT(rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] > 0.0) + << "Interface normal stiffness must be positive, but got " + << rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] << std::endl; + + KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_SHEAR_STIFFNESS)) + << "No interface shear stiffness is defined" << std::endl; + + KRATOS_ERROR_IF_NOT(rMaterialProperties[INTERFACE_SHEAR_STIFFNESS] > 0.0) + << "Interface shear stiffness must be positive, but got " + << rMaterialProperties[INTERFACE_SHEAR_STIFFNESS] << std::endl; + + KRATOS_ERROR_IF_NOT(dynamic_cast(&rElementGeometry)) + << "Expected a line interface geometry, but got " << rElementGeometry.Info() << std::endl; + + return result; +} + +void GeoIncrementalLinearElasticInterfaceLaw::save(Serializer& rSerializer) const +{ + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseType) + rSerializer.save("PreviousRelativeDisplacement", mPreviousRelativeDisplacement); + rSerializer.save("PreviousTraction", mPreviousTraction); +} + +void GeoIncrementalLinearElasticInterfaceLaw::load(Serializer& rSerializer) +{ + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseType) + rSerializer.load("PreviousRelativeDisplacement", mPreviousRelativeDisplacement); + rSerializer.load("PreviousTraction", mPreviousTraction); +} + +Matrix GeoIncrementalLinearElasticInterfaceLaw::MakeConstitutiveMatrix(double NormalStiffness, + double ShearStiffness) const +{ + auto result = Matrix{ZeroMatrix{GetStrainSize(), GetStrainSize()}}; + result(0, 0) = NormalStiffness; + result(1, 1) = ShearStiffness; + return result; +} + +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.h b/applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.h new file mode 100644 index 000000000000..b3ea3e660982 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.h @@ -0,0 +1,53 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Wijtze Pieter Kikstra +// Anne van de Graaf +// + +#pragma once + +#include "includes/constitutive_law.h" +#include "includes/kratos_export_api.h" + +namespace Kratos +{ + +class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoIncrementalLinearElasticInterfaceLaw : public ConstitutiveLaw +{ +public: + using BaseType = ConstitutiveLaw; + + [[nodiscard]] Pointer Clone() const override; + SizeType WorkingSpaceDimension() override; + [[nodiscard]] SizeType GetStrainSize() const override; + Vector& GetValue(const Variable& rThisVariable, Vector& rValue) override; + using BaseType::GetValue; + StressMeasure GetStressMeasure() override; + bool IsIncremental() override; + void InitializeMaterial(const Properties&, const GeometryType&, const Vector&) override; + void CalculateMaterialResponseCauchy(Parameters& rValues) override; + bool RequiresInitializeMaterialResponse() override; + void FinalizeMaterialResponseCauchy(Parameters& rValues) override; + int Check(const Properties& rMaterialProperties, + const GeometryType& rElementGeometry, + const ProcessInfo& rCurrentProcessInfo) const override; + +private: + friend class Serializer; + void save(Serializer& rSerializer) const override; + void load(Serializer& rSerializer) override; + + [[nodiscard]] Matrix MakeConstitutiveMatrix(double NormalStiffness, double ShearStiffness) const; + + Vector mPreviousRelativeDisplacement; + Vector mPreviousTraction; +}; + +} // namespace Kratos diff --git a/applications/GeoMechanicsApplication/geo_mechanics_application.cpp b/applications/GeoMechanicsApplication/geo_mechanics_application.cpp index 79a23d09b896..6c4e81230415 100644 --- a/applications/GeoMechanicsApplication/geo_mechanics_application.cpp +++ b/applications/GeoMechanicsApplication/geo_mechanics_application.cpp @@ -591,5 +591,8 @@ void KratosGeoMechanicsApplication::Register() KRATOS_REGISTER_VARIABLE(STRAINS_OF_PIECEWISE_LINEAR_LAW) KRATOS_REGISTER_VARIABLE(STRESSES_OF_PIECEWISE_LINEAR_LAW) + + KRATOS_REGISTER_VARIABLE(INTERFACE_NORMAL_STIFFNESS) + KRATOS_REGISTER_VARIABLE(INTERFACE_SHEAR_STIFFNESS) } } // namespace Kratos. diff --git a/applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp b/applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp index da495caa7fb4..a5b2a9409fd7 100644 --- a/applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp +++ b/applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp @@ -243,4 +243,7 @@ KRATOS_CREATE_VARIABLE(double, STATE_VARIABLE_50) KRATOS_CREATE_VARIABLE(Vector, STRAINS_OF_PIECEWISE_LINEAR_LAW) KRATOS_CREATE_VARIABLE(Vector, STRESSES_OF_PIECEWISE_LINEAR_LAW) +KRATOS_CREATE_VARIABLE(double, INTERFACE_NORMAL_STIFFNESS) +KRATOS_CREATE_VARIABLE(double, INTERFACE_SHEAR_STIFFNESS) + } // namespace Kratos diff --git a/applications/GeoMechanicsApplication/geo_mechanics_application_variables.h b/applications/GeoMechanicsApplication/geo_mechanics_application_variables.h index 83d0d7a2df55..67585a3bdc32 100644 --- a/applications/GeoMechanicsApplication/geo_mechanics_application_variables.h +++ b/applications/GeoMechanicsApplication/geo_mechanics_application_variables.h @@ -259,6 +259,10 @@ KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, STATE_VARI KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, STRAINS_OF_PIECEWISE_LINEAR_LAW) KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, STRESSES_OF_PIECEWISE_LINEAR_LAW) + +KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, INTERFACE_NORMAL_STIFFNESS) +KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, INTERFACE_SHEAR_STIFFNESS) + } // namespace Kratos #endif /* KRATOS_GEO_MECHANICS_APPLICATION_VARIABLES_H_INCLUDED */ diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_interface_law.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_interface_law.cpp new file mode 100644 index 000000000000..28381be1a623 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/custom_constitutive/test_incremental_linear_elastic_interface_law.cpp @@ -0,0 +1,281 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Wijtze Pieter Kikstra +// Anne van de Graaf +// + +#include "custom_constitutive/incremental_linear_elastic_interface_law.h" +#include "custom_geometries/line_interface_geometry.h" +#include "geo_mechanics_application_variables.h" +#include "geometries/line_3d_2.h" +#include "includes/checks.h" +#include "includes/serializer.h" +#include "tests/cpp_tests/geo_mechanics_fast_suite.h" + +#include +#include +#include + +using namespace Kratos; +using namespace std::string_literals; + +namespace +{ + +constexpr auto absolute_tolerance = 1.0e-6; +constexpr auto relative_tolerance = 1.0e-6; + +} // namespace + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesHas2DWorkingSpace, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + + KRATOS_EXPECT_EQ(law.WorkingSpaceDimension(), 2); +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesHasStrainSizeOfTwo, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + + KRATOS_EXPECT_EQ(law.GetStrainSize(), 2); +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesUsesCauchyStressMeasure, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + + KRATOS_EXPECT_EQ(law.GetStressMeasure(), ConstitutiveLaw::StressMeasure_Cauchy); +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesIsIncremental, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + + KRATOS_EXPECT_TRUE(law.IsIncremental()) +} + +KRATOS_TEST_CASE_IN_SUITE(CloneOfLinearElasticLawForInterfacesIsIndependentOfOriginal, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto original_law = GeoIncrementalLinearElasticInterfaceLaw{}; + auto p_cloned_law = original_law.Clone(); + + KRATOS_EXPECT_NE(p_cloned_law.get(), nullptr); + KRATOS_EXPECT_NE(p_cloned_law.get(), &original_law); + KRATOS_EXPECT_NE(dynamic_cast(p_cloned_law.get()), nullptr); +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesDoesNotRequireInitializationOfMaterialResponse, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + + KRATOS_EXPECT_FALSE(law.RequiresInitializeMaterialResponse()) +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesChecksForCorrectMaterialProperties, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + auto properties = Properties{}; + const auto geometry = LineInterfaceGeometry{}; + const auto process_info = ProcessInfo{}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN(law.Check(properties, geometry, process_info), + "No interface normal stiffness is defined") + + properties[INTERFACE_NORMAL_STIFFNESS] = -5.0; + KRATOS_EXPECT_EXCEPTION_IS_THROWN(law.Check(properties, geometry, process_info), + "Interface normal stiffness must be positive, but got -5") + + properties[INTERFACE_NORMAL_STIFFNESS] = 0.0; + KRATOS_EXPECT_EXCEPTION_IS_THROWN(law.Check(properties, geometry, process_info), + "Interface normal stiffness must be positive, but got 0") + + properties[INTERFACE_NORMAL_STIFFNESS] = 5.0; + KRATOS_EXPECT_EXCEPTION_IS_THROWN(law.Check(properties, geometry, process_info), + "No interface shear stiffness is defined") + + properties[INTERFACE_SHEAR_STIFFNESS] = -2.5; + KRATOS_EXPECT_EXCEPTION_IS_THROWN(law.Check(properties, geometry, process_info), + "Interface shear stiffness must be positive, but got -2.5") + + properties[INTERFACE_SHEAR_STIFFNESS] = 0.0; + KRATOS_EXPECT_EXCEPTION_IS_THROWN(law.Check(properties, geometry, process_info), + "Interface shear stiffness must be positive, but got 0") + + properties[INTERFACE_SHEAR_STIFFNESS] = 2.5; + KRATOS_EXPECT_EQ(law.Check(properties, geometry, process_info), 0); +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesChecksForCorrectGeometry, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + const auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + auto properties = Properties{}; + properties[INTERFACE_NORMAL_STIFFNESS] = 5.0; + properties[INTERFACE_SHEAR_STIFFNESS] = 2.5; + const auto node1 = make_intrusive(1, 0.0, 0.0, 0.0); + const auto node2 = make_intrusive(2, 5.0, 5.0, 0.0); + const auto line_3d_geometry = Line3D2{node1, node2}; + const auto process_info = ProcessInfo{}; + + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + law.Check(properties, line_3d_geometry, process_info), + "Expected a line interface geometry, but got 1 dimensional line with 2 nodes in 3D space") +} + +KRATOS_TEST_CASE_IN_SUITE(WhenNoInitialStateIsGivenStartWithZeroRelativeDisplacementAndZeroTraction, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + + const auto dummy_properties = Properties{}; + const auto dummy_geometry = Geometry{}; + const auto dummy_shape_function_values = Vector{}; + law.InitializeMaterial(dummy_properties, dummy_geometry, dummy_shape_function_values); + + auto value = Vector{}; + law.GetValue(STRAIN, value); + const auto zero_vector = Vector{ZeroVector{2}}; + KRATOS_EXPECT_VECTOR_NEAR(value, zero_vector, absolute_tolerance) + law.GetValue(CAUCHY_STRESS_VECTOR, value); + KRATOS_EXPECT_VECTOR_NEAR(value, zero_vector, absolute_tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(WhenAnInitialStateIsGivenStartFromThereAfterMaterialInitialization, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + const auto initial_relative_displacement = Vector{ScalarVector{2, 0.5}}; + const auto initial_traction = Vector{ScalarVector{2, 30.0}}; + auto p_initial_state = make_intrusive(initial_relative_displacement, initial_traction); + law.SetInitialState(p_initial_state); + + const auto dummy_properties = Properties{}; + const auto dummy_geometry = Geometry{}; + const auto dummy_shape_function_values = Vector{}; + law.InitializeMaterial(dummy_properties, dummy_geometry, dummy_shape_function_values); + + auto value = Vector{}; + law.GetValue(STRAIN, value); + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(value, initial_relative_displacement, relative_tolerance) + law.GetValue(CAUCHY_STRESS_VECTOR, value); + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(value, initial_traction, relative_tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(ComputedIncrementalTractionIsProductOfIncrementalRelativeDisplacementAndStiffness, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law_parameters = ConstitutiveLaw::Parameters{}; + auto relative_displacement = Vector{2}; + relative_displacement <<= 0.1, 0.3; + law_parameters.SetStrainVector(relative_displacement); + auto traction = Vector{ZeroVector{2}}; + law_parameters.SetStressVector(traction); + auto properties = Properties{}; + properties[INTERFACE_NORMAL_STIFFNESS] = 20.0; + properties[INTERFACE_SHEAR_STIFFNESS] = 10.0; + law_parameters.SetMaterialProperties(properties); + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + const auto dummy_geometry = Geometry{}; + const auto dummy_shape_function_values = Vector{}; + law.InitializeMaterial(properties, dummy_geometry, dummy_shape_function_values); + + law.CalculateMaterialResponseCauchy(law_parameters); + + auto expected_traction = Vector{2}; + expected_traction <<= 2.0, 3.0; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(law_parameters.GetStressVector(), expected_traction, relative_tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(ComputedTractionIsSumOfPreviousTractionAndTractionIncrement, KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law_parameters = ConstitutiveLaw::Parameters{}; + auto relative_displacement = Vector{ZeroVector{2}}; + law_parameters.SetStrainVector(relative_displacement); + auto traction = Vector{ZeroVector{2}}; + law_parameters.SetStressVector(traction); + auto properties = Properties{}; + properties[INTERFACE_NORMAL_STIFFNESS] = 20.0; + properties[INTERFACE_SHEAR_STIFFNESS] = 10.0; + law_parameters.SetMaterialProperties(properties); + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + const auto initial_relative_displacement = Vector{ScalarVector{2, 5.0}}; + const auto initial_traction = Vector{ScalarVector{2, 30.0}}; + auto p_initial_state = make_intrusive(initial_relative_displacement, initial_traction); + law.SetInitialState(p_initial_state); + const auto dummy_geometry = Geometry{}; + const auto dummy_shape_function_values = Vector{}; + law.InitializeMaterial(properties, dummy_geometry, dummy_shape_function_values); + + // First step + relative_displacement <<= 5.1, 5.3; + law.CalculateMaterialResponseCauchy(law_parameters); + law.FinalizeMaterialResponseCauchy(law_parameters); + + // Second step + relative_displacement <<= 5.2, 5.6; + law.CalculateMaterialResponseCauchy(law_parameters); + law.FinalizeMaterialResponseCauchy(law_parameters); + + auto expected_traction = Vector{2}; + expected_traction <<= 30.0 + (5.1 - 5.0) * 20.0 + (5.2 - 5.1) * 20.0, 30.0 + (5.3 - 5.0) * 10.0 + (5.6 - 5.3) * 10.0; + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(law_parameters.GetStressVector(), expected_traction, relative_tolerance) + auto expected_relative_displacement = Vector{2}; + expected_relative_displacement <<= 5.0 + (5.1 - 5.0) + (5.2 - 5.1), 5.0 + (5.3 - 5.0) + (5.6 - 5.3); + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(law_parameters.GetStrainVector(), expected_relative_displacement, relative_tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(LinearElasticLawForInterfacesCanBeSavedToAndLoadedFromASerializer, + KratosGeoMechanicsFastSuiteWithoutKernel) +{ + auto law = GeoIncrementalLinearElasticInterfaceLaw{}; + const auto initial_relative_displacement = Vector{ScalarVector{2, 0.5}}; + const auto initial_traction = Vector{ScalarVector{2, 30.0}}; + auto p_initial_state = make_intrusive(initial_relative_displacement, initial_traction); + law.SetInitialState(p_initial_state); + + const auto dummy_properties = Properties{}; + const auto dummy_geometry = Geometry{}; + const auto dummy_shape_function_values = Vector{}; + law.InitializeMaterial(dummy_properties, dummy_geometry, dummy_shape_function_values); + + auto law_parameters = ConstitutiveLaw::Parameters{}; + auto relative_displacement = Vector{2}; + relative_displacement <<= 0.1, 0.3; + law_parameters.SetStrainVector(relative_displacement); + auto traction = Vector{2}; + traction <<= 20.0, 45.0; + law_parameters.SetStressVector(traction); + law.FinalizeMaterialResponseCauchy(law_parameters); + + auto serializer = Serializer{new std::stringstream{}}; + const auto tag = "test_tag"s; + serializer.save(tag, law); + + auto restored_law = GeoIncrementalLinearElasticInterfaceLaw{}; + serializer.load(tag, restored_law); + + auto value = Vector{}; + restored_law.GetValue(STRAIN, value); + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(value, relative_displacement, relative_tolerance) + restored_law.GetValue(CAUCHY_STRESS_VECTOR, value); + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(value, traction, relative_tolerance) + KRATOS_EXPECT_TRUE(restored_law.HasInitialState()) + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(restored_law.GetInitialState().GetInitialStrainVector(), + initial_relative_displacement, relative_tolerance) + KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(restored_law.GetInitialState().GetInitialStressVector(), + initial_traction, relative_tolerance) +} + +} // namespace Kratos::Testing