diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestHex8.h b/Tests/NumLib/CoordinatesMappingTestData/TestHex8.h
new file mode 100644
index 0000000000000000000000000000000000000000..c164026df6e2508b622c99d975d9925e241f9cd9
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestHex8.h
@@ -0,0 +1,126 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTHEX8_H_
+#define COORDINATESMAPPINGTESTDATA_TESTHEX8_H_
+
+#include "MeshLib/Elements/Hex.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestHex8
+{
+ public:
+ // Element information
+ typedef MeshLib::Hex ElementType;
+ typedef NumLib::ShapeHex8 ShapeFunctionType;
+ static const unsigned dim = 3; //ElementType::dimension;
+ static const unsigned e_nnodes = ElementType::n_all_nodes;
+ // Coordinates where shape functions are evaluated
+ static const double r[dim];
+ // Expected results for natural shape
+ static const double nat_exp_N[e_nnodes];
+ static const double nat_exp_dNdr[e_nnodes*dim];
+ // Expected results for irregular shape
+ static const double ir_exp_J[dim*dim];
+ static const double ir_exp_invJ[dim*dim];
+ static const double ir_exp_detJ;
+ static const double ir_exp_dNdx[e_nnodes*dim];
+ // Expected results for clock-wise node ordering
+ static const double cl_exp_J[dim*dim];
+ // Expected results for zero volume
+ static const double cl_exp_detJ;
+ static const double ze_exp_J[dim*dim];
+
+ // element shape identical to that in natural coordinates (see ShapeHex8.h)
+ MeshLib::Hex* createNaturalShape()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
+ nodes[1] = new MeshLib::Node( 1.0, -1.0, -1.0);
+ nodes[2] = new MeshLib::Node( 1.0, 1.0, -1.0);
+ nodes[3] = new MeshLib::Node(-1.0, 1.0, -1.0);
+ nodes[4] = new MeshLib::Node(-1.0, -1.0, 1.0);
+ nodes[5] = new MeshLib::Node( 1.0, -1.0, 1.0);
+ nodes[6] = new MeshLib::Node( 1.0, 1.0, 1.0);
+ nodes[7] = new MeshLib::Node(-1.0, 1.0, 1.0);
+ return new MeshLib::Hex(nodes);
+ }
+
+ // element having irregular or skew shape
+ MeshLib::Hex* createIrregularShape()
+ {
+ // two times longer in z direction than the natural
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
+ nodes[1] = new MeshLib::Node( 1.0, -1.0, -1.0);
+ nodes[2] = new MeshLib::Node( 1.0, 1.0, -1.0);
+ nodes[3] = new MeshLib::Node(-1.0, 1.0, -1.0);
+ nodes[4] = new MeshLib::Node(-1.0, -1.0, 3.0);
+ nodes[5] = new MeshLib::Node( 1.0, -1.0, 3.0);
+ nodes[6] = new MeshLib::Node( 1.0, 1.0, 3.0);
+ nodes[7] = new MeshLib::Node(-1.0, 1.0, 3.0);
+ return new MeshLib::Hex(nodes);
+ }
+
+ // invalid case: clock wise node ordering
+ MeshLib::Hex* createClockWise()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
+ nodes[3] = new MeshLib::Node( 1.0, -1.0, -1.0);
+ nodes[2] = new MeshLib::Node( 1.0, 1.0, -1.0);
+ nodes[1] = new MeshLib::Node(-1.0, 1.0, -1.0);
+ nodes[4] = new MeshLib::Node(-1.0, -1.0, 1.0);
+ nodes[7] = new MeshLib::Node( 1.0, -1.0, 1.0);
+ nodes[6] = new MeshLib::Node( 1.0, 1.0, 1.0);
+ nodes[5] = new MeshLib::Node(-1.0, 1.0, 1.0);
+ return new MeshLib::Hex(nodes);
+ }
+
+ // invalid case: zero volume
+ MeshLib::Hex* createZeroVolume()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-1.0, -1.0, 1.0);
+ nodes[1] = new MeshLib::Node( 1.0, -1.0, 1.0);
+ nodes[2] = new MeshLib::Node( 1.0, 1.0, 1.0);
+ nodes[3] = new MeshLib::Node(-1.0, 1.0, 1.0);
+ nodes[4] = new MeshLib::Node(-1.0, -1.0, 1.0);
+ nodes[5] = new MeshLib::Node( 1.0, -1.0, 1.0);
+ nodes[6] = new MeshLib::Node( 1.0, 1.0, 1.0);
+ nodes[7] = new MeshLib::Node(-1.0, 1.0, 1.0);
+ return new MeshLib::Hex(nodes);
+ }
+};
+
+const double TestHex8::r[dim] = {0.5, 0.5, 0.5};
+const double TestHex8::nat_exp_N[e_nnodes]
+ = {0.015625, 0.046875, 0.140625, 0.046875, 0.046875, 0.140625, 0.421875, 0.140625};
+const double TestHex8::nat_exp_dNdr[e_nnodes*dim]
+ = {-0.03125, 0.03125, 0.09375, -0.09375, -0.09375, 0.09375, 0.28125, -0.28125,
+ -0.03125, -0.09375, 0.09375, 0.03125, -0.09375, -0.28125, 0.28125, 0.09375,
+ -0.03125, -0.09375, -0.28125, -0.09375, 0.03125, 0.09375, 0.28125, 0.09375};
+const double TestHex8::ir_exp_J[dim*dim]= {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0};
+const double TestHex8::ir_exp_detJ = 2.;
+const double TestHex8::ir_exp_invJ[dim*dim] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1./2.0};
+const double TestHex8::ir_exp_dNdx[dim*e_nnodes]
+ = {-0.03125, 0.03125, 0.09375, -0.09375, -0.09375, 0.09375, 0.28125, -0.28125,
+ -0.03125, -0.09375, 0.09375, 0.03125, -0.09375, -0.28125, 0.28125, 0.09375,
+ -0.015625, -0.046875, -0.140625, -0.046875, 0.015625, 0.046875, 0.140625, 0.046875};
+const double TestHex8::cl_exp_J[dim*dim] = {0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0};
+const double TestHex8::cl_exp_detJ = -1.;
+const double TestHex8::ze_exp_J[dim*dim] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0};
+
+}
+
+#endif //OGS_USE_EIGEN
+
+
diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestLine2.h b/Tests/NumLib/CoordinatesMappingTestData/TestLine2.h
new file mode 100644
index 0000000000000000000000000000000000000000..048380919c5ca7d314a9c8088aec266ab6a44829
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestLine2.h
@@ -0,0 +1,93 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTLINE2_H_
+#define COORDINATESMAPPINGTESTDATA_TESTLINE2_H_
+
+#include "MeshLib/Elements/Line.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestLine2
+{
+public:
+ // Element information
+ typedef MeshLib::Line ElementType;
+ typedef NumLib::ShapeLine2 ShapeFunctionType;
+ static const unsigned dim = ElementType::dimension;
+ static const unsigned e_nnodes = ElementType::n_all_nodes;
+ // Coordinates where shape functions are evaluated
+ static const double r[dim];
+ // Expected results for natural shape
+ static const double nat_exp_N[e_nnodes];
+ static const double nat_exp_dNdr[e_nnodes*dim];
+ // Expected results for irregular shape
+ static const double ir_exp_J[dim*dim];
+ static const double ir_exp_invJ[dim*dim];
+ static const double ir_exp_detJ;
+ static const double ir_exp_dNdx[e_nnodes*dim];
+ // Expected results for clock-wise node ordering
+ static const double cl_exp_J[dim*dim];
+ // Expected results for zero volume
+ static const double cl_exp_detJ;
+ static const double ze_exp_J[dim*dim];
+
+ // element shape identical to that in natural coordinates (see ShapeLine2.h)
+ static MeshLib::Line* createNaturalShape()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-1.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ // element having irregular or skew shape
+ MeshLib::Line* createIrregularShape()
+ {
+ // two times longer than the natural
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-2.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node( 2.0, 0.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ // invalid case: clock wise node ordering
+ MeshLib::Line* createClockWise()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[1] = new MeshLib::Node(-1.0, 0.0, 0.0);
+ nodes[0] = new MeshLib::Node( 1.0, 0.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ // invalid case: zero volume
+ MeshLib::Line* createZeroVolume()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(0.0, 0.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+};
+
+const double TestLine2::r[dim] = {0.5};
+const double TestLine2::nat_exp_N[e_nnodes] = {0.25, 0.75};
+const double TestLine2::nat_exp_dNdr[e_nnodes*dim] = {-0.5, 0.5};
+const double TestLine2::ir_exp_J[dim*dim] = {2.0};
+const double TestLine2::ir_exp_invJ[dim*dim] = {0.5};
+const double TestLine2::ir_exp_detJ = 2.;
+const double TestLine2::ir_exp_dNdx[dim*e_nnodes] = {-0.25, 0.25};
+const double TestLine2::cl_exp_J[dim*dim] = {-1.};
+const double TestLine2::cl_exp_detJ = -1;
+const double TestLine2::ze_exp_J[dim*dim] = {0.0};
+
+}
+
+#endif
diff --git a/Tests/NumLib/CoordinatesMappingTestData/TestQuad4.h b/Tests/NumLib/CoordinatesMappingTestData/TestQuad4.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2548cab02fccca6ea911c86ca8327cd565c8327
--- /dev/null
+++ b/Tests/NumLib/CoordinatesMappingTestData/TestQuad4.h
@@ -0,0 +1,105 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#ifndef COORDINATESMAPPINGTESTDATA_TESTQUAD4_H_
+#define COORDINATESMAPPINGTESTDATA_TESTQUAD4_H_
+
+#include "MeshLib/Elements/Quad.h"
+#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
+
+namespace CoordinatesMappingTestData
+{
+
+class TestQuad4
+{
+ public:
+ // Element information
+ typedef MeshLib::Quad ElementType;
+ typedef NumLib::ShapeQuad4 ShapeFunctionType;
+ static const unsigned dim = 2; //ElementType::dimension;
+ static const unsigned e_nnodes = ElementType::n_all_nodes;
+ // Coordinates where shape functions are evaluated
+ static const double r[dim];
+ // Expected results for natural shape
+ static const double nat_exp_N[e_nnodes];
+ static const double nat_exp_dNdr[e_nnodes*dim];
+ // Expected results for irregular shape
+ static const double ir_exp_J[dim*dim];
+ static const double ir_exp_invJ[dim*dim];
+ static const double ir_exp_detJ;
+ static const double ir_exp_dNdx[e_nnodes*dim];
+ // Expected results for clock-wise node ordering
+ static const double cl_exp_J[dim*dim];
+ // Expected results for zero volume
+ static const double cl_exp_detJ;
+ static const double ze_exp_J[dim*dim];
+
+ // element shape identical to that in natural coordinates
+ MeshLib::Quad* createNaturalShape()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 1.0, 1.0, 0.0);
+ nodes[1] = new MeshLib::Node(-1.0, 1.0, 0.0);
+ nodes[2] = new MeshLib::Node(-1.0, -1.0, 0.0);
+ nodes[3] = new MeshLib::Node( 1.0, -1.0, 0.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+ // element having irregular or skew shape
+ MeshLib::Quad* createIrregularShape()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(-0.5, -0.5, 0.0);
+ nodes[1] = new MeshLib::Node( 0.6, -0.6, 0.0);
+ nodes[2] = new MeshLib::Node( 0.5, 0.4, 0.0);
+ nodes[3] = new MeshLib::Node(-0.3, 0.1, 0.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+ // invalid case: clock wise node ordering
+ MeshLib::Quad* createClockWise()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 1.0, 1.0, 0.0);
+ nodes[3] = new MeshLib::Node(-1.0, 1.0, 0.0);
+ nodes[2] = new MeshLib::Node(-1.0, -1.0, 0.0);
+ nodes[1] = new MeshLib::Node( 1.0, -1.0, 0.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+ // invalid case: zero area
+ MeshLib::Quad* createZeroVolume()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 1.0, 1.0, 0.0);
+ nodes[1] = new MeshLib::Node(-1.0, 1.0, 0.0);
+ nodes[2] = new MeshLib::Node(-1.0, 1.0, 0.0);
+ nodes[3] = new MeshLib::Node( 1.0, 1.0, 0.0);
+ return new MeshLib::Quad(nodes);
+ }
+};
+
+const double TestQuad4::r[dim] = {0.5, 0.5};
+const double TestQuad4::nat_exp_N[e_nnodes] = {0.5625, 0.1875, 0.0625, 0.1875};
+const double TestQuad4::nat_exp_dNdr[e_nnodes*dim]
+ = {0.375, -0.375, -0.125, 0.125, 0.375, 0.125, -0.125, -0.375};
+const double TestQuad4::ir_exp_J[dim*dim] = {-0.5125, 0.0, -0.0625, -0.35};
+const double TestQuad4::ir_exp_detJ = 0.179375;
+const double TestQuad4::ir_exp_invJ[dim*dim]
+ = {-1.9512195121951219, 0.0, 0.3484320557491290, -2.8571428571428572};
+const double TestQuad4::ir_exp_dNdx[dim*e_nnodes]
+ = {-0.73170731707317072, 0.73170731707317072, 0.243902439024390, -0.24390243902439029,
+ -0.940766550522648, -0.48780487804878048, 0.313588850174216, 1.1149825783972125};
+const double TestQuad4::cl_exp_J[dim*dim] = {0.0, 1.0, 1.0, 0.0};
+const double TestQuad4::cl_exp_detJ = -1.;
+const double TestQuad4::ze_exp_J[dim*dim] = {1.0, 0.0, 0.0, 0.0};
+
+}
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/MatrixTools.h b/Tests/NumLib/FeTestData/MatrixTools.h
new file mode 100644
index 0000000000000000000000000000000000000000..b4f986c09c6856b1d68bf43303708205d6cef5ff
--- /dev/null
+++ b/Tests/NumLib/FeTestData/MatrixTools.h
@@ -0,0 +1,39 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATRIXTOOLS_H_
+#define MATRIXTOOLS_H_
+
+namespace FeTestData
+{
+
+// copy matrix entries in upper triangle to lower triangle
+template <class T_MATRIX, typename ID_TYPE=signed>
+inline void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m)
+{
+ for (ID_TYPE i=0; i<dim; i++)
+ for (ID_TYPE j=0; j<i; j++)
+ m(i,j) = m(j,i);
+}
+
+// set an identity matrix
+template <class T_MATRIX, typename ID_TYPE=signed>
+inline void setIdentityMatrix(unsigned dim, T_MATRIX &m)
+{
+ for (unsigned i=0; i<dim; i++)
+ for (unsigned j=0; j<dim; j++)
+ m(i,j) = 0.0;
+ for (unsigned i=0; i<dim; i++)
+ m(i,i) = 1.0;
+}
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/TestFeHEX8.h b/Tests/NumLib/FeTestData/TestFeHEX8.h
new file mode 100644
index 0000000000000000000000000000000000000000..e566087f9864633e5f98603f09d544a0d6c2b8fe
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeHEX8.h
@@ -0,0 +1,93 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef TESTFEHEX8_H_
+#define TESTFEHEX8_H_
+
+#include "MeshLib/Elements/Hex.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for HEX8
+class TestFeHEX8
+{
+public:
+ // Fe type information
+ template <class T_MATRIX_TYPES>
+ struct FeType
+ {
+ typedef NumLib::FeHEX8<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+ };
+ typedef MeshLib::Hex MeshElementType;
+ static const unsigned dim = 3; //MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+ static const unsigned n_sample_pt_order2 = 2*2*2;
+ static const unsigned n_sample_pt_order3 = 3*3*3;
+
+ /// create a mesh element
+ MeshLib::Hex* createMeshElement()
+ {
+ // cubic
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
+ nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
+ nodes[4] = new MeshLib::Node(0.0, 0.0, 1.0);
+ nodes[5] = new MeshLib::Node(1.0, 0.0, 1.0);
+ nodes[6] = new MeshLib::Node(1.0, 1.0, 1.0);
+ nodes[7] = new MeshLib::Node(0.0, 1.0, 1.0);
+ return new MeshLib::Hex(nodes);
+ }
+
+ /// set an expected mass matrix
+ template <class T_MATRIX>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2; m(0,4) = 1./2; m(0,5) = 1./4; m(0,6) = 1./8; m(0,7) = 1./4;
+ m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4; m(1,4) = 1./4; m(1,5) = 1./2; m(1,6) = 1./4; m(1,7) = 1./8;
+ m(2,2) = 1.0; m(2,3) = 1./2; m(2,4) = 1./8; m(2,5) = 1./4; m(2,6) = 1./2; m(2,7) = 1./4;
+ m(3,3) = 1.0; m(3,4) = 1./4; m(3,5) = 1./8; m(3,6) = 1./4; m(3,7) = 1./2;
+ m(4,4) = 1.0; m(4,5) = 1./2; m(4,6) = 1./4; m(4,7) = 1./2;
+ m(5,5) = 1.0; m(5,6) = 1./2; m(5,7) = 1./4;
+ m(6,6) = 1.0; m(6,7) = 1./2;
+ m(7,7) = 1.0;
+ // make symmetric
+ copyUpperToLower(e_nnodes, m);
+ m *= 1./27.;
+ }
+
+ /// set an expected laplace matrix
+ template <class T_MATRIX>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 2./3; m(0,1) = 0; m(0,2) = -1./6; m(0,3) = 0; m(0,4) = 0; m(0,5) = -1./6; m(0,6) = -1./6; m(0,7) = -1./6;
+ m(1,1) = 2./3; m(1,2) = 0; m(1,3) = -1./6; m(1,4) = -1./6; m(1,5) = 0; m(1,6) = -1./6; m(1,7) = -1./6;
+ m(2,2) = 2./3; m(2,3) = 0; m(2,4) = -1./6; m(2,5) = -1./6; m(2,6) = 0; m(2,7) = -1./6;
+ m(3,3) = 2./3; m(3,4) = -1./6; m(3,5) = -1./6; m(3,6) = -1./6; m(3,7) = 0;
+ m(4,4) = 2./3; m(4,5) = 0; m(4,6) = -1./6; m(4,7) = 0;
+ m(5,5) = 2./3; m(5,6) = 0; m(5,7) = -1./6;
+ m(6,6) = 2./3; m(6,7) = 0;
+ m(7,7) = 2./3;
+ // make symmetric
+ copyUpperToLower(e_nnodes, m);
+ m *= k/2.;
+ }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/TestFeLINE2.h b/Tests/NumLib/FeTestData/TestFeLINE2.h
new file mode 100644
index 0000000000000000000000000000000000000000..ce7b5e68d36dc79822a1591ace1c8211844cc982
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeLINE2.h
@@ -0,0 +1,73 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef TESTFELINE2_H_
+#define TESTFELINE2_H_
+
+#include "MeshLib/Elements/Line.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for LINE2
+class TestFeLINE2
+{
+public:
+ // Fe type information
+ template <class T_MATRIX_TYPES>
+ struct FeType
+ {
+ typedef NumLib::FeLINE2<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+ };
+ typedef MeshLib::Line MeshElementType;
+ static const unsigned dim = MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+ static const unsigned n_sample_pt_order2 = 2;
+ static const unsigned n_sample_pt_order3 = 3;
+
+ /// create a mesh element
+ MeshLib::Line* createMeshElement()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ /// set an expected mass matrix
+ template <class T_MATRIX>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1./3.; m(0,1) = 1./6.;
+ m(1,1) = 1./3.;
+ // make symmetric
+ copyUpperToLower(2, m);
+ }
+
+ /// set an expected laplace matrix
+ template <class T_MATRIX>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = -1.0;
+ m(1,1) = 1.0;
+ // make symmetric
+ copyUpperToLower(2, m);
+ m *= k;
+ }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/FeTestData/TestFeQUAD4.h b/Tests/NumLib/FeTestData/TestFeQUAD4.h
new file mode 100644
index 0000000000000000000000000000000000000000..009e531ad0afd3d1908240a683ec05410da49598
--- /dev/null
+++ b/Tests/NumLib/FeTestData/TestFeQUAD4.h
@@ -0,0 +1,81 @@
+/**
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef TESTFEQUAD4_H_
+#define TESTFEQUAD4_H_
+
+#include "MeshLib/Elements/Quad.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+
+#include "MatrixTools.h"
+
+namespace FeTestData
+{
+
+// Test case for QUAD4
+class TestFeQUAD4
+{
+public:
+ // Fe type information
+ template <class T_MATRIX_TYPES>
+ struct FeType
+ {
+ typedef NumLib::FeQUAD4<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
+ };
+ typedef MeshLib::Quad MeshElementType;
+ static const unsigned dim = 2; //MeshElementType::dimension;
+ static const unsigned e_nnodes = MeshElementType::n_all_nodes;
+ static const unsigned n_sample_pt_order2 = 2*2;
+ static const unsigned n_sample_pt_order3 = 3*3;
+
+ /// create a mesh element
+ MeshLib::Quad* createMeshElement()
+ {
+ // square
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
+ nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
+ nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+ /// set an expected mass matrix for 1m x 1m
+ template <class T_MATRIX>
+ void setExpectedMassMatrix(T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
+ m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
+ m(2,2) = 1.0; m(2,3) = 1./2;
+ m(3,3) = 1.0;
+ // make symmetric
+ copyUpperToLower(4, m);
+ m *= 1./9.;
+ }
+
+ /// set an expected laplace matrix for 1m x 1m
+ template <class T_MATRIX>
+ void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
+ {
+ // set upper triangle entries
+ m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
+ m(1,1) = 4.0; m(1,2) = -1.0; m(1,3) = -2.0;
+ m(2,2) = 4.0; m(2,3) = -1.0;
+ m(3,3) = 4.0;
+ // make symmetric
+ copyUpperToLower(4, m);
+ m *= k/6.;
+ }
+};
+
+} // namespace
+
+#endif
+
diff --git a/Tests/NumLib/TestCoordinatesMapping.cpp b/Tests/NumLib/TestCoordinatesMapping.cpp
index e3a881bc46cd5760e6289ab2b77eb77589530862..f4a5bf67f8985269b7a3478f84992a79418810e8 100644
--- a/Tests/NumLib/TestCoordinatesMapping.cpp
+++ b/Tests/NumLib/TestCoordinatesMapping.cpp
@@ -16,267 +16,29 @@
#include <Eigen/Eigen>
#endif
-#include "MeshLib/Elements/Line.h"
-#include "MeshLib/Elements/Quad.h"
-#include "MeshLib/Elements/Hex.h"
-#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
-#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
-#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h"
+#include "CoordinatesMappingTestData/TestLine2.h"
+#include "CoordinatesMappingTestData/TestQuad4.h"
+#include "CoordinatesMappingTestData/TestHex8.h"
+
#include "../TestTools.h"
using namespace NumLib;
+using namespace CoordinatesMappingTestData;
-#ifdef OGS_USE_EIGEN
namespace
{
-class Line2Test
-{
-public:
- typedef MeshLib::Line ElementType;
- typedef ShapeLine2 ShapeFunctionType;
- static const unsigned dim = ElementType::dimension;
- static const unsigned e_nnodes = ElementType::n_all_nodes;
- static const double r[dim];
- static const double nat_exp_N[e_nnodes];
- static const double nat_exp_dNdr[e_nnodes*dim];
- static const double ir_exp_J[dim*dim];
- static const double ir_exp_invJ[dim*dim];
- static const double ir_exp_detJ;
- static const double ir_exp_dNdx[e_nnodes*dim];
- static const double cl_exp_J[dim*dim];
- static const double cl_exp_detJ;
- static const double ze_exp_J[dim*dim];
-
- // element shape identical to that in natural coordinates (see ShapeLine2.h)
- static MeshLib::Line* createNaturalShape()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-1.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node( 1.0, 0.0, 0.0);
- return new MeshLib::Line(nodes);
- }
-
- // element having irregular or skew shape
- MeshLib::Line* createIrregularShape()
- {
- // two times longer than the natural
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-2.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node( 2.0, 0.0, 0.0);
- return new MeshLib::Line(nodes);
- }
-
- // invalid case: clock wise node ordering
- MeshLib::Line* createClockWise()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[1] = new MeshLib::Node(-1.0, 0.0, 0.0);
- nodes[0] = new MeshLib::Node( 1.0, 0.0, 0.0);
- return new MeshLib::Line(nodes);
- }
-
- // invalid case: zero volume
- MeshLib::Line* createZeroVolume()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node(0.0, 0.0, 0.0);
- return new MeshLib::Line(nodes);
- }
-};
-
-const double Line2Test::r[dim] = {0.5};
-const double Line2Test::nat_exp_N[e_nnodes] = {0.25, 0.75};
-const double Line2Test::nat_exp_dNdr[e_nnodes*dim] = {-0.5, 0.5};
-const double Line2Test::ir_exp_J[dim*dim]= {2.0};
-const double Line2Test::ir_exp_invJ[dim*dim]= {0.5};
-const double Line2Test::ir_exp_detJ = 2.;
-const double Line2Test::ir_exp_dNdx[dim*e_nnodes]= {-0.25, 0.25};
-const double Line2Test::cl_exp_J[dim*dim]= {-1.};
-const double Line2Test::cl_exp_detJ = -1;
-const double Line2Test::ze_exp_J[dim*dim]= {0.0};
-
-class Quad4Test
-{
- public:
- // Matrix types
- typedef MeshLib::Quad ElementType;
- typedef ShapeQuad4 ShapeFunctionType;
- static const unsigned dim = 2; //ElementType::dimension;
- static const unsigned e_nnodes = ElementType::n_all_nodes;
-
- // quad having shape identical to that in natural coordinates
- MeshLib::Quad* createNaturalShape()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node( 1.0, 1.0, 0.0);
- nodes[1] = new MeshLib::Node(-1.0, 1.0, 0.0);
- nodes[2] = new MeshLib::Node(-1.0, -1.0, 0.0);
- nodes[3] = new MeshLib::Node( 1.0, -1.0, 0.0);
- return new MeshLib::Quad(nodes);
- }
-
- // quad having irregular or skew shape
- MeshLib::Quad* createIrregularShape()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-0.5, -0.5, 0.0);
- nodes[1] = new MeshLib::Node( 0.6, -0.6, 0.0);
- nodes[2] = new MeshLib::Node( 0.5, 0.4, 0.0);
- nodes[3] = new MeshLib::Node(-0.3, 0.1, 0.0);
- return new MeshLib::Quad(nodes);
- }
-
- // invalid case: clock wise node ordering
- MeshLib::Quad* createClockWise()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node( 1.0, 1.0, 0.0);
- nodes[3] = new MeshLib::Node(-1.0, 1.0, 0.0);
- nodes[2] = new MeshLib::Node(-1.0, -1.0, 0.0);
- nodes[1] = new MeshLib::Node( 1.0, -1.0, 0.0);
- return new MeshLib::Quad(nodes);
- }
-
- // invalid case: zero area
- MeshLib::Quad* createZeroVolume()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node( 1.0, 1.0, 0.0);
- nodes[1] = new MeshLib::Node(-1.0, 1.0, 0.0);
- nodes[2] = new MeshLib::Node(-1.0, 1.0, 0.0);
- nodes[3] = new MeshLib::Node( 1.0, 1.0, 0.0);
- return new MeshLib::Quad(nodes);
- }
-
- static const double r[dim];
- static const double nat_exp_N[e_nnodes];
- static const double nat_exp_dNdr[e_nnodes*dim];
- static const double ir_exp_J[dim*dim];
- static const double ir_exp_invJ[dim*dim];
- static const double ir_exp_detJ;
- static const double ir_exp_dNdx[e_nnodes*dim];
- static const double cl_exp_J[dim*dim];
- static const double cl_exp_detJ;
- static const double ze_exp_J[dim*dim];
-};
-
-const double Quad4Test::r[dim] = {0.5, 0.5};
-const double Quad4Test::nat_exp_N[e_nnodes] = {0.5625, 0.1875, 0.0625, 0.1875};
-const double Quad4Test::nat_exp_dNdr[e_nnodes*dim] = {0.375, -0.375, -0.125, 0.125, 0.375, 0.125, -0.125, -0.375};
-const double Quad4Test::ir_exp_J[dim*dim]= {-0.5125, 0.0, -0.0625, -0.35};
-const double Quad4Test::ir_exp_detJ= 0.179375;
-const double Quad4Test::ir_exp_invJ[dim*dim]= {-1.9512195121951219, 0.0, 0.3484320557491290, -2.8571428571428572};
-const double Quad4Test::ir_exp_dNdx[dim*e_nnodes]= {-0.73170731707317072, 0.73170731707317072, 0.243902439024390, -0.24390243902439029, -0.940766550522648, -0.48780487804878048, 0.313588850174216, 1.1149825783972125};
-const double Quad4Test::cl_exp_J[dim*dim]= {0.0, 1.0, 1.0, 0.0};
-const double Quad4Test::cl_exp_detJ= -1.;
-const double Quad4Test::ze_exp_J[dim*dim]= {1.0, 0.0, 0.0, 0.0};
-
-class Hex8Test
-{
- public:
- // Matrix types
- typedef MeshLib::Hex ElementType;
- typedef ShapeHex8 ShapeFunctionType;
- static const unsigned dim = 3; //ElementType::dimension;
- static const unsigned e_nnodes = ElementType::n_all_nodes;
-
- // element shape identical to that in natural coordinates (see ShapeHex8.h)
- MeshLib::Hex* createNaturalShape()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
- nodes[1] = new MeshLib::Node( 1.0, -1.0, -1.0);
- nodes[2] = new MeshLib::Node( 1.0, 1.0, -1.0);
- nodes[3] = new MeshLib::Node(-1.0, 1.0, -1.0);
- nodes[4] = new MeshLib::Node(-1.0, -1.0, 1.0);
- nodes[5] = new MeshLib::Node( 1.0, -1.0, 1.0);
- nodes[6] = new MeshLib::Node( 1.0, 1.0, 1.0);
- nodes[7] = new MeshLib::Node(-1.0, 1.0, 1.0);
- return new MeshLib::Hex(nodes);
- }
-
- // element having irregular or skew shape
- MeshLib::Hex* createIrregularShape()
- {
- // two times longer in z direction than the natural
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
- nodes[1] = new MeshLib::Node( 1.0, -1.0, -1.0);
- nodes[2] = new MeshLib::Node( 1.0, 1.0, -1.0);
- nodes[3] = new MeshLib::Node(-1.0, 1.0, -1.0);
- nodes[4] = new MeshLib::Node(-1.0, -1.0, 3.0);
- nodes[5] = new MeshLib::Node( 1.0, -1.0, 3.0);
- nodes[6] = new MeshLib::Node( 1.0, 1.0, 3.0);
- nodes[7] = new MeshLib::Node(-1.0, 1.0, 3.0);
- return new MeshLib::Hex(nodes);
- }
-
- // invalid case: clock wise node ordering
- MeshLib::Hex* createClockWise()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-1.0, -1.0, -1.0);
- nodes[3] = new MeshLib::Node( 1.0, -1.0, -1.0);
- nodes[2] = new MeshLib::Node( 1.0, 1.0, -1.0);
- nodes[1] = new MeshLib::Node(-1.0, 1.0, -1.0);
- nodes[4] = new MeshLib::Node(-1.0, -1.0, 1.0);
- nodes[7] = new MeshLib::Node( 1.0, -1.0, 1.0);
- nodes[6] = new MeshLib::Node( 1.0, 1.0, 1.0);
- nodes[5] = new MeshLib::Node(-1.0, 1.0, 1.0);
- return new MeshLib::Hex(nodes);
- }
-
- // invalid case: zero volume
- MeshLib::Hex* createZeroVolume()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(-1.0, -1.0, 1.0);
- nodes[1] = new MeshLib::Node( 1.0, -1.0, 1.0);
- nodes[2] = new MeshLib::Node( 1.0, 1.0, 1.0);
- nodes[3] = new MeshLib::Node(-1.0, 1.0, 1.0);
- nodes[4] = new MeshLib::Node(-1.0, -1.0, 1.0);
- nodes[5] = new MeshLib::Node( 1.0, -1.0, 1.0);
- nodes[6] = new MeshLib::Node( 1.0, 1.0, 1.0);
- nodes[7] = new MeshLib::Node(-1.0, 1.0, 1.0);
- return new MeshLib::Hex(nodes);
- }
-
- static const double r[dim];
- static const double nat_exp_N[e_nnodes];
- static const double nat_exp_dNdr[e_nnodes*dim];
- static const double ir_exp_J[dim*dim];
- static const double ir_exp_invJ[dim*dim];
- static const double ir_exp_detJ;
- static const double ir_exp_dNdx[e_nnodes*dim];
- static const double cl_exp_J[dim*dim];
- static const double cl_exp_detJ;
- static const double ze_exp_J[dim*dim];
-};
-
-const double Hex8Test::r[dim] = {0.5, 0.5, 0.5};
-const double Hex8Test::nat_exp_N[e_nnodes]
- = {0.015625, 0.046875, 0.140625, 0.046875, 0.046875, 0.140625, 0.421875, 0.140625};
-const double Hex8Test::nat_exp_dNdr[e_nnodes*dim]
- = {-0.03125, 0.03125, 0.09375, -0.09375, -0.09375, 0.09375, 0.28125, -0.28125,
- -0.03125, -0.09375, 0.09375, 0.03125, -0.09375, -0.28125, 0.28125, 0.09375,
- -0.03125, -0.09375, -0.28125, -0.09375, 0.03125, 0.09375, 0.28125, 0.09375};
-const double Hex8Test::ir_exp_J[dim*dim]= {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 2.0};
-const double Hex8Test::ir_exp_detJ= 2.;
-const double Hex8Test::ir_exp_invJ[dim*dim]= {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1./2.0};
-const double Hex8Test::ir_exp_dNdx[dim*e_nnodes]
- = {-0.03125, 0.03125, 0.09375, -0.09375, -0.09375, 0.09375, 0.28125, -0.28125,
- -0.03125, -0.09375, 0.09375, 0.03125, -0.09375, -0.28125, 0.28125, 0.09375,
- -0.015625, -0.046875, -0.140625, -0.046875, 0.015625, 0.046875, 0.140625, 0.046875};
-const double Hex8Test::cl_exp_J[dim*dim]= {0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0};
-const double Hex8Test::cl_exp_detJ= -1.;
-const double Hex8Test::ze_exp_J[dim*dim]= {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0};
-
+// Element types to be tested
+typedef ::testing::Types<
+ TestLine2,
+ TestQuad4,
+ TestHex8
+ > TestTypes;
+}
template <class T_TEST>
class NumLibFemNaturalCoordinatesMappingTest : public ::testing::Test, public T_TEST
@@ -287,16 +49,19 @@ public:
static const unsigned dim = T_TEST::dim;
static const unsigned e_nnodes = T_TEST::e_nnodes;
// Matrix types
+#ifdef OGS_USE_EIGEN
typedef Eigen::Matrix<double, e_nnodes, 1> NodalVector;
typedef Eigen::Matrix<double, dim, e_nnodes, Eigen::RowMajor> DimNodalMatrix;
typedef Eigen::Matrix<double, dim, dim, Eigen::RowMajor> DimMatrix;
+#endif
// Shape data type
typedef ShapeMatrices<NodalVector,DimNodalMatrix,DimMatrix> ShapeMatricesType;
// Natural coordinates mapping type
typedef NaturalCoordinatesMapping<ElementType, ShapeFunctionType, ShapeMatricesType> NaturalCoordsMappingType;
public:
- NumLibFemNaturalCoordinatesMappingTest() : eps(std::numeric_limits<double>::epsilon())
+ NumLibFemNaturalCoordinatesMappingTest()
+ : eps(std::numeric_limits<double>::epsilon())
{
// create four elements used for testing
naturalEle = this->createNaturalShape();
@@ -331,15 +96,7 @@ public:
ElementType* zeroVolumeEle;
};
-} // namespace
-
-typedef ::testing::Types<
- Line2Test,
- Quad4Test,
- Hex8Test
- > ElementTypes;
-
-TYPED_TEST_CASE(NumLibFemNaturalCoordinatesMappingTest, ElementTypes);
+TYPED_TEST_CASE(NumLibFemNaturalCoordinatesMappingTest, TestTypes);
TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_N)
{
@@ -468,7 +225,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckIrregularShape)
// irregular shape
ShapeMatricesType shape(this->dim, this->e_nnodes);
NaturalCoordsMappingType::computeShapeMatrices(*this->irregularEle, this->r, shape);
- //std::cout << shape;
+ //std::cout << std::setprecision(16) << shape;
ASSERT_ARRAY_NEAR(this->nat_exp_N, shape.N.data(), shape.N.size(), this->eps);
ASSERT_ARRAY_NEAR(this->nat_exp_dNdr, shape.dNdr.data(), shape.dNdr.size(), this->eps);
@@ -486,7 +243,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckClockwise)
// clockwise node ordering, which is invalid)
ShapeMatricesType shape(this->dim, this->e_nnodes);
NaturalCoordsMappingType::computeShapeMatrices(*this->clockwiseEle, this->r, shape);
- //std::cout << shape;
+ //std::cout << std::setprecision(16) << shape;
// Inverse of the Jacobian matrix doesn't exist
double exp_invJ[TestFixture::dim*TestFixture::dim]= {0.0};
double exp_dNdx[TestFixture::dim*TestFixture::e_nnodes]= {0.0};
@@ -506,7 +263,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckZeroVolume)
ShapeMatricesType shape(this->dim, this->e_nnodes);
NaturalCoordsMappingType::computeShapeMatrices(*this->zeroVolumeEle, this->r, shape);
- //std::cout << shape;
+ //std::cout << std::setprecision(16) << shape;
// Inverse of the Jacobian matrix doesn't exist
double exp_invJ[TestFixture::dim*TestFixture::dim]= {0.0};
double exp_dNdx[TestFixture::dim*TestFixture::e_nnodes]= {0.0};
@@ -519,6 +276,4 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckZeroVolume)
ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), this->eps);
}
-#endif //OGS_USE_EIGEN
-
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index c2bccba1a31b8b14bfc3ba494b0f1923c9286edc..e53e50dfab05993efc0a24c87f978aafc1620726 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -15,200 +15,61 @@
#include <Eigen/Eigen>
#endif
-#include "MeshLib/Elements/Line.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+#include "FeTestData/TestFeLINE2.h"
+#include "FeTestData/TestFeQUAD4.h"
+#include "FeTestData/TestFeHEX8.h"
+
#include "Tests/TestTools.h"
using namespace NumLib;
+using namespace FeTestData;
namespace
{
-
-// copy matrix entries in upper triangle to lower triangle
-template <class T_MATRIX, typename ID_TYPE=signed>
-void copyUpperToLower(const ID_TYPE dim, T_MATRIX &m)
-{
- for (ID_TYPE i=0; i<dim; i++)
- for (ID_TYPE j=0; j<i; j++)
- m(i,j) = m(j,i);
-}
-
-// set an identity matrix
-template <class T_MATRIX, typename ID_TYPE=signed>
-void setIdentityMatrix(unsigned dim, T_MATRIX &m)
-{
- for (unsigned i=0; i<dim; i++)
- for (unsigned j=0; j<dim; j++)
- m(i,j) = 0.0;
- for (unsigned i=0; i<dim; i++)
- m(i,i) = 1.0;
-}
-
-// Test case for LINE2
-class TestFeLINE2
+#ifdef OGS_USE_EIGEN
+template <typename T_FE>
+struct EigenFixedMatrixTypes
{
-public:
- template <class T_MATRIX_TYPES>
- struct FeType
- {
- typedef NumLib::FeLINE2<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
- };
- typedef MeshLib::Line MeshElementType;
- static const unsigned dim = MeshElementType::dimension;
- static const unsigned e_nnodes = MeshElementType::n_all_nodes;
-
- MeshLib::Line* createMeshElement()
- {
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
- return new MeshLib::Line(nodes);
- }
-
- // set an expected mass matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedMassMatrix(T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1./3.; m(0,1) = 1./6.;
- m(1,1) = 1./3.;
- // make symmetric
- copyUpperToLower(2, m);
- }
-
- // set an expected laplace matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1.0; m(0,1) = -1.0;
- m(1,1) = 1.0;
- // make symmetric
- copyUpperToLower(2, m);
- m *= k;
- }
+ typedef Eigen::Matrix<double, T_FE::e_nnodes, T_FE::e_nnodes, Eigen::RowMajor> NodalMatrixType;
+ typedef Eigen::Matrix<double, T_FE::e_nnodes, 1> NodalVectorType;
+ typedef Eigen::Matrix<double, T_FE::dim, T_FE::e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
+ typedef Eigen::Matrix<double, T_FE::dim, T_FE::dim, Eigen::RowMajor> DimMatrixType;
};
-// Test case for QUAD4
-class TestFeQUAD4
+struct EigenDynamicMatrixTypes
{
-public:
- template <class T_MATRIX_TYPES>
- struct FeType
- {
- typedef NumLib::FeQUAD4<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
- };
- typedef MeshLib::Quad MeshElementType;
- static const unsigned dim = 2; //MeshElementType::dimension;
- static const unsigned e_nnodes = MeshElementType::n_all_nodes;
-
- MeshLib::Quad* createMeshElement()
- {
- // square
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
- nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
- nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
- return new MeshLib::Quad(nodes);
- }
-
- // set an expected mass matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedMassMatrix(T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2;
- m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4;
- m(2,2) = 1.0; m(2,3) = 1./2;
- m(3,3) = 1.0;
- // make symmetric
- copyUpperToLower(4, m);
- m *= 1./9.;
- }
-
- // set an expected laplace matrix for 1m x 1m
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 4.0; m(0,1) = -1.0; m(0,2) = -2.0; m(0,3) = -1.0;
- m(1,1) = 4.0; m(1,2) = -1.0; m(1,3) = -2.0;
- m(2,2) = 4.0; m(2,3) = -1.0;
- m(3,3) = 4.0;
- // make symmetric
- copyUpperToLower(4, m);
- m *= k/6.;
- }
+ typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
+ typedef NodalMatrixType DimNodalMatrixType;
+ typedef NodalMatrixType DimMatrixType;
+ typedef Eigen::VectorXd NodalVectorType;
};
-// Test case for HEX8
-class TestFeHEX8
-{
-public:
- template <class T_MATRIX_TYPES>
- struct FeType
- {
- typedef NumLib::FeHEX8<typename T_MATRIX_TYPES::NodalVectorType, typename T_MATRIX_TYPES::DimNodalMatrixType, typename T_MATRIX_TYPES::DimMatrixType> type;
- };
- typedef MeshLib::Hex MeshElementType;
- static const unsigned dim = 3; //MeshElementType::dimension;
- static const unsigned e_nnodes = MeshElementType::n_all_nodes;
-
- MeshLib::Hex* createMeshElement()
- {
- // cubic
- MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
- nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
- nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0);
- nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0);
- nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0);
- nodes[4] = new MeshLib::Node(0.0, 0.0, 1.0);
- nodes[5] = new MeshLib::Node(1.0, 0.0, 1.0);
- nodes[6] = new MeshLib::Node(1.0, 1.0, 1.0);
- nodes[7] = new MeshLib::Node(0.0, 1.0, 1.0);
- return new MeshLib::Hex(nodes);
- }
-
- // set an expected mass matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedMassMatrix(T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 1.0; m(0,1) = 1./2; m(0,2) = 1./4; m(0,3) = 1./2; m(0,4) = 1./2; m(0,5) = 1./4; m(0,6) = 1./8; m(0,7) = 1./4;
- m(1,1) = 1.0; m(1,2) = 1./2; m(1,3) = 1./4; m(1,4) = 1./4; m(1,5) = 1./2; m(1,6) = 1./4; m(1,7) = 1./8;
- m(2,2) = 1.0; m(2,3) = 1./2; m(2,4) = 1./8; m(2,5) = 1./4; m(2,6) = 1./2; m(2,7) = 1./4;
- m(3,3) = 1.0; m(3,4) = 1./4; m(3,5) = 1./8; m(3,6) = 1./4; m(3,7) = 1./2;
- m(4,4) = 1.0; m(4,5) = 1./2; m(4,6) = 1./4; m(4,7) = 1./2;
- m(5,5) = 1.0; m(5,6) = 1./2; m(5,7) = 1./4;
- m(6,6) = 1.0; m(6,7) = 1./2;
- m(7,7) = 1.0;
- // make symmetric
- copyUpperToLower(e_nnodes, m);
- m *= 1./27.;
- }
+typedef EigenDynamicMatrixTypes DefaultMatrixType;
+#endif // OGS_USE_EIGEN
- // set an expected laplace matrix
- template <class T_MATRIX, typename ID_TYPE=signed>
- void setExpectedLaplaceMatrix(double k, T_MATRIX &m)
- {
- // set upper triangle entries
- m(0,0) = 2./3; m(0,1) = 0; m(0,2) = -1./6; m(0,3) = 0; m(0,4) = 0; m(0,5) = -1./6; m(0,6) = -1./6; m(0,7) = -1./6;
- m(1,1) = 2./3; m(1,2) = 0; m(1,3) = -1./6; m(1,4) = -1./6; m(1,5) = 0; m(1,6) = -1./6; m(1,7) = -1./6;
- m(2,2) = 2./3; m(2,3) = 0; m(2,4) = -1./6; m(2,5) = -1./6; m(2,6) = 0; m(2,7) = -1./6;
- m(3,3) = 2./3; m(3,4) = -1./6; m(3,5) = -1./6; m(3,6) = -1./6; m(3,7) = 0;
- m(4,4) = 2./3; m(4,5) = 0; m(4,6) = -1./6; m(4,7) = 0;
- m(5,5) = 2./3; m(5,6) = 0; m(5,7) = -1./6;
- m(6,6) = 2./3; m(6,7) = 0;
- m(7,7) = 2./3;
- // make symmetric
- copyUpperToLower(e_nnodes, m);
- m *= k/2.;
- }
+// test cases
+template <class T_FE_TYPE, class T_MAT_TYPES=DefaultMatrixType>
+struct TestCase
+{
+ typedef T_FE_TYPE T_FE;
+ typedef T_MAT_TYPES T_MATRIX_TYPES;
};
+typedef ::testing::Types<
+ TestCase<TestFeLINE2>,
+ TestCase<TestFeQUAD4>,
+ TestCase<TestFeHEX8>
+#ifdef OGS_USE_EIGEN
+ ,TestCase<TestFeLINE2, EigenFixedMatrixTypes<TestFeLINE2> >
+ ,TestCase<TestFeQUAD4, EigenFixedMatrixTypes<TestFeQUAD4> >
+ ,TestCase<TestFeHEX8, EigenFixedMatrixTypes<TestFeHEX8> >
+#endif
+ > TestTypes;
+}
+
template <class T>
class NumLibFemIsoTest : public ::testing::Test, public T::T_FE
{
@@ -229,6 +90,8 @@ class NumLibFemIsoTest : public ::testing::Test, public T::T_FE
static const unsigned dim = T_FE::dim;
static const unsigned e_nnodes = T_FE::e_nnodes;
+ static const unsigned n_sample_pt_order2 = T_FE::n_sample_pt_order2;
+ static const unsigned n_sample_pt_order3 = T_FE::n_sample_pt_order3;
public:
NumLibFemIsoTest() :
@@ -289,48 +152,11 @@ const unsigned NumLibFemIsoTest<T>::dim;
template <class T>
const unsigned NumLibFemIsoTest<T>::e_nnodes;
-} // namespace
-
-#ifdef OGS_USE_EIGEN
-template <typename T_FE>
-struct EigenFixedMatrixTypes
-{
- typedef Eigen::Matrix<double, T_FE::e_nnodes, T_FE::e_nnodes, Eigen::RowMajor> NodalMatrixType;
- typedef Eigen::Matrix<double, T_FE::e_nnodes, 1> NodalVectorType;
- typedef Eigen::Matrix<double, T_FE::dim, T_FE::e_nnodes, Eigen::RowMajor> DimNodalMatrixType;
- typedef Eigen::Matrix<double, T_FE::dim, T_FE::dim, Eigen::RowMajor> DimMatrixType;
-};
-
-struct EigenDynamicMatrixTypes
-{
- typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> NodalMatrixType;
- typedef NodalMatrixType DimNodalMatrixType;
- typedef NodalMatrixType DimMatrixType;
- typedef Eigen::VectorXd NodalVectorType;
-};
-
-typedef EigenDynamicMatrixTypes DefaultMatrixType;
-#endif // OGS_USE_EIGEN
-
-// test cases
-template <class T_FE_TYPE, class T_MAT_TYPES=DefaultMatrixType>
-struct TestCase
-{
- typedef T_FE_TYPE T_FE;
- typedef T_MAT_TYPES T_MATRIX_TYPES;
-};
-
-typedef ::testing::Types<
- TestCase<TestFeLINE2>,
- TestCase<TestFeQUAD4>,
- TestCase<TestFeHEX8>
-#ifdef OGS_USE_EIGEN
- ,TestCase<TestFeLINE2, EigenFixedMatrixTypes<TestFeLINE2> >
- ,TestCase<TestFeQUAD4, EigenFixedMatrixTypes<TestFeQUAD4> >
- ,TestCase<TestFeHEX8, EigenFixedMatrixTypes<TestFeHEX8> >
-#endif
- > TestTypes;
+template <class T>
+const unsigned NumLibFemIsoTest<T>::n_sample_pt_order2;
+template <class T>
+const unsigned NumLibFemIsoTest<T>::n_sample_pt_order3;
TYPED_TEST_CASE(NumLibFemIsoTest, TestTypes);
@@ -417,7 +243,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckGaussIntegrationLevel)
// evaluate a mass matrix
NodalMatrix M(this->e_nnodes, this->e_nnodes);
ShapeMatricesType shape(this->dim, this->e_nnodes);
- ASSERT_EQ(std::pow(2u, this->dim), this->integration_method.getNPoints());
+ ASSERT_EQ(TestFixture::n_sample_pt_order2, this->integration_method.getNPoints());
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
@@ -429,7 +255,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckGaussIntegrationLevel)
// Change gauss quadrature level to 3
this->integration_method.setIntegrationOrder(3);
M *= .0;
- ASSERT_EQ(std::pow(3u, this->dim), this->integration_method.getNPoints());
+ ASSERT_EQ(TestFixture::n_sample_pt_order3, this->integration_method.getNPoints());
for (std::size_t i=0; i < this->integration_method.getNPoints(); i++) {
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);