Merge pull request #12917 from KratosMultiphysics/correcting-hexa-qua…

…drature-lobatto [CORE] Correcting hexahedron quadrature Lobatto
KratosMultiphysics · Dec 9, 2024 · 888e7b3 · 888e7b3
2 parents 88697ed + 66df33a
commit 888e7b3
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 12 deletions.
diff --git a/kratos/integration/hexahedron_gauss_lobatto_integration_points.h b/kratos/integration/hexahedron_gauss_lobatto_integration_points.h
@@ -50,10 +50,10 @@ class HexahedronGaussLobattoIntegrationPoints1
     static const IntegrationPointsArrayType& IntegrationPoints()
     {
         static const IntegrationPointsArrayType s_integration_points{{
-            IntegrationPointType( -1.00 , -1.00, 0.00, 1.00 ),
-            IntegrationPointType(  1.00 , -1.00, 0.00, 1.00 ),
-            IntegrationPointType(  1.00 ,  1.00, 0.00, 1.00 ),
-            IntegrationPointType( -1.00 ,  1.00, 0.00, 1.00 )
+            IntegrationPointType( -1.0 , -1.0, 0.0, 1.0 ),
+            IntegrationPointType(  1.0 , -1.0, 0.0, 1.0 ),
+            IntegrationPointType(  1.0 ,  1.0, 0.0, 1.0 ),
+            IntegrationPointType( -1.0 ,  1.0, 0.0, 1.0 )
         }};
         return s_integration_points;
     }
@@ -90,14 +90,14 @@ class HexahedronGaussLobattoIntegrationPoints2
     static const IntegrationPointsArrayType& IntegrationPoints()
     {
         static const IntegrationPointsArrayType s_integration_points{{
-            IntegrationPointType( -1.00 , -1.00, -1.00, 0.50 ),
-            IntegrationPointType(  1.00 , -1.00, -1.00, 0.50 ),
-            IntegrationPointType(  1.00 ,  1.00, -1.00, 0.50 ),
-            IntegrationPointType( -1.00 ,  1.00, -1.00, 0.50 ),
-            IntegrationPointType( -1.00 , -1.00,  1.00, 0.50 ),
-            IntegrationPointType(  1.00 , -1.00,  1.00, 0.50 ),
-            IntegrationPointType(  1.00 ,  1.00,  1.00, 0.50 ),
-            IntegrationPointType( -1.00 ,  1.00,  1.00, 0.50 )
+            IntegrationPointType( -1.0 , -1.0, -1.0, 1.0 ),
+            IntegrationPointType(  1.0 , -1.0, -1.0, 1.0 ),
+            IntegrationPointType(  1.0 ,  1.0, -1.0, 1.0 ),
+            IntegrationPointType( -1.0 ,  1.0, -1.0, 1.0 ),
+            IntegrationPointType( -1.0 , -1.0,  1.0, 1.0 ),
+            IntegrationPointType(  1.0 , -1.0,  1.0, 1.0 ),
+            IntegrationPointType(  1.0 ,  1.0,  1.0, 1.0 ),
+            IntegrationPointType( -1.0 ,  1.0,  1.0, 1.0 )
         }};
         return s_integration_points;
     }

diff --git a/kratos/tests/cpp_tests/integration/test_integration_quadratures.cpp b/kratos/tests/cpp_tests/integration/test_integration_quadratures.cpp
@@ -17,6 +17,7 @@
 // Project includes
 #include "testing/testing.h"
 #include "integration/quadrilateral_gauss_lobatto_integration_points.h"
+#include "integration/hexahedron_gauss_lobatto_integration_points.h"
 
 namespace Kratos::Testing
 {
@@ -75,5 +76,31 @@ KRATOS_TEST_CASE_IN_SUITE(GaussLobattoQuadrilateralQuadraturesTest, KratosCoreFa
 
 }
 
+KRATOS_TEST_CASE_IN_SUITE(GaussLobattoHexaQuadraturesTest, KratosCoreFastSuite)
+{
+
+    // In This test we evaluate the Gauss-Lobatto quadratures for integrating
+    // f = (x+y+z+5) over a [-1, 1] hexa
+
+    const auto& r_lobatto = HexahedronGaussLobattoIntegrationPoints2();
+
+    // Analytical results, reference
+    const double integral_f = 40.0;
+
+    double quadrature_integral_f = 0.0;
+
+    // Integral for f with Lobatto 1
+    for (IndexType IP = 0; IP < r_lobatto.IntegrationPoints().size(); ++IP) {
+        const auto& r_IP = r_lobatto.IntegrationPoints()[IP];
+        const double X = r_IP.X();
+        const double Y = r_IP.Y();
+        const double Z = r_IP.Z();
+        const double w = r_IP.Weight();
+
+        quadrature_integral_f += w * (X + Y + Z + 5.0);
+    }
+
+    KRATOS_CHECK_NEAR(quadrature_integral_f, integral_f, 1.0E-6);
+}
 
 } // namespace Kratos::Testing