[GeoMechanicsApplication] Extend class LineInterfaceElement such that it also includes planar interface elements (#13634)

* renamed LineInterfaceElement as InterfaceElement

* added LocalSpaceDimension check and Calculate2DRotationMatrixForPlaneGeometry

* updated element_hierarchy.svg

* added unit tests and 8 still red

* added LumpedIntegrationScheme

* changed 2D triangle into 3D triangle

* column based rotation matrix

* added pStressStatePolicy in constructors

* fixed unit tests

* extracted MakeIntegrationSheme

* removed code smells

* joint CreateLinearElasticMaterialProperties and CreatePlanarElasticMaterialProperties

* joint Create2DInterfaceElementWithUDofs and Create3DInterfaceElementWithUDofs

* introduced serializer

* extracted CreateAndInitializeElement

* used emplace_back

* added static_assert

* added function pointer mfpCalculateRotationMatrix

* response to the review 1

* response to review 2

* created CreateExpectedLeftHandSideForTriangleElement

* a bit nicer

* one more use of CreateExpectedLeftHandSideForTriangleElement

* reverted expected_left_hand_side vlaues

* renamed as ExpectedLeftHandSideForTriangleElement

Author: markelov208

applications/GeoMechanicsApplication/custom_elements/element_hierarchy.puml

@@ -47,8 +47,8 @@ TransientPwInterfaceElement <|-- SteadyStatePwInterfaceElement
 SteadyStatePwInterfaceElement <|-- SteadyStatePwPipingElement
 UPwSmallStrainInterfaceElement <|-- UPwSmallStrainLinkInterfaceElement
 
-class LineInterfaceElement
-Element <|-- LineInterfaceElement
+class InterfaceElement
+Element <|-- InterfaceElement
 
 'Pw Elements
 class TransientPwLineElement

applications/GeoMechanicsApplication/custom_elements/element_hierarchy.svg

@@ -8,15 +8,17 @@
 //  License:         geo_mechanics_application/license.txt
 //
 //  Main authors:    Richard Faasse
+//                   Gennady Markelov
 //
-#include "line_interface_element.h"
+#include "interface_element.h"
 
 #include "custom_utilities/dof_utilities.h"
 #include "custom_utilities/element_utilities.hpp"
 #include "custom_utilities/equation_of_motion_utilities.h"
 #include "custom_utilities/geometry_utilities.h"
 #include "interface_stress_state.h"
 #include "lobatto_integration_scheme.h"
+#include "lumped_integration_scheme.h"
 
 namespace
 {
@@ -60,42 +62,54 @@ std::vector<Matrix> CalculateConstitutiveMatricesAtIntegrationPoints(
 namespace Kratos
 {
 
-LineInterfaceElement::LineInterfaceElement(IndexType NewId,
-                                           const Geometry<GeometricalObject::NodeType>::Pointer& rGeometry,
-                                           const Properties::Pointer& rProperties)
-    : Element(NewId, rGeometry, rProperties),
-      mIntegrationScheme(std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2)),
-      mStressStatePolicy(std::make_unique<Line2DInterfaceStressState>())
+InterfaceElement::InterfaceElement(IndexType                                             NewId,
+                                   const Geometry<GeometricalObject::NodeType>::Pointer& rGeometry,
+                                   const Properties::Pointer&         rProperties,
+                                   std::unique_ptr<StressStatePolicy> pStressStatePolicy)
+    : Element(NewId, rGeometry, rProperties), mpStressStatePolicy(std::move(pStressStatePolicy))
 {
+    MakeIntegrationSchemeAndAssignFunction();
 }
 
-LineInterfaceElement::LineInterfaceElement(IndexType NewId, const GeometryType::Pointer& rGeometry)
-    : Element(NewId, rGeometry),
-      mIntegrationScheme(std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2)),
-      mStressStatePolicy(std::make_unique<Line2DInterfaceStressState>())
+InterfaceElement::InterfaceElement(IndexType                          NewId,
+                                   const GeometryType::Pointer&       rGeometry,
+                                   std::unique_ptr<StressStatePolicy> pStressStatePolicy)
+    : Element(NewId, rGeometry), mpStressStatePolicy(std::move(pStressStatePolicy))
 {
+    MakeIntegrationSchemeAndAssignFunction();
 }
 
-Element::Pointer LineInterfaceElement::Create(IndexType               NewId,
-                                              const NodesArrayType&   rNodes,
-                                              PropertiesType::Pointer pProperties) const
+void InterfaceElement::MakeIntegrationSchemeAndAssignFunction()
+{
+    if (GetGeometry().LocalSpaceDimension() == 1) {
+        mIntegrationScheme = std::make_unique<LobattoIntegrationScheme>(GetGeometry().PointsNumber() / 2);
+        mfpCalculateRotationMatrix = GeometryUtilities::Calculate2DRotationMatrixForLineGeometry;
+    } else {
+        mIntegrationScheme = std::make_unique<LumpedIntegrationScheme>(GetGeometry().PointsNumber() / 2);
+        mfpCalculateRotationMatrix = GeometryUtilities::Calculate3DRotationMatrixForPlaneGeometry;
+    }
+}
+
+Element::Pointer InterfaceElement::Create(IndexType               NewId,
+                                          const NodesArrayType&   rNodes,
+                                          PropertiesType::Pointer pProperties) const
 {
     return Create(NewId, this->GetGeometry().Create(rNodes), pProperties);
 }
 
-Element::Pointer LineInterfaceElement::Create(IndexType               NewId,
-                                              GeometryType::Pointer   pGeometry,
-                                              PropertiesType::Pointer pProperties) const
+Element::Pointer InterfaceElement::Create(IndexType               NewId,
+                                          GeometryType::Pointer   pGeometry,
+                                          PropertiesType::Pointer pProperties) const
 {
-    return make_intrusive<LineInterfaceElement>(NewId, pGeometry, pProperties);
+    return make_intrusive<InterfaceElement>(NewId, pGeometry, pProperties, mpStressStatePolicy->Clone());
 }
 
-void LineInterfaceElement::EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo&) const
+void InterfaceElement::EquationIdVector(EquationIdVectorType& rResult, const ProcessInfo&) const
 {
     rResult = Geo::DofUtilities::ExtractEquationIdsFrom(GetDofs());
 }
 
-void LineInterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo&)
+void InterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, const ProcessInfo&)
 {
     // Currently, the left-hand side matrix only includes the stiffness matrix. In the future, it
     // will also include water pressure contributions and coupling terms.
@@ -104,7 +118,7 @@ void LineInterfaceElement::CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix
         CalculateConstitutiveMatricesAtIntegrationPoints(), CalculateIntegrationCoefficients());
 }
 
-void LineInterfaceElement::CalculateRightHandSide(Element::VectorType& rRightHandSideVector, const ProcessInfo&)
+void InterfaceElement::CalculateRightHandSide(Element::VectorType& rRightHandSideVector, const ProcessInfo&)
 {
     // Currently, the right-hand side only includes the internal force vector. In the future, it
     // will also include water pressure contributions and coupling terms.
@@ -116,17 +130,17 @@ void LineInterfaceElement::CalculateRightHandSide(Element::VectorType& rRightHan
         local_b_matrices, tractions, integration_coefficients);
 }
 
-void LineInterfaceElement::CalculateLocalSystem(MatrixType&        rLeftHandSideMatrix,
-                                                VectorType&        rRightHandSideVector,
-                                                const ProcessInfo& rCurrentProcessInfo)
+void InterfaceElement::CalculateLocalSystem(MatrixType&        rLeftHandSideMatrix,
+                                            VectorType&        rRightHandSideVector,
+                                            const ProcessInfo& rCurrentProcessInfo)
 {
     CalculateLeftHandSide(rLeftHandSideMatrix, rCurrentProcessInfo);
     CalculateRightHandSide(rRightHandSideVector, rCurrentProcessInfo);
 }
 
-void LineInterfaceElement::CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
-                                                        std::vector<Vector>&    rOutput,
-                                                        const ProcessInfo&      rCurrentProcessInfo)
+void InterfaceElement::CalculateOnIntegrationPoints(const Variable<Vector>& rVariable,
+                                                    std::vector<Vector>&    rOutput,
+                                                    const ProcessInfo&      rCurrentProcessInfo)
 {
     if (rVariable == STRAIN) {
         const auto local_b_matrices = CalculateLocalBMatricesAtIntegrationPoints();
@@ -140,22 +154,22 @@ void LineInterfaceElement::CalculateOnIntegrationPoints(const Variable<Vector>&
     }
 }
 
-void LineInterfaceElement::CalculateOnIntegrationPoints(const Variable<ConstitutiveLaw::Pointer>& rVariable,
-                                                        std::vector<ConstitutiveLaw::Pointer>& rOutput,
-                                                        const ProcessInfo&)
+void InterfaceElement::CalculateOnIntegrationPoints(const Variable<ConstitutiveLaw::Pointer>& rVariable,
+                                                    std::vector<ConstitutiveLaw::Pointer>& rOutput,
+                                                    const ProcessInfo&)
 {
     KRATOS_ERROR_IF_NOT(rVariable == CONSTITUTIVE_LAW)
         << "Cannot calculate on integration points: got unexpected variable " << rVariable.Name() << "\n";
 
     rOutput = mConstitutiveLaws;
 }
 
-void LineInterfaceElement::GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo&) const
+void InterfaceElement::GetDofList(DofsVectorType& rElementalDofList, const ProcessInfo&) const
 {
     rElementalDofList = GetDofs();
 }
 
-void LineInterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
+void InterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
 {
     Element::Initialize(rCurrentProcessInfo);
 
@@ -170,7 +184,7 @@ void LineInterfaceElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
     }
 }
 
-int LineInterfaceElement::Check(const ProcessInfo& rCurrentProcessInfo) const
+int InterfaceElement::Check(const ProcessInfo& rCurrentProcessInfo) const
 {
     int error = Element::Check(rCurrentProcessInfo);
     if (error != 0) return error;
@@ -189,30 +203,29 @@ int LineInterfaceElement::Check(const ProcessInfo& rCurrentProcessInfo) const
     return 0;
 }
 
-Element::DofsVectorType LineInterfaceElement::GetDofs() const
+Element::DofsVectorType InterfaceElement::GetDofs() const
 {
     const auto no_Pw_geometry_yet = Geometry<Node>{};
     return Geo::DofUtilities::ExtractUPwDofsFromNodes(GetGeometry(), no_Pw_geometry_yet,
                                                       GetGeometry().WorkingSpaceDimension());
 }
 
-std::vector<Matrix> LineInterfaceElement::CalculateLocalBMatricesAtIntegrationPoints() const
+std::vector<Matrix> InterfaceElement::CalculateLocalBMatricesAtIntegrationPoints() const
 {
     const auto& r_integration_points = mIntegrationScheme->GetIntegrationPoints();
     const auto  shape_function_values_at_integration_points =
         GeoElementUtilities::EvaluateShapeFunctionsAtIntegrationPoints(r_integration_points, GetGeometry());
 
     auto result = std::vector<Matrix>{};
     result.reserve(shape_function_values_at_integration_points.size());
-    auto calculate_local_b_matrix = [&r_geometry = GetGeometry(), p_policy = mStressStatePolicy.get()](
-                                        const auto& rShapeFunctionValuesAtIntegrationPoint,
-                                        const auto& rIntegrationPoint) {
+    auto calculate_local_b_matrix = [&r_geometry = GetGeometry(), this](const auto& rShapeFunctionValuesAtIntegrationPoint,
+                                                                        const auto& rIntegrationPoint) {
         // For interface elements, the shape function gradients are not used, since these are
         // non-continuum elements. Therefore, we pass an empty matrix.
+        const auto rotation_matrix = mfpCalculateRotationMatrix(r_geometry, rIntegrationPoint);
         const auto dummy_gradients = Matrix{};
-        return Matrix{prod(
-            GeometryUtilities::Calculate2DRotationMatrixForLineGeometry(r_geometry, rIntegrationPoint),
-            p_policy->CalculateBMatrix(dummy_gradients, rShapeFunctionValuesAtIntegrationPoint, r_geometry))};
+        return Matrix{prod(rotation_matrix, mpStressStatePolicy->CalculateBMatrix(dummy_gradients, rShapeFunctionValuesAtIntegrationPoint,
+                                                                                  r_geometry))};
     };
     std::transform(shape_function_values_at_integration_points.begin(),
                    shape_function_values_at_integration_points.end(), r_integration_points.begin(),
@@ -221,21 +234,20 @@ std::vector<Matrix> LineInterfaceElement::CalculateLocalBMatricesAtIntegrationPo
     return result;
 }
 
-std::vector<double> LineInterfaceElement::CalculateIntegrationCoefficients() const
+std::vector<double> InterfaceElement::CalculateIntegrationCoefficients() const
 {
     const auto determinants_of_jacobian = CalculateDeterminantsOfJacobiansAtIntegrationPoints(
         mIntegrationScheme->GetIntegrationPoints(), GetGeometry());
     return mIntegrationCoefficientsCalculator.Run<>(mIntegrationScheme->GetIntegrationPoints(),
                                                     determinants_of_jacobian, this);
 }
 
-std::vector<Matrix> LineInterfaceElement::CalculateConstitutiveMatricesAtIntegrationPoints()
+std::vector<Matrix> InterfaceElement::CalculateConstitutiveMatricesAtIntegrationPoints()
 {
     return ::CalculateConstitutiveMatricesAtIntegrationPoints(mConstitutiveLaws, GetProperties());
 }
 
-std::vector<Vector> LineInterfaceElement::CalculateRelativeDisplacementsAtIntegrationPoints(
-    const std::vector<Matrix>& rLocalBMatrices) const
+std::vector<Vector> InterfaceElement::CalculateRelativeDisplacementsAtIntegrationPoints(const std::vector<Matrix>& rLocalBMatrices) const
 {
     const Geometry<Node> no_Pw_geometry;
     const auto           dofs = Geo::DofUtilities::ExtractUPwDofsFromNodes(
@@ -255,7 +267,7 @@ std::vector<Vector> LineInterfaceElement::CalculateRelativeDisplacementsAtIntegr
     return result;
 }
 
-std::vector<Vector> LineInterfaceElement::CalculateTractionsAtIntegrationPoints(const std::vector<Vector>& rRelativeDisplacements)
+std::vector<Vector> InterfaceElement::CalculateTractionsAtIntegrationPoints(const std::vector<Vector>& rRelativeDisplacements)
 {
     // We have to make a copy of each relative displacement vector, since setting it at the
     // constitutive law parameters requires a reference to a _mutable_ object!
@@ -276,4 +288,9 @@ std::vector<Vector> LineInterfaceElement::CalculateTractionsAtIntegrationPoints(
     return result;
 }
 
+// Instances of this class can not be copied but can be moved. Check that at compile time.
+static_assert(!std::is_copy_constructible_v<InterfaceElement>);
+static_assert(!std::is_copy_assignable_v<InterfaceElement>);
+static_assert(std::is_move_constructible_v<InterfaceElement>);
+static_assert(std::is_move_assignable_v<InterfaceElement>);
 } // namespace Kratos

…stom_elements/line_interface_element.cpp …on/custom_elements/interface_element.cppapplications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp renamed to applications/GeoMechanicsApplication/custom_elements/interface_element.cpp applications/GeoMechanicsApplication/custom_elements/line_interface_element.cpp renamed to applications/GeoMechanicsApplication/custom_elements/interface_element.cpp

@@ -8,6 +8,7 @@
 //  License:         geo_mechanics_application/license.txt
 //
 //  Main authors:    Richard Faasse
+//                   Gennady Markelov
 //
 
 #pragma once
@@ -21,28 +22,30 @@
 namespace Kratos
 {
 
-class KRATOS_API(GEO_MECHANICS_APPLICATION) LineInterfaceElement : public Element
+class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceElement : public Element
 {
 public:
-    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(LineInterfaceElement);
+    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(InterfaceElement);
 
     using Element::GeometryType;
     using Element::PropertiesType;
 
     // Following the UPwBaseElement example, we will follow the rule of 5
     // to avoid the noexcept code smell.
-    LineInterfaceElement()                                           = default;
-    ~LineInterfaceElement() override                                 = default;
-    LineInterfaceElement(const LineInterfaceElement&)                = delete;
-    LineInterfaceElement& operator=(const LineInterfaceElement&)     = delete;
-    LineInterfaceElement(LineInterfaceElement&&) noexcept            = default;
-    LineInterfaceElement& operator=(LineInterfaceElement&&) noexcept = default;
-
-    LineInterfaceElement(IndexType                      NewId,
-                         const GeometryType::Pointer&   rGeometry,
-                         const PropertiesType::Pointer& rProperties);
-
-    LineInterfaceElement(IndexType NewId, const GeometryType::Pointer& rGeometry);
+    ~InterfaceElement() override                             = default;
+    InterfaceElement(const InterfaceElement&)                = delete;
+    InterfaceElement& operator=(const InterfaceElement&)     = delete;
+    InterfaceElement(InterfaceElement&&) noexcept            = default;
+    InterfaceElement& operator=(InterfaceElement&&) noexcept = default;
+
+    InterfaceElement(IndexType                          NewId,
+                     const GeometryType::Pointer&       rGeometry,
+                     const PropertiesType::Pointer&     rProperties,
+                     std::unique_ptr<StressStatePolicy> pStressStatePolicy);
+
+    InterfaceElement(IndexType                          NewId,
+                     const GeometryType::Pointer&       rGeometry,
+                     std::unique_ptr<StressStatePolicy> pStressStatePolicy);
     Element::Pointer Create(IndexType NewId, const NodesArrayType& rNodes, PropertiesType::Pointer pProperties) const override;
     Element::Pointer Create(IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties) const override;
 
@@ -65,18 +68,24 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) LineInterfaceElement : public Elemen
     int  Check(const ProcessInfo& rCurrentProcessInfo) const override;
 
 private:
+    InterfaceElement() = default;
+
     Element::DofsVectorType GetDofs() const;
 
     std::vector<Matrix> CalculateLocalBMatricesAtIntegrationPoints() const;
     std::vector<double> CalculateIntegrationCoefficients() const;
     std::vector<Matrix> CalculateConstitutiveMatricesAtIntegrationPoints();
     std::vector<Vector> CalculateRelativeDisplacementsAtIntegrationPoints(const std::vector<Matrix>& rLocalBMatrices) const;
     std::vector<Vector> CalculateTractionsAtIntegrationPoints(const std::vector<Vector>& rRelativeDisplacements);
+    void MakeIntegrationSchemeAndAssignFunction();
+    std::function<Matrix(const Geometry<Node>&, const array_1d<double, 3>&)> mfpCalculateRotationMatrix;
 
     std::unique_ptr<IntegrationScheme>    mIntegrationScheme;
-    std::unique_ptr<StressStatePolicy>    mStressStatePolicy;
+    std::unique_ptr<StressStatePolicy>    mpStressStatePolicy;
     std::vector<ConstitutiveLaw::Pointer> mConstitutiveLaws;
     IntegrationCoefficientsCalculator     mIntegrationCoefficientsCalculator;
+
+    friend class Serializer;
 };
 
 } // namespace Kratos

…custom_elements/line_interface_element.h …tion/custom_elements/interface_element.happlications/GeoMechanicsApplication/custom_elements/line_interface_element.h renamed to applications/GeoMechanicsApplication/custom_elements/interface_element.h applications/GeoMechanicsApplication/custom_elements/line_interface_element.h renamed to applications/GeoMechanicsApplication/custom_elements/interface_element.h

@@ -35,4 +35,27 @@ Matrix GeometryUtilities::Calculate2DRotationMatrixForLineGeometry(const Geometr
     return result;
 }
 
+Matrix GeometryUtilities::Calculate3DRotationMatrixForPlaneGeometry(const Geometry<Node>& rGeometry,
+                                                                    const array_1d<double, 3>& rLocalCoordinate)
+{
+    Matrix jacobian;
+    rGeometry.Jacobian(jacobian, rLocalCoordinate);
+    const auto tangential_vector_1 = GeoMechanicsMathUtilities::Normalized(Vector{column(jacobian, 0)});
+    const auto tangential_vector_2 = GeoMechanicsMathUtilities::Normalized(
+        Vector{Vector{column(jacobian, 1)} -
+               inner_prod(Vector{column(jacobian, 1)}, tangential_vector_1) * tangential_vector_1});
+    Vector normal_vector(3);
+    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;
+}
+
 } // namespace Kratos

applications/GeoMechanicsApplication/custom_utilities/geometry_utilities.cpp

@@ -26,6 +26,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeometryUtilities
 public:
     static Matrix Calculate2DRotationMatrixForLineGeometry(const Geometry<Node>& rGeometry,
                                                            const array_1d<double, 3>& rLocalCoordinate);
+    static Matrix Calculate3DRotationMatrixForPlaneGeometry(const Geometry<Node>& rGeometry,
+                                                            const array_1d<double, 3>& rLocalCoordinate);
 };
 
 } // namespace Kratos