Geo/14021 interface fluid body flow calculator

applications/GeoMechanicsApplication/custom_elements/Pw_element.cpp

@@ -363,11 +363,6 @@ template <unsigned int TDim, unsigned int TNumNodes>
 std::vector<Vector> PwElement<TDim, TNumNodes>::CalculateProjectedGravityAtIntegrationPoints(const Matrix& rNContainer) const
 {
     const auto number_integration_points = GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
-    GeometryType::JacobiansType J_container{number_integration_points};
-    for (auto& j : J_container) {
-        j.resize(GetGeometry().WorkingSpaceDimension(), GetGeometry().LocalSpaceDimension(), false);
-    }
-    GetGeometry().Jacobian(J_container, this->GetIntegrationMethod());
 
     array_1d<double, TNumNodes * TDim> volume_acceleration;
     GeoElementUtilities::GetNodalVariableVector<TDim, TNumNodes>(volume_acceleration, GetGeometry(),
@@ -377,6 +372,11 @@ std::vector<Vector> PwElement<TDim, TNumNodes>::CalculateProjectedGravityAtInteg
     std::vector<Vector> projected_gravity;
     projected_gravity.reserve(number_integration_points);
     if (GetGeometry().LocalSpaceDimension() == 1) {
+        GeometryType::JacobiansType J_container{number_integration_points};
+        for (auto& j : J_container) {
+            j.resize(GetGeometry().WorkingSpaceDimension(), GetGeometry().LocalSpaceDimension(), false);
+        }
+        GetGeometry().Jacobian(J_container, this->GetIntegrationMethod());
         for (unsigned int integration_point_index = 0;
              integration_point_index < number_integration_points; ++integration_point_index) {
             GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(

applications/GeoMechanicsApplication/custom_elements/U_Pw_interface_element.h

@@ -14,6 +14,8 @@
 #pragma once
 
 #include "contribution_calculators/calculation_contribution.h"
+#include "contribution_calculators/fluid_body_flow_calculator.hpp"
+#include "contribution_calculators/permeability_calculator.hpp"
 #include "contribution_calculators/stiffness_calculator.hpp"
 #include "custom_elements/contribution_calculators/pu_coupling_calculator.hpp"
 #include "custom_elements/contribution_calculators/up_coupling_calculator.hpp"
@@ -30,8 +32,6 @@
 
 namespace Kratos
 {
-class RetentionLaw;
-
 class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwInterfaceElement : public Element
 {
 public:
@@ -65,13 +65,16 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwInterfaceElement : public Element
 
     void EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo&) const override;
     void CalculateAndAssignStiffnessMatrix(Element::MatrixType& rLeftHandSideMatrix, const ProcessInfo& rProcessInfo);
+    void CalculateAndAssignPermeabilityMatrix(Element::MatrixType& rLeftHandSideMatrix);
     void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo&) override;
     void CalculateAndAssembleStiffnessForceVector(Element::VectorType& rRightHandSideVector,
                                                   const ProcessInfo&   rProcessInfo);
     void CalculateAndAssignUPCouplingMatrix(MatrixType& rLeftHandSideMatrix) const;
     void CalculateAndAssignPUCouplingMatrix(MatrixType& rLeftHandSideMatrix) const;
     void CalculateAndAssembleUPCouplingForceVector(Element::VectorType& rRightHandSideVector) const;
     void CalculateAndAssemblePUCouplingForceVector(Element::VectorType& rRightHandSideVector) const;
+    void CalculateAndAssemblePermeabilityFlowVector(Element::VectorType& rRightHandSideVector);
+    void CalculateAndAssembleFluidBodyFlowVector(Element::VectorType& rRightHandSideVector);
     void CalculateRightHandSide(VectorType& rRightHandSideVector, const ProcessInfo&) override;
     void CalculateLocalSystem(MatrixType&        rLeftHandSideMatrix,
                               VectorType&        rRightHandSideVector,
@@ -117,6 +120,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwInterfaceElement : public Element
     Matrix RotateStressToLocalCoordinates(const Geo::IntegrationPointType& rIntegrationPoint,
                                           const Vector& rGlobalStressVector) const;
     Vector ConvertLocalStressToTraction(const Matrix& rLocalStress) const;
+    Matrix CalculatePwBMatrix(const Geo::IntegrationPointType& rIntegrationPoint) const;
+    Geometry<Node>::ShapeFunctionsGradientsType CalculateLocalPwBMatricesAtIntegrationPoints() const;
+    std::vector<double> CalculateIntegrationPointFluidPressures() const;
+    std::vector<Vector> CalculateProjectedGravity() const;
 
     Geo::BMatricesGetter                 CreateBMatricesGetter() const;
     Geo::StrainVectorsGetter             CreateRelativeDisplacementsGetter() const;
@@ -127,29 +134,34 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwInterfaceElement : public Element
     std::function<Vector()>              CreateVoigtVectorGetter() const;
     std::function<std::vector<double>()> CreateBiotCoefficientsGetter() const;
     std::function<std::vector<double>()> CreateBishopCoefficientsGetter() const;
-    std::function<Vector()>              CreateNodalPressuresGetter() const;
+    Geo::NodalValuesGetter               CreateWaterPressureGeometryNodalVariableGetter() const;
     std::function<Vector()>              CreateNodalVelocitiesGetter() const;
     std::function<std::vector<double>()> CreateDegreesOfSaturationGetter() const;
 
-    std::vector<double> CalculateBiotCoefficients() const;
-    std::vector<double> CalculateBishopCoefficients() const;
-    Vector              CalculateIntegrationPointFluidPressures() const;
-    Vector              GetWaterPressureGeometryNodalVariable() const;
-    Vector              GetGeometryVelocityValues() const;
-    Matrix              GetNpContainer() const;
-    std::vector<double> GetDegreesOfSaturationValues() const;
+    std::vector<double>                   CalculateBiotCoefficients() const;
+    std::vector<double>                   CalculateBishopCoefficients() const;
+    Vector                                GetWaterPressureGeometryNodalVariable() const;
+    Vector                                GetGeometryVelocityValues() const;
+    Matrix                                GetNpContainer() const;
+    std::vector<double>                   GetDegreesOfSaturationValues() const;
+    Geo::RetentionLawsGetter              CreateRetentionLawsGetter() const;
+    Geo::MaterialPermeabilityMatrixGetter CreateMaterialPermeabilityGetter() const;
+    Geo::ShapeFunctionGradientsGetter     CreatePwBMatricesGetter() const;
+    Geo::IntegrationPointValuesGetter     CreateFluidPressureCalculator() const;
+    std::function<std::vector<Vector>()>  CreateProjectedGravityCalculator() const;
 
     template <unsigned int MatrixSize>
-    typename StiffnessCalculator<MatrixSize>::InputProvider CreateStiffnessInputProvider(const ProcessInfo& rProcessInfo);
+    typename StiffnessCalculator<MatrixSize>::InputProvider CreateStiffnessInputProvider(const ProcessInfo& rProcessInfo) const;
 
     template <unsigned int MatrixSize>
-    auto CreateStiffnessCalculator(const ProcessInfo& rProcessInfo);
+    auto CreateStiffnessCalculator(const ProcessInfo& rProcessInfo) const;
 
     template <unsigned int MatrixSize>
-    void CalculateAndAssignStiffnessMatrix(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rProcessInfo);
+    void CalculateAndAssignStiffnessMatrix(MatrixType& rLeftHandSideMatrix, const ProcessInfo& rProcessInfo) const;
 
     template <unsigned int MatrixSize>
-    void CalculateAndAssembleStiffnessForceVector(VectorType& rRightHandSideVector, const ProcessInfo& rProcessInfo);
+    void CalculateAndAssembleStiffnessForceVector(VectorType&        rRightHandSideVector,
+                                                  const ProcessInfo& rProcessInfo) const;
 
     template <unsigned int NumberOfRows, unsigned int NumberOfColumns>
     typename UPCouplingCalculator<NumberOfRows, NumberOfColumns>::InputProvider CreateUPCouplingInputProvider() const;
@@ -175,15 +187,36 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwInterfaceElement : public Element
     template <unsigned int NumberOfRows, unsigned int NumberOfColumns>
     void CalculateAndAssemblePUCouplingForceVector(VectorType& rRightHandSideVector) const;
 
+    template <unsigned int TNumNodes>
+    typename PermeabilityCalculator<TNumNodes>::InputProvider CreatePermeabilityInputProvider() const;
+
+    template <unsigned int TNumNodes>
+    auto CreatePermeabilityCalculator() const;
+
+    template <unsigned int TnumNodes>
+    void CalculateAndAssignPermeabilityMatrix(MatrixType& rLeftHandSideMatrix) const;
+
+    template <unsigned int TnumNodes>
+    void CalculateAndAssemblePermeabilityFlowVector(VectorType& rRightHandSideVector) const;
+
+    template <unsigned int TNumNodes>
+    typename FluidBodyFlowCalculator<TNumNodes>::InputProvider CreateFluidBodyFlowInputProvider() const;
+
+    template <unsigned int TNumNodes>
+    auto CreateFluidBodyFlowCalculator() const;
+
+    template <unsigned int TnumNodes>
+    void CalculateAndAssembleFluidBodyFlowVector(VectorType& rRightHandSideVector) const;
+
     std::function<Matrix(const Geometry<Node>&, const array_1d<double, 3>&)> mfpCalculateRotationMatrix;
 
     std::unique_ptr<IntegrationScheme>    mpIntegrationScheme;
     std::unique_ptr<StressStatePolicy>    mpStressStatePolicy;
     std::vector<ConstitutiveLaw::Pointer> mConstitutiveLaws;
+    std::vector<RetentionLaw::Pointer>    mRetentionLaws;
     IntegrationCoefficientsCalculator     mIntegrationCoefficientsCalculator;
     Geo::GeometryUniquePtr                mpOptionalPressureGeometry;
     std::vector<CalculationContribution>  mContributions;
-    std::vector<RetentionLaw::Pointer>    mRetentionLawVector;
 
     friend class Serializer;
 };

applications/GeoMechanicsApplication/custom_elements/contribution_calculators/permeability_calculator.hpp

@@ -22,7 +22,7 @@ namespace Kratos
 {
 
 template <unsigned int TNumNodes>
-class PermeabilityCalculator : public ContributionCalculator<TNumNodes> ///
+class PermeabilityCalculator : public ContributionCalculator<TNumNodes>
 {
 public:
     struct InputProvider {
@@ -95,7 +95,6 @@ class PermeabilityCalculator : public ContributionCalculator<TNumNodes> ///
             const auto relative_permeability =
                 mInputProvider.GetRetentionLaws()[integration_point_index]->CalculateRelativePermeability(
                     retention_parameters);
-
             noalias(result) += GeoTransportEquationUtilities::CalculatePermeabilityMatrix(
                 r_shape_function_gradients[integration_point_index], dynamic_viscosity_inverse, material_permeability,
                 relative_permeability, r_integration_coefficients[integration_point_index]);

applications/GeoMechanicsApplication/custom_elements/contribution_calculators/pu_coupling_calculator.hpp

@@ -81,11 +81,11 @@ class PUCouplingCalculator : public ContributionCalculator<NumberOfRows, NumberO
         // Also note that instead of the Bishop coefficients, the degrees of saturation are passed.
         // Finally, the nodal water pressure getter is not needed for the matrix calculation, so
         // we pass a dummy function.
-        auto dummy_nodal_water_pressure_function = []() { return Vector(); };
+        auto dummy_nodal_values_getter = [](const Variable<double>&) { return Vector(); };
         typename UPCouplingCalculator<NumberOfColumns, NumberOfRows>::InputProvider input_provider(
             mInputProvider.GetNpContainer, mInputProvider.GetBMatrices, mInputProvider.GetVoigtVector,
             mInputProvider.GetIntegrationCoefficients, mInputProvider.GetBiotCoefficients,
-            mInputProvider.GetDegreesOfSaturation, dummy_nodal_water_pressure_function);
+            mInputProvider.GetDegreesOfSaturation, dummy_nodal_values_getter);
         UPCouplingCalculator<NumberOfColumns, NumberOfRows> up_coupling_calculator(input_provider);
         const auto up_coupling_matrix = up_coupling_calculator.LHSContribution().value();
 

applications/GeoMechanicsApplication/custom_elements/contribution_calculators/up_coupling_calculator.hpp

@@ -23,30 +23,30 @@ class UPCouplingCalculator : public ContributionCalculator<NumberOfRows, NumberO
 {
 public:
     struct InputProvider {
-        InputProvider(std::function<const Matrix()>        GetNpContainer,
-                      Geo::BMatricesGetter                 GetBMatrices,
-                      std::function<Vector()>              GetVoigtVector,
-                      Geo::IntegrationCoefficientsGetter   GetIntegrationCoefficients,
-                      std::function<std::vector<double>()> GetBiotCoefficients,
-                      std::function<std::vector<double>()> GetBishopCoefficients,
-                      std::function<Vector()>              GetNodalWaterPressures)
+        InputProvider(std::function<const Matrix()>                  GetNpContainer,
+                      Geo::BMatricesGetter                           GetBMatrices,
+                      std::function<Vector()>                        GetVoigtVector,
+                      Geo::IntegrationCoefficientsGetter             GetIntegrationCoefficients,
+                      std::function<std::vector<double>()>           GetBiotCoefficients,
+                      std::function<std::vector<double>()>           GetBishopCoefficients,
+                      std::function<Vector(const Variable<double>&)> GetNodalValues)
             : GetNpContainer(std::move(GetNpContainer)),
               GetBMatrices(std::move(GetBMatrices)),
               GetVoigtVector(std::move(GetVoigtVector)),
               GetIntegrationCoefficients(std::move(GetIntegrationCoefficients)),
               GetBiotCoefficients(std::move(GetBiotCoefficients)),
               GetBishopCoefficients(std::move(GetBishopCoefficients)),
-              GetNodalWaterPressures(std::move(GetNodalWaterPressures))
+              GetNodalValues(std::move(GetNodalValues))
         {
         }
 
-        std::function<const Matrix()>        GetNpContainer;
-        Geo::BMatricesGetter                 GetBMatrices;
-        std::function<Vector()>              GetVoigtVector;
-        Geo::IntegrationCoefficientsGetter   GetIntegrationCoefficients;
-        std::function<std::vector<double>()> GetBiotCoefficients;
-        std::function<std::vector<double>()> GetBishopCoefficients;
-        std::function<Vector()>              GetNodalWaterPressures;
+        std::function<const Matrix()>                  GetNpContainer;
+        Geo::BMatricesGetter                           GetBMatrices;
+        std::function<Vector()>                        GetVoigtVector;
+        Geo::IntegrationCoefficientsGetter             GetIntegrationCoefficients;
+        std::function<std::vector<double>()>           GetBiotCoefficients;
+        std::function<std::vector<double>()>           GetBishopCoefficients;
+        std::function<Vector(const Variable<double>&)> GetNodalValues;
     };
 
     explicit UPCouplingCalculator(InputProvider CouplingInputProvider)
@@ -63,7 +63,7 @@ class UPCouplingCalculator : public ContributionCalculator<NumberOfRows, NumberO
 
     typename BaseType::RHSVectorType RHSContribution() override
     {
-        return prod(CalculateCouplingMatrix(), mInputProvider.GetNodalWaterPressures());
+        return prod(CalculateCouplingMatrix(), mInputProvider.GetNodalValues(WATER_PRESSURE));
     }
 
     std::pair<std::optional<typename BaseType::LHSMatrixType>, typename BaseType::RHSVectorType> LocalSystemContribution() override

applications/GeoMechanicsApplication/custom_utilities/element_utilities.cpp

@@ -96,6 +96,14 @@ Matrix GeoElementUtilities::FillPermeabilityMatrix(const Element::PropertiesType
     }
 }
 
+Matrix GeoElementUtilities::FillInterfacePermeabilityMatrix(const Element::PropertiesType& rProperties,
+                                                            std::size_t Dimension)
+{
+    auto result  = Matrix{Dimension, Dimension, 0.0};
+    result(0, 0) = rProperties[TRANSVERSAL_PERMEABILITY];
+    return result;
+}
+
 double GeoElementUtilities::CalculateRadius(const Vector& rN, const GeometryType& rGeometry)
 {
     auto radius = 0.0;

applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp

@@ -128,6 +128,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoElementUtilities
 
     static Matrix FillPermeabilityMatrix(const Element::PropertiesType& rProperties, std::size_t Dimension);
 
+    static Matrix FillInterfacePermeabilityMatrix(const Element::PropertiesType& rProperties, std::size_t Dimension);
+
     template <typename MatrixType1, typename MatrixType2>
     static void AssembleUUBlockMatrix(MatrixType1& rLeftHandSideMatrix, const MatrixType2& rUUBlockMatrix)
     {

applications/GeoMechanicsApplication/custom_utilities/geometry_utilities.cpp

@@ -12,8 +12,7 @@
 #include "geometry_utilities.h"
 #include "geometries/geometry_data.h"
 #include "math_utilities.hpp"
-
-#include <boost/numeric/ublas/assignment.hpp>
+#include "ublas_utilities.h"
 
 namespace
 {
@@ -96,23 +95,17 @@ void ReverseNodes(InputIt                               Begin,
 
 namespace Kratos
 {
-
 Matrix GeometryUtilities::Calculate2DRotationMatrixForLineGeometry(const Geometry<Node>& rGeometry,
                                                                    const array_1d<double, 3>& rLocalCoordinate)
 {
     // Since the shape functions depend on one coordinate only
-    // for lines, the jacobian only has one column.
+    // for lines, the Jacobian only has one column.
     Matrix jacobian;
     rGeometry.Jacobian(jacobian, rLocalCoordinate);
     const auto tangential_vector = GeoMechanicsMathUtilities::Normalized(Vector{column(jacobian, 0)});
 
-    // clang-format off
-    Matrix result(2, 2);
-    result <<= tangential_vector[0], -tangential_vector[1],
-               tangential_vector[1],  tangential_vector[0];
-    // clang-format on
-
-    return result;
+    return UblasUtilities::CreateMatrix({{tangential_vector[0], -tangential_vector[1]},
+                                         {tangential_vector[1], tangential_vector[0]}});
 }
 
 Matrix GeometryUtilities::Calculate3DRotationMatrixForPlaneGeometry(const Geometry<Node>& rGeometry,
@@ -128,14 +121,10 @@ Matrix GeometryUtilities::Calculate3DRotationMatrixForPlaneGeometry(const Geomet
     MathUtils<>::CrossProduct(normal_vector, tangential_vector_1, tangential_vector_2);
     normal_vector = GeoMechanicsMathUtilities::Normalized(normal_vector);
 
-    // clang-format off
-    Matrix result(3, 3);
-    result <<= tangential_vector_1[0], tangential_vector_2[0], normal_vector[0],
-               tangential_vector_1[1], tangential_vector_2[1], normal_vector[1],
-               tangential_vector_1[2], tangential_vector_2[2], normal_vector[2];
-    // clang-format on
-
-    return result;
+    return UblasUtilities::CreateMatrix(
+        {{tangential_vector_1[0], tangential_vector_2[0], normal_vector[0]},
+         {tangential_vector_1[1], tangential_vector_2[1], normal_vector[1]},
+         {tangential_vector_1[2], tangential_vector_2[2], normal_vector[2]}});
 }
 
 void GeometryUtilities::ReverseNodes(PointerVector<Node>&                  rNodes,

applications/GeoMechanicsApplication/geo_aliases.h

@@ -43,11 +43,14 @@ using KappaDependentFunction = std::function<double(double)>;
 using GeometryUniquePtr = std::unique_ptr<Geometry<Node>>;
 
 using BMatricesGetter               = std::function<std::vector<Matrix>()>;
+using ShapeFunctionGradientsGetter  = std::function<Geometry<Node>::ShapeFunctionsGradientsType()>;
 using StrainVectorsGetter           = std::function<std::vector<Vector>()>;
 using IntegrationCoefficientsGetter = std::function<std::vector<double>()>;
 using PropertiesGetter              = std::function<const Properties&()>;
 using ProcessInfoGetter             = std::function<const ProcessInfo&()>;
 using ConstitutiveLawsGetter        = std::function<const std::vector<ConstitutiveLaw::Pointer>&()>;
 using RetentionLawsGetter           = std::function<const std::vector<RetentionLaw::Pointer>&()>;
 using MaterialPermeabilityMatrixGetter = std::function<Matrix()>;
+using NodalValuesGetter                = std::function<Vector(const Variable<double>&)>;
+using IntegrationPointValuesGetter     = std::function<std::vector<double>()>;
 } // namespace Kratos::Geo

applications/GeoMechanicsApplication/tests/cpp_tests/custom_elements/contribution_calculators/test_coupling_calculators.cpp

@@ -47,8 +47,8 @@ typename UPCouplingCalculator<NumberOfRows, NumberOfColumns>::InputProvider Crea
     auto get_bishop_coefficients = [BishopCoefficient, NumberOfIntegrationPoints]() {
         return std::vector(NumberOfIntegrationPoints, BishopCoefficient);
     };
-    auto get_voigt_vector    = [rVoigtVector]() { return rVoigtVector; };
-    auto get_fluid_pressures = [rFluidPressures]() { return rFluidPressures; };
+    auto get_voigt_vector = [rVoigtVector]() { return rVoigtVector; };
+    auto get_fluid_pressures = [rFluidPressures](const Variable<double>&) { return rFluidPressures; };
 
     return typename UPCouplingCalculator<NumberOfRows, NumberOfColumns>::InputProvider(
         get_np_container, get_b_matrices, get_voigt_vector, get_integration_coefficients,