diff --git a/NumLib/Fem/ShapeFunction/ShapeHex20-impl.h b/NumLib/Fem/ShapeFunction/ShapeHex20-impl.h
index 2ad6dd2b597746733fbdcae4732b9501b066c34f..ca24efa2d0eebec3f0e926ea776f940fd2ea4e23 100644
--- a/NumLib/Fem/ShapeFunction/ShapeHex20-impl.h
+++ b/NumLib/Fem/ShapeFunction/ShapeHex20-impl.h
@@ -12,12 +12,12 @@ namespace
inline double ShapeFunctionHexHQ_Corner(const double r, const double s, const double t)
{
- return 0.125 * (1 + r) * (1 + s) * (1 + t) * ( r + s + t - 2.0);
+ return 0.125 * (r - 1) * (s - 1) * (t - 1) * (r + s + t + 2.0);
}
inline double ShapeFunctionHexHQ_Middle(const double r, const double s, const double t)
{
- return 0.25 * (1 - r * r) * (1 + s) * (1 + t);
+ return 0.25 * (1 - r * r) * (s - 1) * (t - 1);
}
inline double dShapeFunctionHexHQ_Corner(const double r, const double s, const double t, const int ty)
@@ -25,11 +25,11 @@ inline double dShapeFunctionHexHQ_Corner(const double r, const double s, const d
switch(ty)
{
case 0:
- return 0.125 * (1 + s) * (1 + t) * (2.0 * r + s + t - 1.0);
+ return 0.125 * (s - 1) * (t - 1) * (2.0 * r + s + t + 1.0);
case 1:
- return 0.125 * (1 + t) * (1 + r) * (2.0 * s + r + t - 1.0);
+ return 0.125 * (t - 1) * (r - 1) * (2.0 * s + r + t + 1.0);
case 2:
- return 0.125 * (1 + r) * (1 + s) * (2.0 * t + s + r - 1.0);
+ return 0.125 * (r - 1) * (s - 1) * (2.0 * t + s + r + 1.0);
}
return 0.0;
}
@@ -39,11 +39,11 @@ inline double dShapeFunctionHexHQ_Middle(const double r, const double s, const d
switch(ty)
{
case 0:
- return -0.5 * r * (1 + s) * (1 + t);
+ return -0.5 * r * (s - 1) * (t - 1);
case 1:
- return 0.25 * (1 - r * r) * (1 + t);
+ return 0.25 * (1 - r * r) * (t - 1);
case 2:
- return 0.25 * (1 - r * r) * (1 + s);
+ return 0.25 * (1 - r * r) * (s - 1);
}
return 0.0;
}
diff --git a/ProcessLib/GroundwaterFlow/CMakeLists.txt b/ProcessLib/GroundwaterFlow/CMakeLists.txt
index 5c95ffbed06489b4f7b72a7487fd2e70a53e02f6..8e122a09ea178043446ef94fbed359721f7d4906 100644
--- a/ProcessLib/GroundwaterFlow/CMakeLists.txt
+++ b/ProcessLib/GroundwaterFlow/CMakeLists.txt
@@ -63,6 +63,19 @@ foreach(mesh_size 1e4 2e4 3e4 4e4 5e4 1e5 1e6)
)
endforeach()
+# Quadratic hex element.
+AddTest(
+ NAME GroundWaterFlowProcess_cube_1x1x1_1e0_QuadraticHex
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e0_quadratic_hex.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1e-15 RELTOL 1e-15
+ DIFF_DATA
+ cube_1x1x1_hex20_1e0.vtu cube_1e0_quadratic_hex_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
+)
+
# SQUARE 1x1 GROUNDWATER FLOW TESTS
foreach(mesh_size 1e0 1e1 1e2 1e3 1e4)
AddTest(
diff --git a/Tests/MathLib/AutoCheckTools.h b/Tests/AutoCheckTools.h
similarity index 89%
rename from Tests/MathLib/AutoCheckTools.h
rename to Tests/AutoCheckTools.h
index eed06985275c2f908c7f9604947dc7a709109dc3..dcd29d2751077e69e75ec644f9c9ad1505bf82df 100644
--- a/Tests/MathLib/AutoCheckTools.h
+++ b/Tests/AutoCheckTools.h
@@ -13,6 +13,22 @@
#include "MathLib/Point3d.h"
+
+namespace std{
+
+/// Output of array enclosed in brackets. Required for autocheck for output of
+/// std::arrays in case of failing test.
+template <typename T, std::size_t N>
+std::ostream& operator<<(std::ostream& os, std::array<T, N> const& array)
+{
+ os << "[";
+ std::copy(std::begin(array), std::end(array),
+ std::ostream_iterator<double>(os, ","));
+ os << "]";
+ return os;
+}
+}
+
namespace autocheck
{
diff --git a/Tests/Data b/Tests/Data
index 25646534c5e16759c301b6509764959929ae479f..9242c3d9125209ad5e585f8cd76f5d8a6315a002 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 25646534c5e16759c301b6509764959929ae479f
+Subproject commit 9242c3d9125209ad5e585f8cd76f5d8a6315a002
diff --git a/Tests/MaterialLib/KelvinVector.cpp b/Tests/MaterialLib/KelvinVector.cpp
index 9e148d48484590a301df719977c8f030f32095d3..88fad24936c8eebe1c8e8323e8ff400b28b96408 100644
--- a/Tests/MaterialLib/KelvinVector.cpp
+++ b/Tests/MaterialLib/KelvinVector.cpp
@@ -10,7 +10,7 @@
#include "MaterialLib/SolidModels/KelvinVector.h"
-#include "Tests/MathLib/AutoCheckTools.h"
+#include "Tests/AutoCheckTools.h"
using namespace MaterialLib::SolidModels;
namespace ac = autocheck;
diff --git a/Tests/MathLib/TestPoint3d.cpp b/Tests/MathLib/TestPoint3d.cpp
index 01e2fc85ffed00a9bb7fbf88b59b8eaff2abf29f..fd9b75b3feaa00018c2eb95519c34dcd63f7354e 100644
--- a/Tests/MathLib/TestPoint3d.cpp
+++ b/Tests/MathLib/TestPoint3d.cpp
@@ -16,7 +16,7 @@
#include "MathLib/Point3d.h"
-#include "AutoCheckTools.h"
+#include "Tests/AutoCheckTools.h"
using namespace MathLib;
namespace ac = autocheck;
diff --git a/Tests/MeshGeoToolsLib/TestGeoMapper.cpp b/Tests/MeshGeoToolsLib/TestGeoMapper.cpp
index 9b9209cf57d73f081b7e7175fcf6d1defec692a5..6098be1fb234398c5fdf8581e104ce38f1ec0669 100644
--- a/Tests/MeshGeoToolsLib/TestGeoMapper.cpp
+++ b/Tests/MeshGeoToolsLib/TestGeoMapper.cpp
@@ -9,14 +9,13 @@
#include <ctime>
#include <memory>
#include <gtest/gtest.h>
-#include <autocheck/autocheck.hpp>
#include "GeoLib/GEOObjects.h"
#include "GeoLib/Point.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshGenerators/MeshGenerator.h"
#include "MeshGeoToolsLib/GeoMapper.h"
-#include "../MathLib/AutoCheckTools.h"
+#include "Tests/AutoCheckTools.h"
namespace ac = autocheck;
diff --git a/Tests/NumLib/TestShapeFunctions.cpp b/Tests/NumLib/TestShapeFunctions.cpp
index 3f53a49bc44c8891ff40050593ff8393050f6984..37aedcb3926052ac38c89ac11bc7dc43fefc4d4a 100644
--- a/Tests/NumLib/TestShapeFunctions.cpp
+++ b/Tests/NumLib/TestShapeFunctions.cpp
@@ -11,16 +11,139 @@
#include <gtest/gtest.h>
+#include <algorithm>
+#include <type_traits>
#include <limits>
+#include <numeric>
#include <valarray>
-#include <algorithm>
+#include "MeshLib/Elements/Elements.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
-
+#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#include "Tests/AutoCheckTools.h"
#include "Tests/TestTools.h"
+namespace autocheck
+{
+
+template <typename ShapeFunction,
+ typename Gen = randomTupleGenerator<double, ShapeFunction::DIM>>
+struct NaturalPointGenerator
+{
+ Gen generator; // The default interval is fixed to [-1, 1].
+
+ using result_type = std::array<double, ShapeFunction::DIM>;
+
+ result_type intervalMap(result_type const& tuple) const
+ {
+ if (std::is_same<ShapeFunction, NumLib::ShapeLine2>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeLine3>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeQuad4>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeQuad8>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeQuad9>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeHex8>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeHex20>::value)
+ return tuple;
+ if (std::is_same<ShapeFunction, NumLib::ShapeTri3>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeTri6>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeTet4>::value ||
+ std::is_same<ShapeFunction, NumLib::ShapeTet10>::value)
+ {
+ // Map square (x, y) \in [-1, 1]^2 to a triangle such that x,y
+ // \in [0, 1], and x+y \in [0, 2/2].
+ // Same for three-d case with x + y + z \in [0, 3/2].
+ result_type mapped_tuple;
+ std::transform(std::begin(tuple), std::end(tuple),
+ std::begin(mapped_tuple),
+ [](double const& v) { return (v + 1.) / 2.; });
+
+ if (std::accumulate(std::begin(mapped_tuple),
+ std::end(mapped_tuple),
+ 0.) > ShapeFunction::DIM / 2.)
+ std::transform(std::begin(mapped_tuple), std::end(mapped_tuple),
+ std::begin(mapped_tuple),
+ [](double const& v) { return 1. - v; });
+ return mapped_tuple;
+ }
+ }
+
+ result_type operator()(std::size_t /*size*/ = 0)
+ {
+ return intervalMap(fix(1, generator)());
+ }
+};
+
+}
+
+namespace ac = autocheck;
using namespace NumLib;
+template <typename ShapeFunction>
+struct ShapeFunctionTest : public ::testing::Test
+{
+ ac::NaturalPointGenerator<ShapeFunction> natural_point_generator;
+ ac::gtest_reporter gtest_reporter;
+};
+
+using ShapeFunctionTestTypes =
+ ::testing::Types<ShapeLine2, ShapeLine3, ShapeTri3, ShapeTri6, ShapeTet4,
+ ShapeTet10, ShapeQuad4, ShapeQuad8, ShapeQuad9, ShapeHex8,
+ ShapeHex20>;
+
+TYPED_TEST_CASE(ShapeFunctionTest, ShapeFunctionTestTypes);
+
+// TypeParam is the type of the ShapeFunction.
+// Access private members via this pointer or TestFixture:: for types
+TYPED_TEST(ShapeFunctionTest, PartitionOfUnity)
+{
+ auto isPartitionOfUnity = [this](
+ std::array<double, TypeParam::DIM>& natural_coordinates_point) -> bool {
+
+ // compute shape functions
+ std::array<double, TypeParam::NPOINTS> N;
+ TypeParam::computeShapeFunction(natural_coordinates_point, N);
+ double const sum = std::accumulate(std::begin(N), std::end(N), 0.);
+
+ return std::abs(1. - sum) < 1e-15;
+ };
+
+ ac::check<std::array<double, TypeParam::DIM>>(
+ isPartitionOfUnity,
+ 1000,
+ ac::make_arbitrary(this->natural_point_generator),
+ this->gtest_reporter);
+}
+
+TYPED_TEST(ShapeFunctionTest, SumOfGradientsIsZero)
+{
+ auto isSumOfGradientsZero = [this](
+ std::array<double, TypeParam::DIM>& natural_coordinates_point) -> bool {
+
+ // compute shape functions
+ std::array<double, TypeParam::DIM * TypeParam::NPOINTS> dNdr;
+ TypeParam::computeGradShapeFunction(natural_coordinates_point, dNdr);
+ double const sum =
+ std::accumulate(std::begin(dNdr), std::end(dNdr), 0.);
+
+ return std::abs(sum) <= 5e-15;
+ };
+
+ ac::check<std::array<double, TypeParam::DIM>>(
+ isSumOfGradientsZero,
+ 1000,
+ ac::make_arbitrary(this->natural_point_generator),
+ this->gtest_reporter);
+}
+
TEST(NumLib, FemShapeQuad4)
{
static const double eps = std::numeric_limits<double>::epsilon();