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 f97dea1a7338eb94ec0315c80e4e2f6daaf317fd..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"
@@ -75,7 +76,8 @@ public:
case 3:
return getWeightedPoint<MathLib::GaussLegendrePyramid<3>>(igp);
}
- return WeightedPoint(std::array<double, 3>(), 0);
+ OGS_FATAL("Integration order {:d} not implemented for pyramids.",
+ order);
}
template <typename Method>
@@ -101,7 +103,8 @@ public:
case 3:
return MathLib::GaussLegendrePyramid<3>::NPoints;
}
- return 0;
+ OGS_FATAL("Integration order {:d} not implemented for pyramids.",
+ order);
}
private:
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreRegular-impl.h b/NumLib/Fem/Integration/IntegrationGaussLegendreRegular-impl.h
index ab20b3848944a52730c3c42bb9864711613ce277..10d2643d5c263341150d92c67209f67115519922 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreRegular-impl.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreRegular-impl.h
@@ -70,9 +70,7 @@ IntegrationGaussLegendreRegular<N_DIM>::getWeightedPoint(unsigned order,
case 4:
return getWeightedPoint<MathLib::GaussLegendre<4>>(pos);
}
-
- return MathLib::TemplateWeightedPoint<double, double, N_DIM>(
- std::array<double, N_DIM>(), 0);
+ OGS_FATAL("Integration order {:d} not implemented.", order);
}
template <unsigned N_DIM>
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 ddec59ff00e3bab8f1ebf52e48cf79114950892f..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,8 +75,11 @@ public:
return getWeightedPoint<MathLib::GaussLegendreTet<2>>(igp);
case 3:
return getWeightedPoint<MathLib::GaussLegendreTet<3>>(igp);
+ case 4:
+ return getWeightedPoint<MathLib::GaussLegendreTet<4>>(igp);
}
- return WeightedPoint(std::array<double, 3>(), 0);
+ OGS_FATAL("Integration order {:d} not implemented for tetrahedrals.",
+ order);
}
template <typename Method>
@@ -100,8 +104,11 @@ public:
return MathLib::GaussLegendreTet<2>::NPoints;
case 3:
return MathLib::GaussLegendreTet<3>::NPoints;
+ case 4:
+ return MathLib::GaussLegendreTet<4>::NPoints;
}
- return 0;
+ OGS_FATAL("Integration order {:d} not implemented for tetrahedrals.",
+ order);
}
private:
diff --git a/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h b/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h
index 5202ac12d289ed68f459eb96ba37140f0b8fa868..62d13e15443408654f615a26e5d7ce31ec4082e3 100644
--- a/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h
+++ b/NumLib/Fem/Integration/IntegrationGaussLegendreTri.h
@@ -92,7 +92,8 @@ public:
case 4:
return getWeightedPoint<MathLib::GaussLegendreTri<4>>(igp);
}
- return WeightedPoint(std::array<double, 2>(), 0);
+ OGS_FATAL("Integration order {:d} not implemented for triangles.",
+ order);
}
template <typename Method>
@@ -120,7 +121,8 @@ public:
case 4:
return MathLib::GaussLegendreTri<4>::NPoints;
}
- return 0;
+ OGS_FATAL("Integration order {:d} not implemented for triangles.",
+ order);
}
private:
diff --git a/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test.prj b/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test.prj
index 9ea010b18c9ec45832c25f82e8291f05cbf39092..9711490ae07e265fc64e9d4fda3bd02a56c235e8 100644
--- a/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test.prj
+++ b/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test.prj
@@ -6,7 +6,7 @@
<process>
<name>LiquidFlow</name>
<type>LIQUID_FLOW</type>
- <integration_order>10</integration_order>
+ <integration_order>3</integration_order>
<process_variables>
<process_variable>pressure</process_variable>
</process_variables>
@@ -140,19 +140,19 @@
<geometrical_set>geometry</geometrical_set>
<geometry>tp3</geometry>
<type>Dirichlet</type>
- <parameter>p_bc_out</parameter>
+ <parameter>p_bc_in</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>geometry</geometrical_set>
<geometry>tp4</geometry>
<type>Dirichlet</type>
- <parameter>p_bc_out</parameter>
+ <parameter>p_bc_in</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>geometry</geometrical_set>
<geometry>tp5</geometry>
<type>Dirichlet</type>
- <parameter>p_bc_out</parameter>
+ <parameter>p_bc_in</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
diff --git a/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test_ts_1_t_1.000000.vtu b/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test_ts_1_t_1.000000.vtu
index 6fe1f50258a78567ecf66e3e49cf41b5bd8c3b75..9774c92c476513d260e52cee6b28f46f8eae9493 100644
--- a/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test_ts_1_t_1.000000.vtu
+++ b/Tests/Data/Utils/GMSH2OGS/quadratic_mesh_assembly_test_ts_1_t_1.000000.vtu
@@ -2,27 +2,27 @@
<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
<UnstructuredGrid>
<FieldData>
- <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="26" format="appended" RangeMin="45" RangeMax="121" offset="0" />
+ <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="21" format="appended" RangeMin="45" RangeMax="103" offset="0" />
</FieldData>
<Piece NumberOfPoints="37" NumberOfCells="8" >
<PointData>
- <DataArray type="Float64" Name="HydraulicFlow" format="appended" RangeMin="-0.01635" RangeMax="0.01635" offset="92" />
- <DataArray type="Float64" Name="pressure" format="appended" RangeMin="95095" RangeMax="2000000" offset="260" />
+ <DataArray type="Float64" Name="HydraulicFlow" format="appended" RangeMin="-3.7327" RangeMax="4.208175" offset="84" />
+ <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1000000" RangeMax="2000000" offset="472" />
</PointData>
<CellData>
- <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0" RangeMax="0" offset="444" />
+ <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0" RangeMax="0" offset="896" />
</CellData>
<Points>
- <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="2.4494897428" offset="504" />
+ <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="2.4494897428" offset="956" />
</Points>
<Cells>
- <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="708" />
- <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="936" />
- <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="1020" />
+ <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="1160" />
+ <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="1388" />
+ <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="1472" />
</Cells>
</Piece>
</UnstructuredGrid>
<AppendedData encoding="base64">
- _AQAAAAAAAAAAgAAAAAAAABoAAAAAAAAAIgAAAAAAAAA=eF4z0zPRM9Q1s9BNNzc0S0s0MDA2NNNLySwqqQQAUHsG4A==AQAAAAAAAAAAgAAAAAAAACgBAAAAAAAAXAAAAAAAAAA=eF5jYACC6EYbEMWgkWFzukfjLe++CfZgfuXSPWD65VSwvIPOYpvZUHnnXgitAaH3g+QbDhyCqGf4YM2ADBY37GGgEuBGsg8GGrw32eBSTwi868E0DxsAAIRYJfQ=AQAAAAAAAAAAgAAAAAAAACgBAAAAAAAAZwAAAAAAAAA=eF7jYmBgKDD/7sAJpbmgNC+QnjDnJ1gcRDNBxf/+//8/H0jD+DxI6hiyfjgwQ2kGKC0ElQdxG1rsHGHiyDRMvKFFzxGfOhD9B2h/P9Q8EP8jkP8/E1MdNhqbuSwMCPfh0w8A4ABO/A==AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAAHgDAAAAAAAAdwAAAAAAAAA=eF6NkjEOACEIBPmZ/Nwv3BMsLS0trzkuYXSj25BdNoCg2Q6j5Ei9CV1B5dmHsQmudIJ9lY91qat9mNsen678Sz3PtlO/f66aOeefYk/UT37mCXW32z4D/gD3s9zPM23wL+/wzNW/Cjyh18z7pR6cdXp5AVnjVtI=AQAAAAAAAAAAgAAAAAAAAKACAAAAAAAAiAAAAAAAAAA=eF6l0ckOgkAURNFWcQDFEVQQRdD//0YX3ruwE8LCtzkJhEd1dQjfmeA0/I7P59H7BS5xhgn6nbPCFN2T4Rr/3ePEuQ84lP+EJRZojg26z/05bnGHezzjFf2v+Rt8YIt3vGCFNd4w7sucTtyD/ZnriUO92Kt5PcfYPR1x7L467PGFb/wAkn0Dig==AQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAAHgAAAAAAAAA=eF5jZoAATigtBKWloLQKlLaA0u5QOgRKAwAbOAEwAQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAAEAAAAAAAAAA=eF4TFZMRl5CUkgYAA18AxQ==
+ _AQAAAAAAAAAAgAAAAAAAABUAAAAAAAAAHQAAAAAAAAA=eF4z0zPRM9Q1NDfQTU8xSElMMzFOMjMAADVzBT8=AQAAAAAAAAAAgAAAAAAAACgBAAAAAAAAAgEAAAAAAAA=eF77MLWm2Mr8hX1IcsTXsmQ+h3fLKiLtu57ab5LZtueQ/DP7g+ul9q6bwe1ww7YyYkUok8PM3dMm8Ff932/do/GWd98E+w9dMlenfV5v/+LTJkOJoy/2V5iYWM5/fcaeAQwCbD72BZeolDM5iCcevqx9VcBhIkibAJtDxXWhT47vWR1aLa4dzT0i4HDhDAi8sc9JAwE+h4xVCSFB6pf2Z0JpI2MQ+Lx/S+TXnbdOMR24bu13ceId3gMgWx6oHNgzi0N6XtxKpgMmRzbq5SkzHcj5/iRxYRjvgbNQdcwQBzn8+w8C7/enQOXZe8H+2L8IrP8fnLYEm/NvvwmUBgDdQI1SAQAAAAAAAAAAgAAAAAAAACgBAAAAAAAAHQEAAAAAAAA=eF6T+Orb4hBt5Djzi77KTU8jx2lXmztuxBg5Jgc9ruVpMXLMc9NZdS7DyDFx+jSdEGYjR0OzW9LdgkaOa8ptnv02NHJc/tGAjXuxkeM1uyWdysVGjrt3OnOKzTByVLAJO2xSbeS4t12qWo3ByJEBCBpa9BxZb9xh/ZVk6HhH/onirmBDuDgqbedYoz8zbNVWQ0dW2/lLwtYZOvLnHDl2OdnQcUW0sp3bc0PHf4fPnBdWM3S8zHW6dMIVQ8eT28/anK00dEwyv8B3IczQ8XppjfEvRqD4/Y7a++8N4OaXs3HP0coydHz+bf3MslMGjutivrNsjjBy7LgQ/HjfBRNHZ+YvVv/Pmzgef5qWu3OfiWP7zIQ2tgMmjgAIhHVFAQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAAHgDAAAAAAAAdwAAAAAAAAA=eF6NkjEOACEIBPmZ/Nwv3BMsLS0trzkuYXSj25BdNoCg2Q6j5Ei9CV1B5dmHsQmudIJ9lY91qat9mNsen678Sz3PtlO/f66aOeefYk/UT37mCXW32z4D/gD3s9zPM23wL+/wzNW/Cjyh18z7pR6cdXp5AVnjVtI=AQAAAAAAAAAAgAAAAAAAAKACAAAAAAAAiAAAAAAAAAA=eF6l0ckOgkAURNFWcQDFEVQQRdD//0YX3ruwE8LCtzkJhEd1dQjfmeA0/I7P59H7BS5xhgn6nbPCFN2T4Rr/3ePEuQ84lP+EJRZojg26z/05bnGHezzjFf2v+Rt8YIt3vGCFNd4w7sucTtyD/ZnriUO92Kt5PcfYPR1x7L467PGFb/wAkn0Dig==AQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAAHgAAAAAAAAA=eF5jZoAATigtBKWloLQKlLaA0u5QOgRKAwAbOAEwAQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAAEAAAAAAAAAA=eF4TFZMRl5CUkgYAA18AxQ==
</AppendedData>
</VTKFile>
diff --git a/Tests/MathLib/TestGaussLegendreIntegration.cpp b/Tests/MathLib/TestGaussLegendreIntegration.cpp
index 5e6ff1b0e8c9a8c1e2072a66a27846ebb6539580..c6b4f49e2653a83af8d715c248856abdf3e49a51 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,15 +564,20 @@ 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})
+ // quadrature must be exact up to the given polynomial degrees
+ // for integration order n = 1, 2, 3 it is 2*n-1 = 1, 3, 5
+ // for integration order 4 it is 5, i.e. the same as for 3rd integration
+ // order
+ const std::vector<unsigned> maximum_polynomial_degrees{0, 1, 3, 5, 5};
+
+ for (unsigned integration_order : {1, 2, 3, 4})
{
DBUG("\n==== integration order: {:d}.\n", integration_order);
GaussLegendreTest::IntegrationTestProcess pcs_tet(*mesh_tet,
integration_order);
for (unsigned polynomial_order = 0;
- // Gauss-Legendre integration is exact up to this order!
- polynomial_order < 2 * integration_order;
+ polynomial_order <= maximum_polynomial_degrees[integration_order];
++polynomial_order)
{
DBUG(" == polynomial order: {:d}.", polynomial_order);