diff --git a/MathLib/Integration/GaussLegendreTet.cpp b/MathLib/Integration/GaussLegendreTet.cpp
index 593b251cb16cbe6ccd800e817994150643729c3b..afa8e7e4672af5928ffecb032b0eb4d7f94cc1dd 100644
--- a/MathLib/Integration/GaussLegendreTet.cpp
+++ b/MathLib/Integration/GaussLegendreTet.cpp
@@ -65,4 +65,46 @@ static const double r = 0.0425460207770812 / 6.;
double const GaussLegendreTet<3>::W[GaussLegendreTet<3>::NPoints] = {
p, p, p, p, q, q, q, q, r, r, r, r, r, r};
+
+static std::array<std::array<double, 3>, GaussLegendreTet<4>::NPoints>
+initGLTet4X()
+{
+ // Cf. Gellert, M., Harbord, R., 1991. Moderate degree cubature formulas for
+ // 3-D tetrahedral finite-element approximations. Communications in Applied
+ // Numerical Methods 7, 487-495. doi:10.1002/cnm.1630070609
+ const double a = 0.3797582452067875;
+ const double b = 0.1202417547932126;
+
+ return {{{{0.0, 1. / 3, 1. / 3}},
+ {{1. / 3, 0.0, 1. / 3}},
+ {{1. / 3, 1. / 3, 0.0}},
+ {{1. / 3, 1. / 3, 1. / 3}},
+ {{8. / 11, 1. / 11, 1. / 11}},
+ {{1. / 11, 8. / 11, 1. / 11}},
+ {{1. / 11, 1. / 11, 8. / 11}},
+ {{1. / 11, 1. / 11, 1. / 11}},
+ {{0.0, 0.0, 0.5}},
+ {{0.0, 0.5, 0.0}},
+ {{0.0, 0.5, 0.5}},
+ {{0.5, 0.0, 0.0}},
+ {{0.5, 0.0, 0.5}},
+ {{0.5, 0.5, 0.0}},
+ {{a, a, b}},
+ {{a, b, a}},
+ {{a, b, b}},
+ {{b, a, a}},
+ {{b, a, b}},
+ {{b, b, a}}}};
+}
+const std::array<std::array<double, 3>, GaussLegendreTet<4>::NPoints>
+ GaussLegendreTet<4>::X = initGLTet4X();
+
+static const double s = 81. / 2240. / 6.;
+static const double t = 161051. / 2304960. / 6.;
+static const double u = 409. / 31395. / 6.;
+static const double v = 2679769. / 32305455. / 6.;
+
+double const GaussLegendreTet<4>::W[GaussLegendreTet<4>::NPoints] = {
+ s, s, s, s, t, t, t, t, u, u, u, u, u, u, v, v, v, v, v, v};
+
} // namespace MathLib
diff --git a/MathLib/Integration/GaussLegendreTet.h b/MathLib/Integration/GaussLegendreTet.h
index 9dd8370ed0941b3daa8f9b6034c5140041200ee1..5897905ec631b9902217326bb5b2266176c08a67 100644
--- a/MathLib/Integration/GaussLegendreTet.h
+++ b/MathLib/Integration/GaussLegendreTet.h
@@ -44,6 +44,15 @@ struct GaussLegendreTet<3> {
static MATHLIB_EXPORT const double W[NPoints];
};
+template <>
+struct GaussLegendreTet<4>
+{
+ static MATHLIB_EXPORT const unsigned Order = 4;
+ static MATHLIB_EXPORT const unsigned NPoints = 20;
+ static MATHLIB_EXPORT const std::array<std::array<double, 3>, NPoints> X;
+ static MATHLIB_EXPORT const double W[NPoints];
+};
+
#ifndef _MSC_VER // The following explicit instantatiation declaration does not
// compile on that particular compiler but is necessary.
template <>
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h b/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h
index 2a679dff4d7050a3d1e6cfd5bbd8ba733d9648c5..088967d5b3bf3e3730f0ccd949e7676b34851c7f 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendrePyramid.h
@@ -10,6 +10,7 @@
#pragma once
+#include "BaseLib/Error.h"
#include "MathLib/Integration/GaussLegendrePyramid.h"
#include "MathLib/TemplateWeightedPoint.h"
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h b/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h
index a9b70ba7bd1a687c4aa284560f3f689d733f632d..b7c4d9cbc5ff4bbb5a13353eeefbc8f5c0d8af74 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreRegular.h
@@ -15,6 +15,7 @@
#include <array>
#include <cmath>
+#include "BaseLib/Error.h"
#include "MathLib/Integration/GaussLegendre.h"
#include "MathLib/TemplateWeightedPoint.h"
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h b/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h
index ebbbcd887c8a06e0cb91f85d4df00ddf833c3dfa..1d911f427cf18ed50e9a929e76fa64aa04972014 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreTet.h
@@ -10,6 +10,7 @@
#pragma once
+#include "BaseLib/Error.h"
#include "MathLib/Integration/GaussLegendreTet.h"
#include "MathLib/TemplateWeightedPoint.h"
@@ -74,6 +75,8 @@ public:
return getWeightedPoint<MathLib::GaussLegendreTet<2>>(igp);
case 3:
return getWeightedPoint<MathLib::GaussLegendreTet<3>>(igp);
+ case 4:
+ return getWeightedPoint<MathLib::GaussLegendreTet<4>>(igp);
}
OGS_FATAL("Integration order {:d} not implemented for tetrahedrals.",
order);
@@ -101,6 +104,8 @@ public:
return MathLib::GaussLegendreTet<2>::NPoints;
case 3:
return MathLib::GaussLegendreTet<3>::NPoints;
+ case 4:
+ return MathLib::GaussLegendreTet<4>::NPoints;
}
OGS_FATAL("Integration order {:d} not implemented for tetrahedrals.",
order);
diff --git a/Tests/MathLib/TestGaussLegendreIntegration.cpp b/Tests/MathLib/TestGaussLegendreIntegration.cpp
index 5e6ff1b0e8c9a8c1e2072a66a27846ebb6539580..d648130f5c60972dc6cb70439fb2a17fb6ae730c 100644
--- a/Tests/MathLib/TestGaussLegendreIntegration.cpp
+++ b/Tests/MathLib/TestGaussLegendreIntegration.cpp
@@ -447,7 +447,7 @@ OGS_DONT_TEST_THIS_IF_PETSC(MathLib, IntegrationGaussLegendreTet)
MeshLib::IO::VtuInterface::readVTUFile(
TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_tet.vtu"));
- for (unsigned integration_order : {1, 2, 3})
+ for (unsigned integration_order : {1, 2, 3, 4})
{
DBUG("\n==== integration order: {:d}.\n", integration_order);
GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
@@ -508,7 +508,7 @@ OGS_DONT_TEST_THIS_IF_PETSC(
MeshLib::IO::VtuInterface::readVTUFile(
TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_tet.vtu"));
- for (unsigned integration_order : {1, 2, 3})
+ for (unsigned integration_order : {1, 2, 3, 4})
{
DBUG("\n==== integration order: {:d}.\n", integration_order);
GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
@@ -564,7 +564,7 @@ OGS_DONT_TEST_THIS_IF_PETSC(MathLib,
MeshLib::IO::VtuInterface::readVTUFile(
TestInfoLib::TestInfo::data_path + "/MathLib/unit_cube_tet.vtu"));
- for (unsigned integration_order : {1, 2, 3})
+ for (unsigned integration_order : {1, 2, 3, 4})
{
DBUG("\n==== integration order: {:d}.\n", integration_order);
GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,