diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 62e3dab9324d15308cdee923a8ec438630fd90b3..1a7ae23ebbe98fbfeab6144ad3a94b04f292e07d 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -124,9 +124,31 @@ if(NOT OGS_USE_MPI)
ABSTOL 1e-14 RELTOL 1e-14
DIFF_DATA line_${mesh_size}_neumann_pcs_0_ts_1_t_1.000000.vtu D1_left_N1_right pressure
)
- endforeach()
+ endforeach()
+ # Some Neumann BC tests
+ AddTest(
+ NAME GroundWaterFlowProcess_cube_top
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e3_top_neumann.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1e-16 RELTOL 1e-16
+ DIFF_DATA cube_1e3_top_neumann_pcs_0_ts_1_t_1.000000.vtu pressure pressure
+ )
+ AddTest(
+ NAME GroundWaterFlowProcess_cube_bottom
+ PATH Elliptic/cube_1x1x1_GroundWaterFlow
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS cube_1e3_bottom_neumann.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1e-16 RELTOL 1e-16
+ DIFF_DATA cube_1e3_bottom_neumann_pcs_0_ts_1_t_1.000000.vtu pressure pressure
+ )
+ # TES tests
AddTest(
NAME TES_zeolite_discharge_small
PATH Parabolic/TES/1D
diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
index a0e159bbef67fc97de4eea77e006c5f7b1a0c8db..758faef0607d32bd561abbf4fbe478ce6b95398d 100644
--- a/GeoLib/AnalyticalGeometry-impl.h
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -118,8 +118,7 @@ void computeRotationMatrixToXY(MathLib::Vector3 const& n,
T_MATRIX & rot_mat)
{
// check if normal points already in the right direction
- if (sqrt(n[0]*n[0]+n[1]*n[1]) == 0) {
- rot_mat(0,0) = 1.0;
+ if (n[0] == 0 && n[1] == 0) {
rot_mat(0,1) = 0.0;
rot_mat(0,2) = 0.0;
rot_mat(1,0) = 0.0;
@@ -127,7 +126,17 @@ void computeRotationMatrixToXY(MathLib::Vector3 const& n,
rot_mat(1,2) = 0.0;
rot_mat(2,0) = 0.0;
rot_mat(2,1) = 0.0;
- rot_mat(2,2) = 1.0;
+
+ if (n[2] > 0) {
+ // identity matrix
+ rot_mat(0,0) = 1.0;
+ rot_mat(2,2) = 1.0;
+ } else {
+ // rotate by pi about the y-axis
+ rot_mat(0,0) = -1.0;
+ rot_mat(2,2) = -1.0;
+ }
+
return;
}
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
index 1c79430075606b19456f78f78de056255de06727..4c4bfbb15c7d60d9fb92073346324b245fbc344c 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.cpp
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -68,24 +68,23 @@ namespace MeshLib
ElementCoordinatesMappingLocal::ElementCoordinatesMappingLocal(
const Element& e,
- const CoordinateSystem &global_coords)
-: _coords(global_coords), _matR2global(3,3)
+ const unsigned global_dim)
+: _global_dim(global_dim), _matR2global(3,3)
{
- assert(e.getDimension() <= global_coords.getDimension());
+ assert(e.getDimension() <= global_dim);
_points.reserve(e.getNumberOfNodes());
for(unsigned i = 0; i < e.getNumberOfNodes(); i++)
_points.emplace_back(*(e.getNode(i)));
- auto const element_dimension = e.getDimension();
- auto const global_dimension = global_coords.getDimension();
+ auto const element_dim = e.getDimension();
- if (global_dimension == element_dimension)
+ if (global_dim == element_dim)
{
_matR2global.setIdentity();
return;
}
- detail::getRotationMatrixToGlobal(element_dimension, global_dimension, _points, _matR2global);
+ detail::getRotationMatrixToGlobal(element_dim, global_dim, _points, _matR2global);
detail::rotateToLocal(_matR2global.transpose(), _points);
}
diff --git a/MeshLib/ElementCoordinatesMappingLocal.h b/MeshLib/ElementCoordinatesMappingLocal.h
index c1d74d07000160b33c0896b9d743c213c962e317..def36698f1e1840c2f98aa8c29becaeb4482ae2a 100644
--- a/MeshLib/ElementCoordinatesMappingLocal.h
+++ b/MeshLib/ElementCoordinatesMappingLocal.h
@@ -10,13 +10,9 @@
#define ELEMENTCOORDINATESMAPPINGLOCAL_H_
#include <vector>
-
#include <Eigen/Eigen>
-
#include "MathLib/Point3d.h"
-#include "MeshLib/CoordinateSystem.h"
-
namespace MeshLib
{
class Element;
@@ -34,13 +30,13 @@ class ElementCoordinatesMappingLocal final
public:
/**
* Constructor
- * \param e Mesh element whose node coordinates are mapped
- * \param global_coord_system Global coordinate system
+ * \param e Mesh element whose node coordinates are mapped
+ * \param global_dim Global dimension
*/
- ElementCoordinatesMappingLocal(const Element &e, const CoordinateSystem &global_coord_system);
+ ElementCoordinatesMappingLocal(const Element &e, const unsigned global_dim);
- /// return the global coordinate system
- const CoordinateSystem getGlobalCoordinateSystem() const { return _coords; }
+ /// return the global dimension
+ unsigned getGlobalDimension() const { return _global_dim; }
/// return mapped coordinates of the node
MathLib::Point3d const& getMappedCoordinates(std::size_t node_id) const
@@ -52,7 +48,7 @@ public:
const RotationMatrix& getRotationMatrixToGlobal() const {return _matR2global;}
private:
- const CoordinateSystem _coords;
+ const unsigned _global_dim;
std::vector<MathLib::Point3d> _points;
RotationMatrix _matR2global;
};
diff --git a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
index 1938678bdc7b218eab88b7e672ab2ce0f7d906b6..7be3bcd18a60a828b87f6200bcb34eb9ae17519b 100644
--- a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
+++ b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
@@ -11,13 +11,10 @@
#include <cassert>
-#include <logog/include/logog.hpp>
+#include "BaseLib/Error.h"
#include "MeshLib/ElementCoordinatesMappingLocal.h"
-#include "MeshLib/CoordinateSystem.h"
-
#include "MeshLib/Elements/TemplateElement.h"
-
#include "MeshLib/Elements/HexRule20.h"
#include "MeshLib/Elements/HexRule8.h"
#include "MeshLib/Elements/LineRule2.h"
@@ -51,9 +48,9 @@
#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "ShapeMatrices.h"
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
namespace NumLib
{
@@ -128,10 +125,8 @@ computeMappingMatrices(
shapemat.detJ = shapemat.J.determinant();
-#ifndef NDEBUG
if (shapemat.detJ<=.0)
- ERR("det|J|=%e is not positive.\n", shapemat.detJ);
-#endif
+ OGS_FATAL("det J = %e is not positive.\n", shapemat.detJ);
}
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline
@@ -173,14 +168,14 @@ computeMappingMatrices(
computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
(ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
- if (shapemat.detJ>.0) {
+ if (shapemat.detJ > 0) {
//J^-1, dshape/dx
shapemat.invJ.noalias() = shapemat.J.inverse();
auto const nnodes(shapemat.dNdr.cols());
auto const ele_dim(shapemat.dNdr.rows());
assert(shapemat.dNdr.rows()==ele.getDimension());
- const unsigned global_dim(ele_local_coord.getGlobalCoordinateSystem().getDimension());
+ const unsigned global_dim = ele_local_coord.getGlobalDimension();
if (global_dim==ele_dim) {
shapemat.dNdx.topLeftCorner(ele_dim, nnodes).noalias() = shapemat.invJ * shapemat.dNdr;
} else {
@@ -189,8 +184,11 @@ computeMappingMatrices(
auto dshape_global = matR.topLeftCorner(3u, ele_dim) * invJ_dNdr; //3 x nnodes
shapemat.dNdx = dshape_global.topLeftCorner(global_dim, nnodes);;
}
+ } else {
+ OGS_FATAL("det J = %e is not positive.\n", shapemat.detJ);
}
}
+
template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
inline
typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
@@ -226,10 +224,10 @@ template <class T_MESH_ELEMENT,
ShapeMatrixType T_SHAPE_MATRIX_TYPE>
void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
const double* natural_pt,
- T_SHAPE_MATRICES& shapemat)
+ T_SHAPE_MATRICES& shapemat,
+ const unsigned global_dim)
{
- const MeshLib::CoordinateSystem coords(ele);
- const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele, coords);
+ const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele, global_dim);
detail::computeMappingMatrices<
T_MESH_ELEMENT,
@@ -251,7 +249,8 @@ void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
ShapeMatrixType::WHICHPART>( \
MeshLib::TemplateElement<MeshLib::RULE> const&, \
double const*, \
- SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&)
+ SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
+ const unsigned global_dim)
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(RULE, SHAPE) \
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
diff --git a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h
index b066c597d9c22761b1b3f95bf60e9041e14c14b6..debc65b6ba8558bdc9e1eab024fb5fe5e829ddf0 100644
--- a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h
+++ b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.h
@@ -28,7 +28,8 @@ template <class T_MESH_ELEMENT,
ShapeMatrixType T_SHAPE_MATRIX_TYPE>
void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
const double* natural_pt,
- T_SHAPE_MATRICES& shapemat);
+ T_SHAPE_MATRICES& shapemat,
+ const unsigned global_dim);
} // namespace detail
/**
@@ -51,12 +52,14 @@ struct NaturalCoordinatesMapping
* @param ele Mesh element object
* @param natural_pt Location in natural coordinates (r,s,t)
* @param shapemat Shape matrix data where calculated shape functions are stored
+ * @param global_dim Global dimension
*/
static void computeShapeMatrices(const T_MESH_ELEMENT& ele,
const double* natural_pt,
- T_SHAPE_MATRICES& shapemat)
+ T_SHAPE_MATRICES& shapemat,
+ const unsigned global_dim)
{
- computeShapeMatrices<ShapeMatrixType::ALL>(ele, natural_pt, shapemat);
+ computeShapeMatrices<ShapeMatrixType::ALL>(ele, natural_pt, shapemat, global_dim);
}
/**
@@ -68,17 +71,19 @@ struct NaturalCoordinatesMapping
* @param natural_pt Location in natural coordinates (r,s,t)
* @param shapemat Shape matrix data where calculated shape
* functions are stored
+ * @param global_dim Global dimension
*/
template <ShapeMatrixType T_SHAPE_MATRIX_TYPE>
static void computeShapeMatrices(const T_MESH_ELEMENT& ele,
const double* natural_pt,
- T_SHAPE_MATRICES& shapemat)
+ T_SHAPE_MATRICES& shapemat,
+ const unsigned global_dim)
{
detail::naturalCoordinatesMappingComputeShapeMatrices<
T_MESH_ELEMENT,
T_SHAPE_FUNC,
T_SHAPE_MATRICES,
- T_SHAPE_MATRIX_TYPE>(ele, natural_pt, shapemat);
+ T_SHAPE_MATRIX_TYPE>(ele, natural_pt, shapemat, global_dim);
}
};
diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index 54931172ea380465f9c92e4ab855c77440a770dc..df6de49cb2875e8765645f5cf622482f8ba346cd 100644
--- a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
+++ b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
@@ -77,16 +77,18 @@ public:
{
this->_ele = &e;
}
-
/**
* compute shape functions
*
* @param natural_pt position in natural coordinates
* @param shape evaluated shape function matrices
+ * @param global_dim global dimension
*/
- void computeShapeFunctions(const double *natural_pt, ShapeMatrices &shape) const
+ void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
+ const unsigned global_dim) const
{
- NaturalCoordsMappingType::computeShapeMatrices(*_ele, natural_pt, shape);
+ NaturalCoordsMappingType::computeShapeMatrices(*_ele, natural_pt, shape,
+ global_dim);
}
/**
@@ -95,14 +97,16 @@ public:
* @tparam T_SHAPE_MATRIX_TYPE shape matrix types to be calculated
* @param natural_pt position in natural coordinates
* @param shape evaluated shape function matrices
+ * @param global_dim global dimension
*/
template <ShapeMatrixType T_SHAPE_MATRIX_TYPE>
- void computeShapeFunctions(const double *natural_pt, ShapeMatrices &shape) const
+ void computeShapeFunctions(const double* natural_pt, ShapeMatrices& shape,
+ const unsigned global_dim) const
{
- NaturalCoordsMappingType::template computeShapeMatrices<T_SHAPE_MATRIX_TYPE>(*_ele, natural_pt, shape);
+ NaturalCoordsMappingType::template computeShapeMatrices<
+ T_SHAPE_MATRIX_TYPE>(*_ele, natural_pt, shape, global_dim);
}
-
private:
const MeshElementType* _ele;
};
diff --git a/ProcessLib/NeumannBcAssembler.h b/ProcessLib/NeumannBcAssembler.h
index 48ad268a67b7fc7883fcd2b5389708a91ce5d880..f731702f24fb2f4226d69c334fade0da4f51f9dd 100644
--- a/ProcessLib/NeumannBcAssembler.h
+++ b/ProcessLib/NeumannBcAssembler.h
@@ -66,7 +66,7 @@ public:
ShapeFunction::NPOINTS);
fe.computeShapeFunctions(
integration_method.getWeightedPoint(ip).getCoords(),
- _shape_matrices[ip]);
+ _shape_matrices[ip], GlobalDim);
}
_neumann_bc_value = value_lookup(e);
diff --git a/ProcessLib/Utils/InitShapeMatrices.h b/ProcessLib/Utils/InitShapeMatrices.h
index 874effec0107580b0e90c3d34f3af5374bc0a3d7..08e1f06d5226af793ea2ab623a3de16f7c65d907 100644
--- a/ProcessLib/Utils/InitShapeMatrices.h
+++ b/ProcessLib/Utils/InitShapeMatrices.h
@@ -40,7 +40,7 @@ initShapeMatrices(MeshLib::Element const& e, unsigned integration_order)
ShapeFunction::NPOINTS);
fe.computeShapeFunctions(
integration_method.getWeightedPoint(ip).getCoords(),
- shape_matrices[ip]);
+ shape_matrices[ip], GlobalDim);
}
return shape_matrices;
diff --git a/Tests/Data b/Tests/Data
index f93f4ce3963eac0be66358eacd5a7094f02807f9..8f42cba79bf9c67460a8fcb1604a6ca777b34f8f 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit f93f4ce3963eac0be66358eacd5a7094f02807f9
+Subproject commit 8f42cba79bf9c67460a8fcb1604a6ca777b34f8f
diff --git a/Tests/MathLib/TestIntegration.cpp b/Tests/MathLib/TestIntegration.cpp
index f4bc9dd331d513a3c1f2395d781f38ed02dd7d74..730ad3001f0ed2685f78f93c00ba92561addef22 100644
--- a/Tests/MathLib/TestIntegration.cpp
+++ b/Tests/MathLib/TestIntegration.cpp
@@ -36,11 +36,15 @@ TEST(MathLib, IntegrationGaussLegendre)
EXPECT_NEAR(2./3, MathLib::WeightedSum<MathLib::GaussLegendre<4>>::add(square),
eps);
- auto const& cube = [](double const x){ return x*x*x; };
- EXPECT_EQ(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<1>>::add(cube));
- EXPECT_EQ(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<2>>::add(cube));
- EXPECT_EQ(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<3>>::add(cube));
- EXPECT_EQ(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<4>>::add(cube));
+ auto const& cube = [](double const x) { return x * x * x; };
+ EXPECT_NEAR(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<1>>::add(cube),
+ eps);
+ EXPECT_NEAR(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<2>>::add(cube),
+ eps);
+ EXPECT_NEAR(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<3>>::add(cube),
+ eps);
+ EXPECT_NEAR(0.0, MathLib::WeightedSum<MathLib::GaussLegendre<4>>::add(cube),
+ eps);
}
diff --git a/Tests/MeshLib/TestCoordinatesMappingLocal.cpp b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
index bc028b1a5a845954ecd2acf42a9f703419d85f52..e06bf5e14f7d138186bb7840a5b2870a39d152c6 100644
--- a/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
+++ b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
@@ -19,13 +19,14 @@
#include "GeoLib/AnalyticalGeometry.h"
#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+#include "MeshLib/CoordinateSystem.h"
#include "MeshLib/Node.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Elements/Line.h"
#include "MeshLib/Elements/Quad.h"
#include "MeshLib/ElementCoordinatesMappingLocal.h"
-#include "../TestTools.h"
+#include "Tests/TestTools.h"
namespace
@@ -153,7 +154,9 @@ void debugOutput(MeshLib::Element *ele, MeshLib::ElementCoordinatesMappingLocal
TEST(MeshLib, CoordinatesMappingLocalLowerDimLineY)
{
auto ele = TestLine2::createY();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(MeshLib::CoordinateSystemType::Y));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele, MeshLib::CoordinateSystem(MeshLib::CoordinateSystemType::Y)
+ .getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
@@ -171,8 +174,10 @@ TEST(MeshLib, CoordinatesMappingLocalLowerDimLineY)
TEST(MeshLib, CoordinatesMappingLocalLowerDimLineZ)
{
auto ele = TestLine2::createZ();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele,
- MeshLib::CoordinateSystem(MeshLib::CoordinateSystemType::Z));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele,
+ MeshLib::CoordinateSystem(MeshLib::CoordinateSystemType::Z)
+ .getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
@@ -188,7 +193,8 @@ TEST(MeshLib, CoordinatesMappingLocalLowerDimLineZ)
TEST(MeshLib, CoordinatesMappingLocalLowerDimLineXY)
{
auto ele = TestLine2::createXY();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele, MeshLib::CoordinateSystem(*ele).getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
@@ -206,7 +212,8 @@ TEST(MeshLib, CoordinatesMappingLocalLowerDimLineXY)
TEST(MeshLib, CoordinatesMappingLocalLowerDimLineXYZ)
{
auto ele = TestLine2::createXYZ();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele, MeshLib::CoordinateSystem(*ele).getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
@@ -224,7 +231,8 @@ TEST(MeshLib, CoordinatesMappingLocalLowerDimLineXYZ)
TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadXZ)
{
auto ele = TestQuad4::createXZ();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele, MeshLib::CoordinateSystem(*ele).getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
@@ -244,7 +252,8 @@ TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadXZ)
TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadYZ)
{
auto ele = TestQuad4::createYZ();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele, MeshLib::CoordinateSystem(*ele).getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
@@ -264,7 +273,8 @@ TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadYZ)
TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadXYZ)
{
auto ele = TestQuad4::createXYZ();
- MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ MeshLib::ElementCoordinatesMappingLocal mapping(
+ *ele, MeshLib::CoordinateSystem(*ele).getDimension());
auto matR(mapping.getRotationMatrixToGlobal());
//debugOutput(ele, mapping);
diff --git a/Tests/NumLib/TestCoordinatesMapping.cpp b/Tests/NumLib/TestCoordinatesMapping.cpp
index 027f54fc5851d67f21a5361c4a0e214d01526653..7700333a229de03bac367f2d02a6a80485847c3c 100644
--- a/Tests/NumLib/TestCoordinatesMapping.cpp
+++ b/Tests/NumLib/TestCoordinatesMapping.cpp
@@ -62,7 +62,7 @@ public:
public:
NumLibFemNaturalCoordinatesMappingTest()
- : eps(std::numeric_limits<double>::epsilon())
+ : eps(2*std::numeric_limits<double>::epsilon())
{
// create four elements used for testing
naturalEle = this->createNaturalShape();
@@ -106,7 +106,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_N)
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
//only N
- NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::N>(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::N>(*this->naturalEle, this->r, shape, this->global_dim);
ASSERT_FALSE(shape.N.isZero());
ASSERT_TRUE(shape.dNdr.isZero());
ASSERT_TRUE(shape.J.isZero());
@@ -122,7 +122,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_DNDR)
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
// dNdr
- NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDR>(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDR>(*this->naturalEle, this->r, shape, this->global_dim);
ASSERT_TRUE(shape.N.isZero());
ASSERT_FALSE(shape.dNdr.isZero());
ASSERT_TRUE(shape.J.isZero());
@@ -139,7 +139,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_N_J)
// N_J
shape.setZero();
- NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::N_J>(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::N_J>(*this->naturalEle, this->r, shape, this->global_dim);
ASSERT_FALSE(shape.N.isZero());
ASSERT_FALSE(shape.dNdr.isZero());
ASSERT_FALSE(shape.J.isZero());
@@ -155,7 +155,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_DNDR_
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
// dNdr, J
- NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDR_J>(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDR_J>(*this->naturalEle, this->r, shape, this->global_dim);
ASSERT_TRUE(shape.N.isZero());
ASSERT_FALSE(shape.dNdr.isZero());
ASSERT_FALSE(shape.J.isZero());
@@ -172,7 +172,9 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_DNDX)
// DNDX
shape.setZero();
- NaturalCoordsMappingType::template computeShapeMatrices<ShapeMatrixType::DNDX>(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::template computeShapeMatrices<
+ ShapeMatrixType::DNDX>(*this->naturalEle, this->r, shape,
+ this->dim);
ASSERT_TRUE(shape.N.isZero());
ASSERT_FALSE(shape.dNdr.isZero());
ASSERT_FALSE(shape.J.isZero());
@@ -189,7 +191,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckFieldSpecification_ALL)
// ALL
shape.setZero();
- NaturalCoordsMappingType::computeShapeMatrices(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::computeShapeMatrices(*this->naturalEle, this->r, shape, this->global_dim);
ASSERT_FALSE(shape.N.isZero());
ASSERT_FALSE(shape.dNdr.isZero());
ASSERT_FALSE(shape.J.isZero());
@@ -205,7 +207,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckNaturalShape)
// identical to natural coordinates
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
- NaturalCoordsMappingType::computeShapeMatrices(*this->naturalEle, this->r, shape);
+ NaturalCoordsMappingType::computeShapeMatrices(*this->naturalEle, this->r, shape, this->global_dim);
double exp_J[TestFixture::dim*TestFixture::dim]= {0.0};
for (unsigned i=0; i<this->dim; i++)
exp_J[i+this->dim*i] = 1.0;
@@ -225,7 +227,7 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckIrregularShape)
// irregular shape
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
- NaturalCoordsMappingType::computeShapeMatrices(*this->irregularEle, this->r, shape);
+ NaturalCoordsMappingType::computeShapeMatrices(*this->irregularEle, this->r, shape, this->global_dim);
//std::cout << std::setprecision(16) << shape;
ASSERT_ARRAY_NEAR(this->nat_exp_N, shape.N.data(), shape.N.size(), this->eps);
@@ -243,18 +245,9 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckClockwise)
// clockwise node ordering, which is invalid)
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
- NaturalCoordsMappingType::computeShapeMatrices(*this->clockwiseEle, this->r, 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};
- 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);
- ASSERT_ARRAY_NEAR(this->cl_exp_J, shape.J.data(), shape.J.size(), this->eps);
- ASSERT_NEAR(this->cl_exp_detJ, shape.detJ, this->eps);
- ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), this->eps);
- ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), this->eps);
+ EXPECT_ANY_THROW(NaturalCoordsMappingType::computeShapeMatrices(
+ *this->clockwiseEle, this->r, shape, this->global_dim));
}
TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckZeroVolume)
@@ -263,18 +256,9 @@ TYPED_TEST(NumLibFemNaturalCoordinatesMappingTest, CheckZeroVolume)
typedef typename TestFixture::NaturalCoordsMappingType NaturalCoordsMappingType;
ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
- NaturalCoordsMappingType::computeShapeMatrices(*this->zeroVolumeEle, this->r, 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};
- 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);
- ASSERT_ARRAY_NEAR(this->ze_exp_J, shape.J.data(), shape.J.size(), this->eps);
- ASSERT_NEAR(0.0, shape.detJ, this->eps);
- ASSERT_ARRAY_NEAR(exp_invJ, shape.invJ.data(), shape.invJ.size(), this->eps);
- ASSERT_ARRAY_NEAR(exp_dNdx, shape.dNdx.data(), shape.dNdx.size(), this->eps);
+ EXPECT_ANY_THROW(NaturalCoordsMappingType::computeShapeMatrices(
+ *this->zeroVolumeEle, this->r, shape, this->global_dim));
}
TEST(NumLib, FemNaturalCoordinatesMappingLineY)
@@ -291,7 +275,7 @@ TEST(NumLib, FemNaturalCoordinatesMappingLineY)
static const unsigned dim = 1;
static const unsigned e_nnodes = 2;
ShapeMatricesType shape(dim, 2, e_nnodes);
- MappingType::computeShapeMatrices(*line, r, shape);
+ MappingType::computeShapeMatrices(*line, r, shape, 2);
double exp_J[dim*dim]= {0.0};
for (unsigned i=0; i<dim; i++)
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 0be3ad9a914db4756e7d647dc289e227c32cc40c..6b031f20d811ea2d7af66e85c115c87e09a97795 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -119,7 +119,7 @@ class NumLibFemIsoTest : public ::testing::Test, public T::TestFeType
setIdentityMatrix(dim, D);
D *= conductivity;
MeshLib::ElementCoordinatesMappingLocal ele_local_coord(
- *mesh_element, MeshLib::CoordinateSystem(*mesh_element));
+ *mesh_element, MeshLib::CoordinateSystem(*mesh_element).getDimension());
auto R = ele_local_coord.getRotationMatrixToGlobal().topLeftCorner(
TestFeType::dim, TestFeType::global_dim);
globalD.noalias() = R.transpose() * (D * R);
@@ -161,7 +161,7 @@ template <class T>
const double NumLibFemIsoTest<T>::conductivity = 1e-11;
template <class T>
-const double NumLibFemIsoTest<T>::eps = std::numeric_limits<double>::epsilon();
+const double NumLibFemIsoTest<T>::eps = 2*std::numeric_limits<double>::epsilon();
template <class T>
const unsigned NumLibFemIsoTest<T>::dim;
@@ -194,7 +194,8 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::N_J>(wp.getCoords(), shape);
+ fe.template computeShapeFunctions<ShapeMatrixType::N_J>(
+ wp.getCoords(), shape, this->global_dim);
M.noalias() += shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
}
//std::cout << "M=\n" << M;
@@ -218,7 +219,8 @@ TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
- fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(wp.getCoords(), shape);
+ fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(
+ wp.getCoords(), shape, this->global_dim);
K.noalias() += shape.dNdx.transpose() * this->globalD * shape.dNdx * shape.detJ * wp.getWeight();
}
//std::cout << "K=\n" << K << std::endl;
@@ -244,7 +246,7 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
shape.setZero();
auto wp = this->integration_method.getWeightedPoint(i);
- fe.computeShapeFunctions(wp.getCoords(), shape);
+ fe.computeShapeFunctions(wp.getCoords(), shape, this->global_dim);
M.noalias() += shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
K.noalias() += shape.dNdx.transpose() * (this->globalD * shape.dNdx) * shape.detJ * wp.getWeight();
}
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
index c2b951dbfe400f899c06127225fcce39218dfc19..7e507407f80d9ad8dd79ca19004d47acf53ef815 100644
--- a/Tests/NumLib/TestFunctionInterpolation.cpp
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -71,7 +71,7 @@ TEST(NumLibFunctionInterpolationTest, Linear1DElement)
finite_element.computeShapeFunctions(
integration_method.getWeightedPoint(0).getCoords(),
- shape_matrix);
+ shape_matrix, ShapeFunction::DIM);
ASSERT_EQ(2, shape_matrix.N.size());
// actual test