[GeoMechanicsApplication] Lumped integration for interface elements (#13590)

Lumped integration scheme for triangular and quadrilateral elements, with unit tests.

Author: WPK4FEM

applications/GeoMechanicsApplication/custom_elements/lobatto_integration_scheme.cpp

@@ -36,10 +36,10 @@ Geo::IntegrationPointVectorType LobattoIntegrationScheme::CreateIntegrationPoint
     // instance, here: https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Lobatto_rules
     switch (NumberOfPoints) {
     case 2:
-        return {{-1.0, 1.0}, {1.0, 1.0}};
+        return {{Point(-1.0), 1.0}, {Point(1.0), 1.0}};
 
     case 3:
-        return {{-1.0, 1.0 / 3.0}, {0.0, 4.0 / 3.0}, {1.0, 1.0 / 3.0}};
+        return {{Point(-1.0), 1.0 / 3.0}, {Point(0.0), 4.0 / 3.0}, {Point(1.0), 1.0 / 3.0}};
 
     default:
         KRATOS_ERROR << "Can't construct Lobatto integration scheme: no support for "

applications/GeoMechanicsApplication/custom_elements/lumped_integration_scheme.cpp

@@ -0,0 +1,101 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Wijtze Pieter Kikstra
+//
+
+#include "lumped_integration_scheme.h"
+
+#include "geometries/geometry.h"
+#include "geometries/quadrilateral_2d_4.h"
+#include "geometries/quadrilateral_2d_8.h"
+#include "geometries/triangle_2d_3.h"
+#include "geometries/triangle_2d_6.h"
+#include "includes/node.h"
+
+namespace
+{
+
+using namespace Kratos;
+
+template <typename GeometryType>
+Geo::IntegrationPointVectorType MakeLumpedIntegrationPointsVector(const std::vector<Point>& rPoints)
+{
+    std::vector<std::size_t> node_ids(rPoints.size());
+    std::iota(node_ids.begin(), node_ids.end(), std::size_t{1});
+    auto node_ptrs = Geometry<Node>::PointsArrayType{};
+    for (auto i = std::size_t{0}; i < node_ids.size(); ++i) {
+        node_ptrs.push_back(make_intrusive<Node>(node_ids[i], rPoints[i]));
+    }
+
+    auto geom            = GeometryType{node_ptrs};
+    auto lumping_factors = Vector{};
+    geom.LumpingFactors(lumping_factors, Geometry<Node>::LumpingMethods::DIAGONAL_SCALING);
+
+    auto make_integration_point = [](const Point& rPoint, double Weight) {
+        return Geo::IntegrationPointType{rPoint, Weight};
+    };
+    Geo::IntegrationPointVectorType result;
+    std::transform(rPoints.begin(), rPoints.end(), lumping_factors.cbegin(),
+                   std::back_inserter(result), make_integration_point);
+    return result;
+}
+
+} // namespace
+
+namespace Kratos
+{
+
+LumpedIntegrationScheme::LumpedIntegrationScheme(std::size_t NumberOfPoints)
+    : mIntegrationPoints(CreateIntegrationPoints(NumberOfPoints))
+{
+}
+
+std::size_t LumpedIntegrationScheme::GetNumberOfIntegrationPoints() const
+{
+    return mIntegrationPoints.size();
+}
+
+const Geo::IntegrationPointVectorType& LumpedIntegrationScheme::GetIntegrationPoints() const
+{
+    return mIntegrationPoints;
+}
+
+Geo::IntegrationPointVectorType LumpedIntegrationScheme::CreateIntegrationPoints(std::size_t NumberOfPoints)
+{
+    // Integration scheme used for planar interface elements in GeoMechanicsApplication
+    // locations coincide with node-pairs, numbering follows the node-pairs too, weights as in lumped mass matrix from diagonal scaling of a surface element ( see geometry.h )
+    switch (NumberOfPoints) {
+    case 3: {
+        const auto points = std::vector{Point(0.0, 0.0), Point(1.0, 0.0), Point(0.0, 1.0)};
+        return MakeLumpedIntegrationPointsVector<Triangle2D3<Node>>(points);
+    }
+    case 4: {
+        const auto points =
+            std::vector{Point(-1.0, -1.0), Point(1.0, -1.0), Point(1.0, 1.0), Point(-1.0, 1.0)};
+        return MakeLumpedIntegrationPointsVector<Quadrilateral2D4<Node>>(points);
+    }
+    case 6: {
+        const auto points = std::vector{Point(0.0, 0.0), Point(1.0, 0.0), Point(0.0, 1.0),
+                                        Point(0.5, 0.0), Point(0.5, 0.5), Point(0.0, 0.5)};
+        return MakeLumpedIntegrationPointsVector<Triangle2D6<Node>>(points);
+    }
+    case 8: {
+        const auto points =
+            std::vector{Point(-1.0, -1.0), Point(1.0, -1.0), Point(1.0, 1.0), Point(-1.0, 1.0),
+                        Point(0.0, -1.0),  Point(1.0, 0.0),  Point(0.0, 1.0), Point(-1.0, 0.0)};
+        return MakeLumpedIntegrationPointsVector<Quadrilateral2D8<Node>>(points);
+    }
+    default:
+        KRATOS_ERROR << "Can't construct Lumped integration scheme: no support for "
+                     << NumberOfPoints << " point(s)" << std::endl;
+    }
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_elements/lumped_integration_scheme.h

@@ -0,0 +1,41 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Wijtze Pieter Kikstra
+//
+
+#pragma once
+
+#include "includes/kratos_export_api.h"
+#include "integration_scheme.h"
+
+namespace Kratos
+{
+/**
+ * @Class LumpedIntegrationScheme
+ * @brief Integration scheme used for planar interface elements in GeoMechanicsApplication
+ * @details Integration point locations coincide with node-pairs, numbering follows the node-pairs too. Integrations weights as in lumped mass matrix from diagonal scaling of a surface element ( see geometry.h )
+ * @author Wijtze Pieter Kikstra
+ */
+class KRATOS_API(GEO_MECHANICS_APPLICATION) LumpedIntegrationScheme : public IntegrationScheme
+{
+public:
+    explicit LumpedIntegrationScheme(std::size_t NumberOfPoints);
+    ~LumpedIntegrationScheme() override = default;
+
+    [[nodiscard]] std::size_t GetNumberOfIntegrationPoints() const override;
+    [[nodiscard]] const Geo::IntegrationPointVectorType& GetIntegrationPoints() const override;
+
+private:
+    static Geo::IntegrationPointVectorType CreateIntegrationPoints(std::size_t NumberOfPoints);
+
+    Geo::IntegrationPointVectorType mIntegrationPoints;
+};
+
+} // namespace Kratos