[CORE] Adding Lobatto quadrature for triangle

kratos/integration/triangle_gauss_legendre_integration_points.h

@@ -11,8 +11,7 @@
 //
 //
 
-#if !defined(KRATOS_TRIANGLE_GAUSS_LEGENDRE_INTEGRATION_POINTS_H_INCLUDED )
-#define  KRATOS_TRIANGLE_GAUSS_LEGENDRE_INTEGRATION_POINTS_H_INCLUDED
+#pragma once
 
 // System includes
 
@@ -87,7 +86,7 @@ class TriangleGaussLegendreIntegrationPoints2
         // Define auxiliary constants
         const double one_sixth = 1.0 / 6.0;
         const double two_thirds = 2.0 / 3.0;
-        const double weight = 1.0 / 6.0;
+        const double weight = one_sixth;
 
         // Integration points using the auxiliary constants
         static const IntegrationPointsArrayType s_integration_points{{
@@ -285,6 +284,3 @@ class TriangleGaussLegendreIntegrationPoints5
 
 }  // namespace Kratos.
 
-#endif // KRATOS_TRIANGLE_GAUSS_LEGENDRE_INTEGRATION_POINTS_H_INCLUDED  defined
-
-

kratos/integration/triangle_gauss_lobatto_integration_points.h

@@ -0,0 +1,78 @@
+//    |  /           |
+//    ' /   __| _` | __|  _ \   __|
+//    . \  |   (   | |   (   |\__ `
+//   _|\_\_|  \__,_|\__|\___/ ____/
+//                   Multi-Physics
+//
+//  License:         BSD License
+//                   Kratos default license: kratos/license.txt
+//
+//  Main authors:    Alejandro Cornejo
+//
+//
+
+#pragma once
+
+// System includes
+
+// External includes
+
+// Project includes
+#include "integration/quadrature.h"
+
+namespace Kratos
+{
+
+class TriangleGaussLobattoIntegrationPoints1
+{
+public:
+    KRATOS_CLASS_POINTER_DEFINITION(TriangleGaussLobattoIntegrationPoints1);
+    typedef std::size_t SizeType;
+
+    static const unsigned int Dimension = 2;
+
+    typedef IntegrationPoint<2> IntegrationPointType;
+
+    typedef std::array<IntegrationPointType, 3> IntegrationPointsArrayType;
+
+    typedef IntegrationPointType::PointType PointType;
+
+    static SizeType IntegrationPointsNumber()
+    {
+        return 3;
+    }
+
+    static const IntegrationPointsArrayType& IntegrationPoints()
+    {
+        const double one_over_six = 1.0 / 6.0;
+        static const IntegrationPointsArrayType s_integration_points{{
+            IntegrationPointType( 0.0, 0.0 , one_over_six ),
+            IntegrationPointType( 1.0, 0.0 , one_over_six ),
+            IntegrationPointType( 0.0, 1.0 , one_over_six )
+        }};
+        return s_integration_points;
+    }
+
+    std::string Info() const
+    {
+        std::stringstream buffer;
+        buffer << "Triangle Gauss-Lobatto quadrature 1 (3 points, degree of precision = 1)";
+        return buffer.str();
+    }
+
+
+}; // Class TriangleGaussLobattoIntegrationPoints1
+
+///@name Type Definitions
+///@{
+
+
+///@}
+///@name Input and output
+///@{
+
+
+///@}
+
+
+}  // namespace Kratos.

kratos/tests/cpp_tests/integration/test_integration_quadratures.cpp

@@ -18,6 +18,7 @@
 #include "testing/testing.h"
 #include "integration/quadrilateral_gauss_lobatto_integration_points.h"
 #include "integration/hexahedron_gauss_lobatto_integration_points.h"
+#include "integration/triangle_gauss_lobatto_integration_points.h"
 
 namespace Kratos::Testing
 {
@@ -103,4 +104,30 @@ KRATOS_TEST_CASE_IN_SUITE(GaussLobattoHexaQuadraturesTest, KratosCoreFastSuite)
     KRATOS_CHECK_NEAR(quadrature_integral_f, integral_f, 1.0E-6);
 }
 
+KRATOS_TEST_CASE_IN_SUITE(GaussLobattoTriangleQuadraturesTest, KratosCoreFastSuite)
+{
+
+    // In This test we evaluate the Gauss-Lobatto quadrature for integrating
+    // f = (x+y) over a [0, 1]x[0, 1-x] triangle
+
+    const auto& r_lobatto_1 = TriangleGaussLobattoIntegrationPoints1();
+
+    // Analytical result, reference
+    const double integral_f = 1.0 / 3.0;
+
+    double quadrature_integral_f = 0.0;
+
+    // Integral for f with Lobatto 1
+    for (IndexType IP = 0; IP < r_lobatto_1.IntegrationPoints().size(); ++IP) {
+        const auto& r_IP = r_lobatto_1.IntegrationPoints()[IP];
+        const double X = r_IP.X();
+        const double Y = r_IP.Y();
+        const double w = r_IP.Weight();
+
+        quadrature_integral_f += w * (X + Y);
+    }
+
+    KRATOS_CHECK_NEAR(quadrature_integral_f, integral_f, 1.0E-6);
+}
+
 } // namespace Kratos::Testing