From cac75d7f2b200244f6cbda718b06c0d063fd020d Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 15 Aug 2017 15:00:23 +0200
Subject: [PATCH] [MaL] fixed Gauss quadrature for Tet 2rd order
---
MathLib/Integration/GaussLegendreTet.cpp | 66 +++++++++++++-----------
MathLib/Integration/GaussLegendreTet.h | 2 +-
2 files changed, 36 insertions(+), 32 deletions(-)
diff --git a/MathLib/Integration/GaussLegendreTet.cpp b/MathLib/Integration/GaussLegendreTet.cpp
index b81cc7272b3..1faa686b15a 100644
--- a/MathLib/Integration/GaussLegendreTet.cpp
+++ b/MathLib/Integration/GaussLegendreTet.cpp
@@ -25,37 +25,41 @@ GaussLegendreTet<2>::X = {{ {{1./4., 1./4., 1./4.}},
{{1./6., 1./6., 1./2.}} }};
double const GaussLegendreTet<2>::W[5] = {-2./15., 0.075, 0.075, 0.075, 0.075};
+std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints> initGLTet3X()
+{
+ // 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.0673422422100983;
+ const double b = 0.3108859192633005;
+ const double c = 0.7217942490673264;
+ const double d = 0.0927352503108912;
+ const double e = 0.4544962958743506;
+ const double f = 0.045503704125649;
+
+ return {{{a, b, b},
+ {b, a, b},
+ {b, b, a},
+ {b, b, b},
+ {c, d, d},
+ {d, c, d},
+ {d, d, c},
+ {d, d, d},
+ {e, e, f},
+ {e, f, e},
+ {e, f, f},
+ {f, e, e},
+ {f, e, f},
+ {f, f, e}}};
+}
+
const std::array<std::array<double, 3>, GaussLegendreTet<3>::NPoints>
-GaussLegendreTet<3>::X = {{ {{1./2., 1./2., 1./2.}},
- {{0.09197107805272303, 0.09197107805272303, 0.09197107805272303}},
- {{0.72408676584183096, 0.09197107805272303, 0.09197107805272303}},
- {{0.09197107805272303, 0.72408676584183096, 0.09197107805272303}},
- {{0.09197107805272303, 0.09197107805272303, 0.72408676584183096}},
- {{0.44364916731037080, 0.05635083268962915, 0.05635083268962915}},
- {{0.05635083268962915, 0.44364916731037080, 0.05635083268962915}},
- {{0.05635083268962915, 0.05635083268962915, 0.44364916731037080}},
- {{0.05635083268962915, 0.44364916731037080, 0.44364916731037080}},
- {{0.44364916731037080, 0.05635083268962915, 0.44364916731037080}},
- {{0.44364916731037080, 0.44364916731037080, 0.05635083268962915}},
- {{0.31979362782962989, 0.31979362782962989, 0.31979362782962989}},
- {{0.04061911651111023, 0.31979362782962989, 0.31979362782962989}},
- {{0.31979362782962989, 0.04061911651111023, 0.31979362782962989}},
- {{0.31979362782962989, 0.31979362782962989, 0.04061911651111023}} }};
-double const GaussLegendreTet<3>::W[15] = {0.019753086419753086,
- 0.011989513963169772,
- 0.011989513963169772,
- 0.011989513963169772,
- 0.011989513963169772,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.008818342151675485,
- 0.011511367871045397,
- 0.011511367871045397,
- 0.011511367871045397,
- 0.011511367871045397
- };
+ GaussLegendreTet<3>::X = initGLTet3X();
+
+static const double p = 0.1126879257180162 / 6.;
+static const double q = 0.0734930431163619 / 6.;
+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};
}
diff --git a/MathLib/Integration/GaussLegendreTet.h b/MathLib/Integration/GaussLegendreTet.h
index 4a6ac7e6c03..110635849f5 100644
--- a/MathLib/Integration/GaussLegendreTet.h
+++ b/MathLib/Integration/GaussLegendreTet.h
@@ -38,7 +38,7 @@ struct GaussLegendreTet<2> {
template <>
struct GaussLegendreTet<3> {
static MATHLIB_EXPORT const unsigned Order = 3;
- static MATHLIB_EXPORT const unsigned NPoints = 15;
+ static MATHLIB_EXPORT const unsigned NPoints = 14;
static MATHLIB_EXPORT const std::array<std::array<double, 3>, NPoints> X;
static MATHLIB_EXPORT const double W[NPoints];
};
--
GitLab