[GeoMechanicsApplication] Implement well constitutive behaviour for line thermal element

applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp

@@ -19,27 +19,16 @@
 namespace Kratos
 {
 
-GeoThermalDispersionLaw::GeoThermalDispersionLaw() : mNumberOfDimensions{2} {}
-
 GeoThermalDispersionLaw::GeoThermalDispersionLaw(std::size_t NumberOfDimensions)
     : mNumberOfDimensions{NumberOfDimensions}
 {
     KRATOS_ERROR_IF(mNumberOfDimensions < 1 || mNumberOfDimensions > 3)
         << "Got invalid number of dimensions: " << mNumberOfDimensions << std::endl;
 }
 
-ConstitutiveLaw::Pointer GeoThermalDispersionLaw::Clone() const
-{
-    return Kratos::make_shared<GeoThermalDispersionLaw>(*this);
-}
-
-SizeType GeoThermalDispersionLaw::WorkingSpaceDimension() { return mNumberOfDimensions; }
-
 Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Properties& rProp,
                                                                  const ProcessInfo& rProcessInfo) const
 {
-    KRATOS_TRY
-
     Matrix result = ZeroMatrix(mNumberOfDimensions, mNumberOfDimensions);
 
     RetentionLaw::Parameters parameters(rProp, rProcessInfo);
@@ -53,7 +42,7 @@ Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Propertie
                    water_fraction * rProp[THERMAL_CONDUCTIVITY_WATER];
 
     if (mNumberOfDimensions >= 2) {
-        const auto y  = static_cast<int>(indexThermalFlux::Y);
+        const auto y = static_cast<int>(indexThermalFlux::Y);
         result(x, y) = solid_fraction * rProp[THERMAL_CONDUCTIVITY_SOLID_XY];
         result(y, y) = solid_fraction * rProp[THERMAL_CONDUCTIVITY_SOLID_YY] +
                        water_fraction * rProp[THERMAL_CONDUCTIVITY_WATER];
@@ -70,8 +59,6 @@ Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Propertie
     }
 
     return result;
-
-    KRATOS_CATCH("")
 }
 
 } // Namespace Kratos

applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h

@@ -14,48 +14,36 @@
 
 #pragma once
 
-#include "includes/constitutive_law.h"
+#include "custom_constitutive/thermal_law.h"
 
 namespace Kratos
 {
 
 class RetentionLaw;
 
-/**
- * @class GeoThermalDispersionLaw
- * @ingroup GeoMechanicsApplication
- * @brief This class defines the thermal dispersion for heat cases
- */
-class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalDispersionLaw : public ConstitutiveLaw
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalDispersionLaw : public GeoThermalLaw
 {
 public:
-    /// Counted pointer of GeoThermalDispersionLaw
     KRATOS_CLASS_POINTER_DEFINITION(GeoThermalDispersionLaw);
 
-    GeoThermalDispersionLaw();
+    GeoThermalDispersionLaw() = default;
 
     explicit GeoThermalDispersionLaw(SizeType NumberOfDimensions);
 
-    ConstitutiveLaw::Pointer Clone() const override;
-
-    SizeType WorkingSpaceDimension() override;
-
-    Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const;
+    Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const override;
 
 private:
-    SizeType mNumberOfDimensions;
+    std::size_t mNumberOfDimensions = 2;
 
     friend class Serializer;
 
     void save(Serializer& rSerializer) const override
     {
-        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw)
         rSerializer.save("NumberOfDimensions", mNumberOfDimensions);
     }
 
     void load(Serializer& rSerializer) override
     {
-        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw)
         rSerializer.load("NumberOfDimensions", mNumberOfDimensions);
     }
 };

applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp

@@ -0,0 +1,25 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   John van Esch
+//
+
+#include "custom_constitutive/thermal_filter_law.h"
+#include "geo_mechanics_application_variables.h"
+
+namespace Kratos
+{
+
+Matrix GeoThermalFilterLaw::CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const
+{
+    return ScalarMatrix(1, 1, rProp[THERMAL_CONDUCTIVITY_WATER]);
+}
+
+} // Namespace Kratos

applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.h

@@ -0,0 +1,41 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   John van Esch
+//
+
+#pragma once
+
+#include "custom_constitutive/thermal_law.h"
+
+namespace Kratos
+{
+
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalFilterLaw : public GeoThermalLaw
+{
+public:
+    KRATOS_CLASS_POINTER_DEFINITION(GeoThermalFilterLaw);
+
+    Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const override;
+
+private:
+    friend class Serializer;
+
+    void save(Serializer&) const override
+    {
+        // Nothing to serialize, since there are no data members, and its base class is abstract
+    }
+
+    void load(Serializer&) override
+    {
+        // Nothing to serialize, since there are no data members, and its base class is abstract
+    }
+};
+} // namespace Kratos.

applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h

@@ -0,0 +1,40 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   John van Esch
+//                   Gennady Markelov
+//
+
+#pragma once
+
+#include "geo_mechanics_application_variables.h"
+#include "includes/serializer.h"
+
+namespace Kratos
+{
+
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalLaw
+{
+public:
+    KRATOS_CLASS_POINTER_DEFINITION(GeoThermalLaw);
+
+    virtual ~GeoThermalLaw() = default;
+
+    virtual Matrix CalculateThermalDispersionMatrix(const Properties&  rProp,
+                                                    const ProcessInfo& rProcessInfo) const = 0;
+
+private:
+    friend class Serializer;
+
+    virtual void save(Serializer& rSerializer) const = 0;
+
+    virtual void load(Serializer& rSerializer) = 0;
+};
+} // namespace Kratos.

applications/GeoMechanicsApplication/custom_elements/transient_thermal_element.h

@@ -15,6 +15,7 @@
 #pragma once
 
 #include "custom_constitutive/thermal_dispersion_law.h"
+#include "custom_constitutive/thermal_filter_law.h"
 #include "custom_retention/retention_law_factory.h"
 #include "custom_utilities/dof_utilities.h"
 #include "geo_mechanics_application_variables.h"
@@ -268,10 +269,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Ele
         const Vector&                                    rIntegrationCoefficients,
         const ProcessInfo&                               rCurrentProcessInfo) const
     {
-        const std::size_t number_of_dimensions = GetGeometry().LocalSpaceDimension();
-
-        GeoThermalDispersionLaw law{number_of_dimensions};
-        const auto constitutive_matrix = law.CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo);
+        const auto law = CreateThermalLaw();
+        const auto constitutive_matrix = law->CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo);
 
         auto result = BoundedMatrix<double, TNumNodes, TNumNodes>{ZeroMatrix{TNumNodes, TNumNodes}};
         for (unsigned int integration_point_index = 0;
@@ -320,6 +319,29 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Ele
         return result;
     }
 
+    std::unique_ptr<GeoThermalLaw> CreateThermalLaw() const
+    {
+        const std::size_t number_of_dimensions = GetGeometry().LocalSpaceDimension();
+        std::unique_ptr<GeoThermalLaw> law;
+
+        if (GetProperties().Has(THERMAL_LAW_NAME)) {
+            const std::string& rThermalLawName = GetProperties()[THERMAL_LAW_NAME];
+            if (rThermalLawName == "GeoThermalDispersionLaw") {
+                law = std::make_unique<GeoThermalDispersionLaw>(number_of_dimensions);
+            } 
+            else if (rThermalLawName == "GeoThermalFilterLaw") {
+                law = std::make_unique<GeoThermalFilterLaw>();
+            } 
+            else {
+                KRATOS_ERROR << "Undefined THERMAL_LAW_NAME! " << rThermalLawName << std::endl;
+            }
+        } 
+        else {
+            law = std::make_unique<GeoThermalDispersionLaw>(number_of_dimensions);
+        }
+        return law;
+    }
+
     [[nodiscard]] DofsVectorType GetDofs() const
     {
         return Geo::DofUtilities::ExtractDofsFromNodes(GetGeometry(), TEMPERATURE);

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -352,8 +352,6 @@ void KratosGeoMechanicsApplication::Register() {
 
     KRATOS_REGISTER_CONSTITUTIVE_LAW("LinearElastic2DBeamLaw",        mLinearElastic2DBeamLaw)
 
-    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoThermalDispersion2DLaw",     mGeoThermalDispersion2DLaw)
-    KRATOS_REGISTER_CONSTITUTIVE_LAW("GeoThermalDispersion3DLaw",     mGeoThermalDispersion3DLaw)
 
     //Register Variables
     KRATOS_REGISTER_VARIABLE( VELOCITY_COEFFICIENT )
@@ -378,6 +376,7 @@ void KratosGeoMechanicsApplication::Register() {
     KRATOS_REGISTER_VARIABLE( DT_TEMPERATURE_COEFFICIENT )
     KRATOS_REGISTER_VARIABLE( DT_TEMPERATURE )
     KRATOS_REGISTER_VARIABLE( NORMAL_HEAT_FLUX )
+    KRATOS_REGISTER_VARIABLE( THERMAL_LAW_NAME)
 
     // Variables for Micro-Climate boundary
     KRATOS_REGISTER_VARIABLE( AIR_TEMPERATURE )

applications/GeoMechanicsApplication/geo_mechanics_application.h

@@ -131,7 +131,6 @@
 #include "custom_constitutive/small_strain_umat_2D_plane_strain_law.hpp"
 #include "custom_constitutive/small_strain_umat_3D_interface_law.hpp"
 #include "custom_constitutive/small_strain_umat_3D_law.hpp"
-#include "custom_constitutive/thermal_dispersion_law.h"
 
 namespace Kratos {
 
@@ -598,9 +597,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoMechanicsApplication : publ
 
     const LinearElastic2DBeamLaw          mLinearElastic2DBeamLaw;
 
-    const GeoThermalDispersionLaw mGeoThermalDispersion2DLaw{ConstitutiveLaw::SizeType(2)};
-    const GeoThermalDispersionLaw mGeoThermalDispersion3DLaw{ConstitutiveLaw::SizeType(3)};
-
     ///@}
 
 }; // Class KratosGeoMechanicsApplication

applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp

@@ -38,6 +38,7 @@ KRATOS_CREATE_VARIABLE( double, SOLID_COMPRESSIBILITY )
 KRATOS_CREATE_VARIABLE( double, DT_TEMPERATURE_COEFFICIENT )
 KRATOS_CREATE_VARIABLE( double, DT_TEMPERATURE )
 KRATOS_CREATE_VARIABLE( double, NORMAL_HEAT_FLUX )
+KRATOS_CREATE_VARIABLE( std::string, THERMAL_LAW_NAME )
 
 // Variables for Micro-Climate boundary
 KRATOS_CREATE_VARIABLE( double, AIR_TEMPERATURE )

applications/GeoMechanicsApplication/geo_mechanics_application_variables.h

@@ -65,6 +65,7 @@ namespace Kratos
     KRATOS_DEFINE_APPLICATION_VARIABLE( GEO_MECHANICS_APPLICATION, double, DT_TEMPERATURE_COEFFICIENT )
     KRATOS_DEFINE_APPLICATION_VARIABLE( GEO_MECHANICS_APPLICATION, double, DT_TEMPERATURE )
     KRATOS_DEFINE_APPLICATION_VARIABLE( GEO_MECHANICS_APPLICATION, double, NORMAL_HEAT_FLUX )
+    KRATOS_DEFINE_APPLICATION_VARIABLE( GEO_MECHANICS_APPLICATION, std::string, THERMAL_LAW_NAME )
 
     // Variables for Micro-Climate boundary
     KRATOS_DEFINE_APPLICATION_VARIABLE( GEO_MECHANICS_APPLICATION, double, AIR_TEMPERATURE )

applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_dispersion_law.cpp

@@ -101,12 +101,9 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateThermalDispersionMatrix3D, KratosGeoMechanics
     }
 }
 
-KRATOS_TEST_CASE_IN_SUITE(GetWorkingSpaceDimension_ReturnsCorrectValue, KratosGeoMechanicsFastSuite)
+KRATOS_TEST_CASE_IN_SUITE(TestDispersionLawThrowsWhenDimensionInvalid, KratosGeoMechanicsFastSuite)
 {
-    constexpr SizeType dimension = 3;
-    GeoThermalDispersionLaw geo_thermal_dispersion_3D_law(dimension);
-
-    KRATOS_EXPECT_EQ(geo_thermal_dispersion_3D_law.WorkingSpaceDimension(), dimension);
+    KRATOS_EXPECT_EXCEPTION_IS_THROWN(GeoThermalDispersionLaw law{0}, "Got invalid number of dimensions: 0")
 }
 
 } // namespace Kratos::Testing

applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_filter_law.cpp

@@ -0,0 +1,41 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//
+
+#include "custom_constitutive/thermal_filter_law.h"
+#include "geo_mechanics_application.h"
+#include "includes/ublas_interface.h"
+#include "testing/testing.h"
+
+namespace Kratos::Testing
+{
+
+KRATOS_TEST_CASE_IN_SUITE(CalculateThermalFilterLawMatrix, KratosGeoMechanicsFastSuite)
+{
+    Model current_model;
+    auto& r_model_part = current_model.CreateModelPart("ModelPart");
+
+    auto p_cond_prop = r_model_part.CreateNewProperties(0);
+    p_cond_prop->SetValue(THERMAL_CONDUCTIVITY_WATER, 1000.0);
+
+    GeoThermalFilterLaw geo_thermal_filter_law;
+    ProcessInfo         info;
+
+    const Matrix thermal_filter_matrix = geo_thermal_filter_law.CalculateThermalDispersionMatrix(*p_cond_prop, info);
+
+    Matrix expected_solution = ScalarMatrix(1, 1, 1000.0);
+
+    constexpr double tolerance{1.0e-6};
+
+    KRATOS_EXPECT_MATRIX_NEAR(thermal_filter_matrix, expected_solution, tolerance)
+}
+
+} // namespace Kratos::Testing