From c61eb5cdec4ee89e014a2da04cb8dc590e640c70 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B6rg=20Buchwald?= <joerg.buchwald@ufz.de>
Date: Thu, 4 Nov 2021 22:25:21 +0100
Subject: [PATCH] add integration order 4 for Tets
---
MathLib/Integration/GaussLegendreTet.cpp | 42 +++++++++++++++++++
MathLib/Integration/GaussLegendreTet.h | 9 ++++
.../IntegrationGaussLegendrePyramid.h | 1 +
.../IntegrationGaussLegendreRegular.h | 1 +
.../Integration/IntegrationGaussLegendreTet.h | 5 +++
.../MathLib/TestGaussLegendreIntegration.cpp | 6 +--
6 files changed, 61 insertions(+), 3 deletions(-)
diff --git a/MathLib/Integration/GaussLegendreTet.cpp b/MathLib/Integration/GaussLegendreTet.cpp
index 593b251cb16..afa8e7e4672 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 9dd8370ed09..5897905ec63 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 2a679dff4d7..088967d5b3b 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 a9b70ba7bd1..b7c4d9cbc5f 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 ebbbcd887c8..1d911f427cf 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 5e6ff1b0e8c..d648130f5c6 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,
--
GitLab