[GeoMechanicsApplication] Implement well constitutive behaviour for line thermal element

applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp

@@ -28,18 +28,9 @@ GeoThermalDispersionLaw::GeoThermalDispersionLaw(std::size_t NumberOfDimensions)
         << "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 +44,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 +61,6 @@ Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Propertie
     }
 
     return result;
-
-    KRATOS_CATCH("")
 }
 
 } // Namespace Kratos

applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h

@@ -14,7 +14,7 @@
 
 #pragma once
 
-#include "includes/constitutive_law.h"
+#include "custom_constitutive/thermal_law.h"
 
 namespace Kratos
 {
@@ -26,7 +26,7 @@ class RetentionLaw;
  * @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
@@ -36,26 +36,24 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalDispersionLaw : public Con
 
     explicit GeoThermalDispersionLaw(SizeType NumberOfDimensions);
 
-    ConstitutiveLaw::Pointer Clone() const override;
+    ~GeoThermalDispersionLaw() override = default;
 
-    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 = 0;
 
     friend class Serializer;
 
     void save(Serializer& rSerializer) const override
     {
-        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw)
+        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, GeoThermalLaw)
         rSerializer.save("NumberOfDimensions", mNumberOfDimensions);
     }
 
     void load(Serializer& rSerializer) override
     {
-        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw)
+        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, GeoThermalLaw)
         rSerializer.load("NumberOfDimensions", mNumberOfDimensions);
     }
 };

applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp

@@ -0,0 +1,35 @@
+// 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
+{
+
+GeoThermalFilterLaw::GeoThermalFilterLaw() : mNumberOfDimensions{1} {}
+
+GeoThermalFilterLaw::GeoThermalFilterLaw(SizeType NumberOfDimensions)
+    : mNumberOfDimensions{NumberOfDimensions}
+{
+    KRATOS_ERROR_IF(mNumberOfDimensions != 1)
+        << "Got invalid number of dimensions. The dimension has to be 1, but got: " << mNumberOfDimensions
+        << std::endl;
+}
+
+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,57 @@
+// 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 GeoThermalFilterLaw
+ * @ingroup GeoMechanicsApplication
+ * @brief This class defines the thermal filter for heat cases
+ */
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalFilterLaw : public GeoThermalLaw
+{
+public:
+    /// Counted pointer of GeoThermalFilterLaw
+    KRATOS_CLASS_POINTER_DEFINITION(GeoThermalFilterLaw);
+
+    GeoThermalFilterLaw();
+
+    explicit GeoThermalFilterLaw(SizeType NumberOfDimensions);
+
+    ~GeoThermalFilterLaw() override = default;
+
+    Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const override;
+
+private:
+    std::size_t mNumberOfDimensions = 0;
+
+    friend class Serializer;
+
+    void save(Serializer& rSerializer) const override
+    {
+        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, GeoThermalLaw)
+        rSerializer.save("NumberOfDimensions", mNumberOfDimensions);
+    }
+
+    void load(Serializer& rSerializer) override
+    {
+        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, GeoThermalLaw)
+        rSerializer.load("NumberOfDimensions", mNumberOfDimensions);
+    }
+};
+} // namespace Kratos.

applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h

@@ -0,0 +1,48 @@
+// 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 GeoThermalDispersionLaw
+ * @ingroup GeoMechanicsApplication
+ * @brief This class defines the thermal dispersion for heat cases
+ */
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalLaw
+{
+public:
+    /// Counted pointer of GeoThermalLaw
+    KRATOS_CLASS_POINTER_DEFINITION(GeoThermalLaw);
+
+    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 {}
+
+    virtual void load(Serializer& rSerializer) {}
+};
+} // 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"
@@ -269,9 +270,25 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Ele
         const ProcessInfo&                               rCurrentProcessInfo) const
     {
         const std::size_t number_of_dimensions = GetGeometry().LocalSpaceDimension();
+        std::unique_ptr<GeoThermalLaw> law;
+
+        if (GetProperties().Has(THERMAL_LAW)) {
+            const std::string& rThermalLawName = GetProperties()[THERMAL_LAW];
+            if (rThermalLawName == "GeoThermalDispersionLaw") {
+                law = std::make_unique<GeoThermalDispersionLaw>(number_of_dimensions);
+            } 
+            else if (rThermalLawName == "GeoThermalFilterLaw") {
+                law = std::make_unique<GeoThermalFilterLaw>(number_of_dimensions);
+            } 
+            else {
+                KRATOS_ERROR << "Undefined THERMAL_LAW! " << rThermalLawName << std::endl;
+            }
+        } 
+        else {
+            law = std::make_unique<GeoThermalDispersionLaw>(number_of_dimensions);
+        }
 
-        GeoThermalDispersionLaw law{number_of_dimensions};
-        const auto constitutive_matrix = law.CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo);
+        const auto constitutive_matrix = law->CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo);
 
         auto result = BoundedMatrix<double, TNumNodes, TNumNodes>{ZeroMatrix{TNumNodes, TNumNodes}};
         for (unsigned int integration_point_index = 0;

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)
 
     // 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 )
 
 // 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 )
 
     // 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(GeoThermalFilterLaw law{0}, "Got invalid number of dimensions: 0")
 }
 
 } // namespace Kratos::Testing

applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_filter_law.cpp

@@ -0,0 +1,48 @@
+// 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 cond_prop = r_model_part.CreateNewProperties(0);
+    cond_prop->SetValue(THERMAL_CONDUCTIVITY_WATER, 1000.0);
+
+    const SizeType      dimension = 1;
+    GeoThermalFilterLaw geo_thermal_filter_law(dimension);
+    ProcessInfo         info;
+
+    const Matrix thermal_filter_matrix = geo_thermal_filter_law.CalculateThermalDispersionMatrix(*cond_prop, info);
+
+    Matrix expected_solution = ScalarMatrix(1, 1, 1000.0);
+
+    constexpr double tolerance{1.0e-6};
+
+    KRATOS_EXPECT_NEAR(thermal_filter_matrix(0, 0), expected_solution(0, 0), tolerance);
+}
+
+KRATOS_TEST_CASE_IN_SUITE(TestFilterLawThrowsWhenDimensionInvalid, KratosGeoMechanicsFastSuite)
+{
+    KRATOS_EXPECT_EXCEPTION_IS_THROWN(GeoThermalFilterLaw law{2},
+        "Got invalid number of dimensions. The dimension has to be 1, but got: 2")
+}
+
+} // namespace Kratos::Testing