diff --git a/Documentation/Doxyfile.in b/Documentation/Doxyfile.in
index c0bf72cf2393804e78933418e9630eb97691c6f8..014864e614694ee7ddb8e26aeff38eb29f733032 100644
--- a/Documentation/Doxyfile.in
+++ b/Documentation/Doxyfile.in
@@ -692,7 +692,7 @@ LAYOUT_FILE = ${CMAKE_SOURCE_DIR}/Documentation/DoxygenLayout.xml
# LATEX_BIB_STYLE. To use this feature you need bibtex and perl available in the
# search path. See also \cite for info how to create references.
-CITE_BIB_FILES =
+CITE_BIB_FILES = ${PROJECT_SOURCE_DIR}/Documentation/bibliography
#---------------------------------------------------------------------------
# Configuration options related to warning and progress messages
diff --git a/Documentation/bibliography.bib b/Documentation/bibliography.bib
new file mode 100644
index 0000000000000000000000000000000000000000..04937bd17cafe59709369cb090be8951688419d3
--- /dev/null
+++ b/Documentation/bibliography.bib
@@ -0,0 +1,8 @@
+@book{Ericson:2004:RCD:1121584,
+ author = {Ericson, Christer},
+ title = {Real-Time Collision Detection (The Morgan Kaufmann Series in Interactive 3-D Technology) (The Morgan Kaufmann Series in Interactive 3D Technology)},
+ year = {2004},
+ isbn = {1558607323},
+ publisher = {Morgan Kaufmann Publishers Inc.},
+ address = {San Francisco, CA, USA},
+}
\ No newline at end of file
diff --git a/GeoLib/AnalyticalGeometry-impl.h b/GeoLib/AnalyticalGeometry-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..3049da7e00a84b7d2eeea43be5ac86b1686d59a1
--- /dev/null
+++ b/GeoLib/AnalyticalGeometry-impl.h
@@ -0,0 +1,142 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+namespace GeoLib
+{
+
+template <class T_POINT>
+void getNewellPlane (const std::vector<T_POINT*>& pnts,
+ MathLib::Vector3 &plane_normal,
+ double& d)
+{
+ d = 0;
+ MathLib::Vector3 centroid;
+ std::size_t n_pnts(pnts.size());
+ for (size_t i(n_pnts - 1), j(0); j < n_pnts; i = j, j++) {
+ plane_normal[0] += ((*(pnts[i]))[1] - (*(pnts[j]))[1])
+ * ((*(pnts[i]))[2] + (*(pnts[j]))[2]); // projection on yz
+ plane_normal[1] += ((*(pnts[i]))[2] - (*(pnts[j]))[2])
+ * ((*(pnts[i]))[0] + (*(pnts[j]))[0]); // projection on xz
+ plane_normal[2] += ((*(pnts[i]))[0] - (*(pnts[j]))[0])
+ * ((*(pnts[i]))[1] + (*(pnts[j]))[1]); // projection on xy
+
+ centroid += *(pnts[j]);
+ }
+
+ plane_normal.normalize();
+ d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts;
+}
+
+template<class T_MATRIX>
+void compute2DRotationMatrixToX(MathLib::Vector3 const& v,
+ T_MATRIX & rot_mat)
+{
+ const double cos_theta = v[0];
+ const double sin_theta = v[1];
+ rot_mat(0,0) = rot_mat(1,1) = cos_theta;
+ rot_mat(0,1) = sin_theta;
+ rot_mat(1,0) = -sin_theta;
+ rot_mat(0,2) = rot_mat(1,2) = rot_mat(2,0) = rot_mat(2,1) = 0.0;
+ rot_mat(2,2) = 1.0;
+}
+
+template <class T_MATRIX>
+void compute3DRotationMatrixToX(MathLib::Vector3 const& v,
+ T_MATRIX & rot_mat)
+{
+ // a vector on the plane
+ MathLib::Vector3 yy(0, 0, 0);
+ if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < std::numeric_limits<double>::epsilon()) {
+ yy[2] = 1.0;
+ } else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < std::numeric_limits<double>::epsilon()) {
+ yy[0] = 1.0;
+ } else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < std::numeric_limits<double>::epsilon()) {
+ yy[1] = 1.0;
+ } else {
+ for (unsigned i = 0; i < 3; i++) {
+ if (std::abs(v[i]) > 0.0) {
+ yy[i] = -v[i];
+ break;
+ }
+ }
+ }
+ // z"_vec
+ MathLib::Vector3 zz(MathLib::crossProduct(v, yy));
+ zz.normalize();
+ // y"_vec
+ yy = MathLib::crossProduct(zz, v);
+ yy.normalize();
+
+ for (unsigned i=0; i<3; ++i) {
+ rot_mat(0, i) = v[i];
+ rot_mat(1, i) = yy[i];
+ rot_mat(2, i) = zz[i];
+ }
+}
+
+template <class T_MATRIX>
+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;
+ rot_mat(0,1) = 0.0;
+ rot_mat(0,2) = 0.0;
+ rot_mat(1,0) = 0.0;
+ rot_mat(1,1) = 1.0;
+ rot_mat(1,2) = 0.0;
+ rot_mat(2,0) = 0.0;
+ rot_mat(2,1) = 0.0;
+ rot_mat(2,2) = 1.0;
+ return;
+ }
+
+ // sqrt (n_1^2 + n_3^2)
+ double const h0(sqrt(n[0]*n[0]+n[2]*n[2]));
+
+ // In case the x and z components of the normal are both zero the rotation
+ // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
+ // required. The angle is either pi/2 or 3/2*pi. Thus the components of
+ // rot_mat are as follows.
+ if (h0 < std::numeric_limits<double>::epsilon()) {
+ rot_mat(0,0) = 1.0;
+ rot_mat(0,1) = 0.0;
+ rot_mat(0,2) = 0.0;
+ rot_mat(1,0) = 0.0;
+ rot_mat(1,1) = 0.0;
+ if (n[1] > 0)
+ rot_mat(1,2) = -1.0;
+ else
+ rot_mat(1,2) = 1.0;
+ rot_mat(2,0) = 0.0;
+ if (n[1] > 0)
+ rot_mat(2,1) = 1.0;
+ else
+ rot_mat(2,1) = -1.0;
+ rot_mat(2,2) = 0.0;
+ return;
+ }
+
+ double h1(1 / n.getLength());
+
+ // general case: calculate entries of rotation matrix
+ rot_mat(0, 0) = n[2] / h0;
+ rot_mat(0, 1) = 0;
+ rot_mat(0, 2) = - n[0] / h0;
+ rot_mat(1, 0) = - n[1]*n[0]/h0 * h1;
+ rot_mat(1, 1) = h0 * h1;
+ rot_mat(1, 2) = - n[1]*n[2]/h0 * h1;
+ rot_mat(2, 0) = n[0] * h1;
+ rot_mat(2, 1) = n[1] * h1;
+ rot_mat(2, 2) = n[2] * h1;
+}
+
+} // end namespace GeoLib
+
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index d651df8b41aaab8f78dcbfabee6853333f3ac36a..67b5f1d6e0f4f1c3e246110d56efd4fe258203dd 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -348,84 +348,6 @@ double calcTetrahedronVolume(MathLib::Point3d const& x1,
return std::abs(GeoLib::scalarTriple(ac, ad, ab)) / 6.0;
}
-// NewellPlane from book Real-Time Collision detection p. 494
-void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, MathLib::Vector3 &plane_normal, double& d)
-{
- d = 0;
- MathLib::Vector3 centroid;
- size_t n_pnts(pnts.size());
- for (size_t i(n_pnts - 1), j(0); j < n_pnts; i = j, j++) {
- plane_normal[0] += ((*(pnts[i]))[1] - (*(pnts[j]))[1])
- * ((*(pnts[i]))[2] + (*(pnts[j]))[2]); // projection on yz
- plane_normal[1] += ((*(pnts[i]))[2] - (*(pnts[j]))[2])
- * ((*(pnts[i]))[0] + (*(pnts[j]))[0]); // projection on xz
- plane_normal[2] += ((*(pnts[i]))[0] - (*(pnts[j]))[0])
- * ((*(pnts[i]))[1] + (*(pnts[j]))[1]); // projection on xy
-
- centroid += *(pnts[j]);
- }
-
- plane_normal *= 1.0 / plane_normal.getLength();
- d = MathLib::scalarProduct(centroid, plane_normal) / n_pnts;
-}
-
-void computeRotationMatrixToXY(MathLib::Vector3 const& n,
- MathLib::DenseMatrix<double> & 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;
- rot_mat(0,1) = 0.0;
- rot_mat(0,2) = 0.0;
- rot_mat(1,0) = 0.0;
- rot_mat(1,1) = 1.0;
- rot_mat(1,2) = 0.0;
- rot_mat(2,0) = 0.0;
- rot_mat(2,1) = 0.0;
- rot_mat(2,2) = 1.0;
- return;
- }
-
- // sqrt (n_1^2 + n_3^2)
- double const h0(sqrt(n[0]*n[0]+n[2]*n[2]));
-
- // In case the x and z components of the normal are both zero the rotation
- // to the x-z-plane is not required, i.e. only the rotation in the z-axis is
- // required. The angle is either pi/2 or 3/2*pi. Thus the components of
- // rot_mat are as follows.
- if (h0 < std::numeric_limits<double>::epsilon()) {
- rot_mat(0,0) = 1.0;
- rot_mat(0,1) = 0.0;
- rot_mat(0,2) = 0.0;
- rot_mat(1,0) = 0.0;
- rot_mat(1,1) = 0.0;
- if (n[1] > 0)
- rot_mat(1,2) = -1.0;
- else
- rot_mat(1,2) = 1.0;
- rot_mat(2,0) = 0.0;
- if (n[1] > 0)
- rot_mat(2,1) = 1.0;
- else
- rot_mat(2,1) = -1.0;
- rot_mat(2,2) = 0.0;
- return;
- }
-
- double h1(1 / n.getLength());
-
- // general case: calculate entries of rotation matrix
- rot_mat(0, 0) = n[2] / h0;
- rot_mat(0, 1) = 0;
- rot_mat(0, 2) = - n[0] / h0;
- rot_mat(1, 0) = - n[1]*n[0]/h0 * h1;
- rot_mat(1, 1) = h0 * h1;
- rot_mat(1, 2) = - n[1]*n[2]/h0 * h1;
- rot_mat(2, 0) = n[0] * h1;
- rot_mat(2, 1) = n[1] * h1;
- rot_mat(2, 2) = n[2] * h1;
-}
-
void computeRotationMatrixToXZ(MathLib::Vector3 const& plane_normal, MathLib::DenseMatrix<double> & rot_mat)
{
// *** some frequently used terms ***
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 1198202138b4b9b7464f52ccb4a3f7814a7f75a8..e07a73b4609a326ac711f3add921e24823d5ef34 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -55,22 +55,41 @@ Orientation getOrientation (const GeoLib::Point* p0,
/**
* compute a supporting plane (represented by plane_normal and the value d) for the polygon
* Let \f$n\f$ be the plane normal and \f$d\f$ a parameter. Then for all points \f$p \in R^3\f$ of the plane
- * it holds \f$ n \cdot p + d = 0\f$
+ * it holds \f$ n \cdot p + d = 0\f$. The Newell algorithm is described in
+ * \cite Ericson:2004:RCD:1121584 .
* @param pnts points of a closed polyline describing a polygon
* @param plane_normal the normal of the plane the polygon is located in
* @param d parameter from the plane equation
*/
-void getNewellPlane (const std::vector<GeoLib::Point*>& pnts,
+template <class T_POINT>
+void getNewellPlane (const std::vector<T_POINT*>& pnts,
MathLib::Vector3 &plane_normal,
double& d);
+/**
+ * Computes a rotation matrix that rotates the given 2D normal vector parallel to X-axis
+ * @param v a 2D normal vector to be rotated
+ * @param rot_mat a 2x2 rotation matrix
+ */
+template<class T_MATRIX>
+void compute2DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat);
+
+/**
+ * Computes a rotation matrix that rotates the given 3D normal vector parallel to X-axis
+ * @param v a 3D normal vector to be rotated
+ * @param rot_mat a 3x3 rotation matrix
+ */
+template <class T_MATRIX>
+void compute3DRotationMatrixToX(MathLib::Vector3 const& v, T_MATRIX & rot_mat);
+
/**
* Method computes the rotation matrix that rotates the given vector parallel to the \f$z\f$ axis.
* @param plane_normal the (3d) vector that is rotated parallel to the \f$z\f$ axis
* @param rot_mat 3x3 rotation matrix
*/
+template <class T_MATRIX>
void computeRotationMatrixToXY(MathLib::Vector3 const& plane_normal,
- MathLib::DenseMatrix<double> & rot_mat);
+ T_MATRIX & rot_mat);
/**
* Method computes the rotation matrix that rotates the given vector parallel to the \f$y\f$ axis.
@@ -296,4 +315,6 @@ GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
} // end namespace GeoLib
+#include "AnalyticalGeometry-impl.h"
+
#endif /* ANALYTICAL_GEOMETRY_H_ */
diff --git a/MathLib/LinAlg/Dense/DenseMatrix-impl.h b/MathLib/LinAlg/Dense/DenseMatrix-impl.h
index 131ab32940df82cba0e4d28f2a502bbd71851cfe..6abda2755fb346fe6ef861c4d06f7401c343f35c 100644
--- a/MathLib/LinAlg/Dense/DenseMatrix-impl.h
+++ b/MathLib/LinAlg/Dense/DenseMatrix-impl.h
@@ -219,6 +219,26 @@ DenseMatrix<FP_TYPE, IDX_TYPE>::transpose() const
return y;
}
+template<typename FP_TYPE, typename IDX_TYPE>
+void
+DenseMatrix<FP_TYPE, IDX_TYPE>::transposeInPlace()
+{
+ if (_n_rows==_n_cols) { // square matrix
+ for (IDX_TYPE i = 0; i < _n_rows; i++)
+ for (IDX_TYPE j = i+1; j < _n_cols; j++)
+ std::swap(_data[address(i, j)], _data[address(j, i)]);
+ } else { // non-square matrix
+ const DenseMatrix<FP_TYPE, IDX_TYPE> org(*this);
+ std::swap(_n_rows, _n_cols);
+ for (IDX_TYPE i = 0; i < _n_rows; i++) {
+ for (IDX_TYPE j = 0; j < _n_cols; j++) {
+ _data[address(i, j)] = org(j, i);
+ }
+ }
+ }
+
+}
+
template<typename FP_TYPE, typename IDX_TYPE>
DenseMatrix<FP_TYPE, IDX_TYPE>*
DenseMatrix<FP_TYPE, IDX_TYPE>::getSubMatrix(
@@ -285,6 +305,16 @@ DenseMatrix<FP_TYPE, IDX_TYPE>::write (std::ostream &out) const
}
}
+template <typename FP_TYPE, typename IDX_TYPE>
+void
+DenseMatrix<FP_TYPE, IDX_TYPE>::setIdentity()
+{
+ (*this) = 0.0;
+ const IDX_TYPE n_square_rows = std::min(_n_rows, _n_cols);
+ for (IDX_TYPE i=0; i<n_square_rows; i++)
+ _data[address(i,i)] = 1.0;
+}
+
template <typename FP_TYPE, typename IDX_TYPE>
FP_TYPE
sqrFrobNrm (const DenseMatrix<FP_TYPE, IDX_TYPE> &mat)
diff --git a/MathLib/LinAlg/Dense/DenseMatrix.h b/MathLib/LinAlg/Dense/DenseMatrix.h
index 6dc6c659a4b7cd138d537c4fd8cfb0048317cc30..a26256876d54370243d873759dbd78dc25d47cf2 100644
--- a/MathLib/LinAlg/Dense/DenseMatrix.h
+++ b/MathLib/LinAlg/Dense/DenseMatrix.h
@@ -57,6 +57,14 @@ public:
return _data;
}
+ /**
+ * Get all entries, which are stored in an array.
+ */
+ const FP_TYPE *data() const
+ {
+ return _data;
+ }
+
/**
* Assignment operator, makes a copy of the internal data of the object.
* @param rhs The DenseMatrix object to the right side of the assignment symbol.
@@ -109,6 +117,11 @@ public:
*/
DenseMatrix* transpose() const; // HB & ZC
+ /**
+ * transpose the matrix in place
+ */
+ void transposeInPlace();
+
DenseMatrix* getSubMatrix (IDX_TYPE b_row, IDX_TYPE b_col, IDX_TYPE e_row, IDX_TYPE e_col) const;
/**
@@ -139,6 +152,16 @@ public:
*/
IDX_TYPE getNCols () const { return _n_cols; }
+ /**
+ * get the number of entries in the matrix
+ */
+ IDX_TYPE size() const { return _n_rows*_n_cols; }
+
+ /**
+ * set the identity matrix
+ */
+ void setIdentity();
+
protected:
/**
* the number of rows
diff --git a/MathLib/Point3d.h b/MathLib/Point3d.h
index 75368621aa822dd8c9df062d82378b795938e4e7..d59fc2ab998c537e6808514599d15bb5c3736945 100644
--- a/MathLib/Point3d.h
+++ b/MathLib/Point3d.h
@@ -42,5 +42,24 @@ bool lessEq(const MathLib::Point3d& p0,
const MathLib::Point3d& p1,
double tol = std::numeric_limits<double>::epsilon());
+/**
+ * rotation of points
+ * @param mat a rotation matrix
+ * @param p a point to be transformed
+ * @return a rotated point
+ */
+template <typename MATRIX>
+inline MathLib::Point3d operator*(MATRIX const& mat, MathLib::Point3d const& p)
+{
+ MathLib::Point3d new_p;
+ for (std::size_t i(0); i<3; ++i) {
+ new_p[i] = 0.0;
+ for (std::size_t j(0); j<3; ++j) {
+ new_p[i] += mat(i,j)*p[j];
+ }
+ }
+ return new_p;
+}
+
#endif /* POINT3D_H_ */
diff --git a/MathLib/TemplatePoint.h b/MathLib/TemplatePoint.h
index d6e3f5a4c1cbade978b8eb455a8d82e21ec3788b..1a58d4f2465bfc033ab639e8894ee71c699238de 100644
--- a/MathLib/TemplatePoint.h
+++ b/MathLib/TemplatePoint.h
@@ -78,6 +78,13 @@ public:
return _x.data();
}
+ /** set the coordinates of the point */
+ void setCoords (const T* coords)
+ {
+ for (unsigned i=0; i<DIM; i++)
+ _x[i] = coords[i];
+ }
+
/** write point coordinates into stream (used from operator<<)
* \param os a standard output stream
*/
diff --git a/MeshLib/CoordinateSystem.cpp b/MeshLib/CoordinateSystem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..03209a90fed6ec1eef254c322644c2025f934fcb
--- /dev/null
+++ b/MeshLib/CoordinateSystem.cpp
@@ -0,0 +1,34 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "CoordinateSystem.h"
+
+#include "MeshLib/Node.h"
+#include "MeshLib/Elements/Element.h"
+
+namespace MeshLib
+{
+
+CoordinateSystem::CoordinateSystem(const Element &ele)
+{
+ const GeoLib::AABB<MeshLib::Node> aabb(ele.getNodes(),
+ ele.getNodes()+ele.getNNodes());
+ const CoordinateSystem bboxCoordSys(getCoordinateSystem(aabb));
+ if (bboxCoordSys.getDimension() >= ele.getDimension()) {
+ _type = bboxCoordSys.getType();
+ } else { // e.g. zero volume elements
+ if (ele.getDimension()>=1)
+ _type = CoordinateSystemType::X;
+ if (ele.getDimension()>=2)
+ _type |= CoordinateSystemType::Y;
+ if (ele.getDimension()==3)
+ _type |= CoordinateSystemType::Z;
+ }
+}
+
+} // end
diff --git a/MeshLib/CoordinateSystem.h b/MeshLib/CoordinateSystem.h
new file mode 100644
index 0000000000000000000000000000000000000000..d8d7af4747d43b3b09a832ba0fefc54a53dba641
--- /dev/null
+++ b/MeshLib/CoordinateSystem.h
@@ -0,0 +1,100 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 COORDINATESYSTEMTYPE_H_
+#define COORDINATESYSTEMTYPE_H_
+
+#include <cmath>
+
+#include "GeoLib/AABB.h"
+
+namespace MeshLib
+{
+
+class Element;
+
+/**
+ * \brief Coordinate system type
+ */
+struct CoordinateSystemType
+{
+ enum type
+ {
+ X = 0x01,
+ Y = 0x02,
+ Z = 0x04
+ };
+};
+
+/**
+ * \brief Coordinate systems
+ *
+ *
+ */
+class CoordinateSystem
+{
+public:
+ /// User provided coordinate system
+ explicit CoordinateSystem(unsigned char coord) : _type (coord) {}
+
+ /// Decides for a given element
+ explicit CoordinateSystem(const Element &ele);
+
+ /// Decides a coordinate system from a bounding box
+ template <class T>
+ explicit CoordinateSystem(const GeoLib::AABB<T> &bbox) : _type(getCoordinateSystem(bbox)) {}
+
+ /// get this coordinate type
+ unsigned char getType() const { return _type; }
+
+ /// get dimension size
+ unsigned getDimension() const {
+ if (hasZ())
+ return 3;
+ else if (hasY())
+ return 2;
+ else
+ return 1;
+ }
+
+ /// has X dimension
+ bool hasX() const { return (_type & CoordinateSystemType::type::X); }
+
+ /// has Y dimension
+ bool hasY() const { return (_type & CoordinateSystemType::type::Y); }
+
+ /// has z dimension
+ bool hasZ() const { return (_type & CoordinateSystemType::type::Z); }
+
+private:
+ template <class T>
+ unsigned char getCoordinateSystem(const GeoLib::AABB<T> &bbox) const;
+
+ unsigned char _type;
+};
+
+template <class T>
+unsigned char CoordinateSystem::getCoordinateSystem(const GeoLib::AABB<T> &bbox) const
+{
+ unsigned char coords = 0;
+
+ const T pt_diff = bbox.getMaxPoint() - bbox.getMinPoint();
+
+ if (std::abs(pt_diff[0]) > .0)
+ coords |= CoordinateSystemType::X;
+ if (std::abs(pt_diff[1]) > .0)
+ coords |= CoordinateSystemType::Y;
+ if (std::abs(pt_diff[2]) > .0)
+ coords |= CoordinateSystemType::Z;
+
+ return coords;
+}
+
+} // MeshLib
+
+#endif // COORDINATESYSTEMTYPE_H_
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ccaf5e57ad619b5ad93326bd3fd1e613d8ed716
--- /dev/null
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -0,0 +1,86 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "ElementCoordinatesMappingLocal.h"
+
+#include <limits>
+#include <cassert>
+
+#include "GeoLib/AnalyticalGeometry.h"
+
+#include "MathLib/MathTools.h"
+#include "MeshLib/Node.h"
+
+namespace MeshLib
+{
+
+ElementCoordinatesMappingLocal::ElementCoordinatesMappingLocal(
+ const Element& e,
+ const CoordinateSystem &global_coords)
+: _coords(global_coords), _matR2global(3,3)
+{
+ assert(e.getDimension() <= global_coords.getDimension());
+ for(unsigned i = 0; i < e.getNNodes(); i++)
+ _vec_nodes.push_back(new MeshLib::Node(*(e.getNode(i))));
+
+ getRotationMatrixToGlobal(e, global_coords, _vec_nodes, _matR2global);
+#ifdef OGS_USE_EIGEN
+ rotateToLocal(_matR2global.transpose(), _vec_nodes);
+#else
+ RotationMatrix* m(_matR2global.transpose());
+ rotateToLocal(*m, _vec_nodes);
+ delete m;
+#endif
+}
+
+ElementCoordinatesMappingLocal::~ElementCoordinatesMappingLocal()
+{
+ for (auto p : _vec_nodes) delete p;
+}
+
+void ElementCoordinatesMappingLocal::rotateToLocal(
+ const RotationMatrix &matR2local,
+ std::vector<MeshLib::Node*> &vec_nodes) const
+{
+ for (MeshLib::Node* node : vec_nodes)
+ node->setCoords((matR2local*static_cast<MathLib::Point3d>(*node)).getCoords());
+}
+
+void ElementCoordinatesMappingLocal::getRotationMatrixToGlobal(
+ const Element &e,
+ const CoordinateSystem &global_coords,
+ const std::vector<MeshLib::Node*> &vec_nodes,
+ RotationMatrix &matR) const
+{
+ const std::size_t global_dim = global_coords.getDimension();
+
+ // compute R in x=R*x' where x are original coordinates and x' are local coordinates
+ if (global_dim == e.getDimension()) {
+ matR.setIdentity();
+ } else if (e.getDimension() == 1) {
+ MathLib::Vector3 xx(*vec_nodes[0], *vec_nodes[1]);
+ xx.normalize();
+ if (global_dim == 2)
+ GeoLib::compute2DRotationMatrixToX(xx, matR);
+ else
+ GeoLib::compute3DRotationMatrixToX(xx, matR);
+ matR.transposeInPlace();
+ } else if (global_dim == 3 && e.getDimension() == 2) {
+ // get plane normal
+ MathLib::Vector3 plane_normal;
+ double d;
+ GeoLib::getNewellPlane (vec_nodes, plane_normal, d);
+ // compute a rotation matrix to XY
+ GeoLib::computeRotationMatrixToXY(plane_normal, matR);
+ // set a transposed matrix
+ matR.transposeInPlace();
+ }
+
+}
+
+} // MeshLib
diff --git a/MeshLib/ElementCoordinatesMappingLocal.h b/MeshLib/ElementCoordinatesMappingLocal.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a46788363ec78c2c5df40eb57b3a3600e4984d2
--- /dev/null
+++ b/MeshLib/ElementCoordinatesMappingLocal.h
@@ -0,0 +1,79 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, 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 ELEMENTCOORDINATESMAPPINGLOCAL_H_
+#define ELEMENTCOORDINATESMAPPINGLOCAL_H_
+
+#include <vector>
+
+#ifdef OGS_USE_EIGEN
+#include <Eigen/Eigen>
+#endif
+
+#include "MathLib/Vector3.h"
+#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/CoordinateSystem.h"
+
+namespace MeshLib
+{
+
+#ifdef OGS_USE_EIGEN
+typedef Eigen::Matrix<double, 3u, 3u, Eigen::RowMajor> RotationMatrix;
+#else
+typedef MathLib::DenseMatrix<double> RotationMatrix;
+#endif
+
+/**
+ * This class maps node coordinates on intrinsic coordinates of the given element.
+ */
+class ElementCoordinatesMappingLocal
+{
+public:
+ /**
+ * Constructor
+ * \param e Mesh element whose node coordinates are mapped
+ * \param global_coord_system Global coordinate system
+ */
+ ElementCoordinatesMappingLocal(const Element &e, const CoordinateSystem &global_coord_system);
+
+ /// Destructor
+ virtual ~ElementCoordinatesMappingLocal();
+
+ /// return the global coordinate system
+ const CoordinateSystem getGlobalCoordinateSystem() const { return _coords; }
+
+ /// return mapped coordinates of the node
+ const MeshLib::Node* getMappedCoordinates(size_t node_id) const
+ {
+ return _vec_nodes[node_id];
+ }
+
+ /// return a rotation matrix converting to global coordinates
+ const RotationMatrix& getRotationMatrixToGlobal() const {return _matR2global;}
+
+private:
+ /// rotate points to local coordinates
+ void rotateToLocal(const RotationMatrix &matR2local, std::vector<MeshLib::Node*> &vec_pt) const;
+
+ /// get a rotation matrix to the global coordinates
+ /// it computes R in x=R*x' where x is original coordinates and x' is local coordinates
+ void getRotationMatrixToGlobal(
+ const Element &e, const CoordinateSystem &coordinate_system,
+ const std::vector<MeshLib::Node*> &vec_pt, RotationMatrix &matR2original) const;
+
+private:
+ const CoordinateSystem _coords;
+ std::vector<MeshLib::Node*> _vec_nodes;
+ RotationMatrix _matR2global;
+};
+
+}
+
+#endif // ELEMENTCOORDINATESMAPPINGLOCAL_H_
diff --git a/Tests/MathLib/TestDenseMatrix.cpp b/Tests/MathLib/TestDenseMatrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..41ad1fe1df845816808f7d489698c3d164fe420c
--- /dev/null
+++ b/Tests/MathLib/TestDenseMatrix.cpp
@@ -0,0 +1,67 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "gtest/gtest.h"
+
+#include "MathLib/LinAlg/Dense/DenseMatrix.h"
+#include "../TestTools.h"
+
+using namespace MathLib;
+
+TEST(MathLib, DenseMatrixTransposeInPlace)
+{
+ const double eps(std::numeric_limits<double>::epsilon());
+
+ // square matrix
+ DenseMatrix<double> m1(3,3);
+ unsigned cnt = 0;
+ for (unsigned i=0; i<m1.getNRows(); i++)
+ for (unsigned j=0; j<m1.getNCols(); j++)
+ m1(i,j) = cnt++;
+ double expected_m1[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+ ASSERT_EQ(3u, m1.getNRows());
+ ASSERT_EQ(3u, m1.getNCols());
+ ASSERT_ARRAY_NEAR(expected_m1, m1.data(), m1.size(), eps);
+ m1.transposeInPlace();
+ ASSERT_EQ(3u, m1.getNRows());
+ ASSERT_EQ(3u, m1.getNCols());
+ double expected_m1t[] = {0, 3, 6, 1, 4, 7, 2, 5, 8};
+ ASSERT_ARRAY_NEAR(expected_m1t, m1.data(), m1.size(), eps);
+
+ // non-square matrix 1
+ DenseMatrix<double> m2(2,3);
+ cnt = 0;
+ for (unsigned i=0; i<m2.getNRows(); i++)
+ for (unsigned j=0; j<m2.getNCols(); j++)
+ m2(i,j) = cnt++;
+ ASSERT_EQ(2u, m2.getNRows());
+ ASSERT_EQ(3u, m2.getNCols());
+ double expected_m2[] = {0, 1, 2, 3, 4, 5};
+ ASSERT_ARRAY_NEAR(expected_m2, m2.data(), m2.size(), eps);
+ m2.transposeInPlace();
+ ASSERT_EQ(3u, m2.getNRows());
+ ASSERT_EQ(2u, m2.getNCols());
+ double expected_m2t[] = {0, 3, 1, 4, 2, 5};
+ ASSERT_ARRAY_NEAR(expected_m2t, m2.data(), m2.size(), eps);
+
+ // non-square matrix 2
+ DenseMatrix<double> m3(3,2);
+ cnt = 0;
+ for (unsigned i=0; i<m3.getNRows(); i++)
+ for (unsigned j=0; j<m3.getNCols(); j++)
+ m3(i,j) = cnt++;
+ ASSERT_EQ(3u, m3.getNRows());
+ ASSERT_EQ(2u, m3.getNCols());
+ double expected_m3[] = {0, 1, 2, 3, 4, 5};
+ ASSERT_ARRAY_NEAR(expected_m3, m3.data(), m3.size(), eps);
+ m3.transposeInPlace();
+ ASSERT_EQ(2u, m3.getNRows());
+ ASSERT_EQ(3u, m3.getNCols());
+ double expected_m3t[] = {0, 2, 4, 1, 3, 5};
+ ASSERT_ARRAY_NEAR(expected_m3t, m3.data(), m3.size(), eps);
+}
diff --git a/Tests/MeshLib/TestCoordinatesMappingLocal.cpp b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e26792226911bb7276b31f1aff54d1c3f0c077c3
--- /dev/null
+++ b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
@@ -0,0 +1,260 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <gtest/gtest.h>
+
+#include <limits>
+#include <algorithm>
+#include <vector>
+#include <cmath>
+
+#ifdef OGS_USE_EIGEN
+#include <Eigen/Eigen>
+#endif
+
+#include "GeoLib/AnalyticalGeometry.h"
+#include "MathLib/LinAlg/Dense/DenseMatrix.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"
+
+
+namespace
+{
+
+class TestLine2
+{
+public:
+ typedef MeshLib::Line ElementType;
+ static const unsigned e_nnodes = ElementType::n_all_nodes;
+
+ static MeshLib::Line* createY()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, -1.0, 0.0);
+ nodes[1] = new MeshLib::Node(0.0, 1.0, 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ static MeshLib::Line* createZ()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, -1.0);
+ nodes[1] = new MeshLib::Node(0.0, 0.0, 1.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ static MeshLib::Line* createXY()
+ {
+ // 45degree inclined
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(2./sqrt(2), 2./sqrt(2), 0.0);
+ return new MeshLib::Line(nodes);
+ }
+
+ static MeshLib::Line* createXYZ()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0);
+ nodes[1] = new MeshLib::Node(2./sqrt(3), 2./sqrt(3), 2./sqrt(3));
+ return new MeshLib::Line(nodes);
+ }
+
+};
+
+class TestQuad4
+{
+ public:
+ // Element information
+ typedef MeshLib::Quad ElementType;
+ static const unsigned e_nnodes = ElementType::n_all_nodes;
+
+ // 2.5D case: inclined
+ static MeshLib::Quad* createXYZ()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ // rotate 45 degree around x axis
+ nodes[0] = new MeshLib::Node( 1.0, 0.7071067811865475, 0.7071067811865475);
+ nodes[1] = new MeshLib::Node(-1.0, 0.7071067811865475, 0.7071067811865475);
+ nodes[2] = new MeshLib::Node(-1.0, -0.7071067811865475, -0.7071067811865475);
+ nodes[3] = new MeshLib::Node( 1.0, -0.7071067811865475, -0.7071067811865475);
+ return new MeshLib::Quad(nodes);
+ }
+
+ // 2.5D case: inclined
+ static MeshLib::Quad* createXZ()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node( 1.0, 0.0, 1.0);
+ nodes[1] = new MeshLib::Node(-1.0, 0.0, 1.0);
+ nodes[2] = new MeshLib::Node(-1.0, 0.0, -1.0);
+ nodes[3] = new MeshLib::Node( 1.0, 0.0, -1.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+ // 2.5D case: inclined
+ static MeshLib::Quad* createYZ()
+ {
+ MeshLib::Node** nodes = new MeshLib::Node*[e_nnodes];
+ nodes[0] = new MeshLib::Node(0.0, 1.0, 1.0);
+ nodes[1] = new MeshLib::Node(0.0, -1.0, 1.0);
+ nodes[2] = new MeshLib::Node(0.0, -1.0, -1.0);
+ nodes[3] = new MeshLib::Node(0.0, 1.0, -1.0);
+ return new MeshLib::Quad(nodes);
+ }
+
+};
+
+#if 0
+// keep this function for debugging
+void debugOutput(MeshLib::Element *ele, MeshLib::ElementCoordinatesMappingLocal &mapping)
+{
+ std::cout.precision(12);
+ std::cout << "original" << std::endl;
+ for (unsigned i=0; i<ele->getNNodes(); i++)
+ std::cout << *ele->getNode(i) << std::endl;
+ std::cout << "local coords=" << std::endl;
+ for (unsigned i=0; i<ele->getNNodes(); i++)
+ std::cout << *mapping.getMappedCoordinates(i) << std::endl;
+ std::cout << "R=\n" << mapping.getRotationMatrixToGlobal() << std::endl;
+ auto matR(mapping.getRotationMatrixToGlobal());
+ std::cout << "global coords=" << std::endl;
+ for (unsigned i=0; i<ele->getNNodes(); i++) {
+ double* raw = const_cast<double*>(&(*mapping.getMappedCoordinates(i))[0]);
+ Eigen::Map<Eigen::Vector3d> v(raw);
+ std::cout << (matR*v).transpose() << std::endl;
+ }
+}
+#endif
+
+// check if using the rotation matrix results in the original coordinates
+#define CHECK_COORDS(ele, mapping)\
+ for (unsigned ii=0; ii<ele->getNNodes(); ii++) {\
+ MathLib::Point3d global(matR*static_cast<MathLib::Point3d>(*mapping.getMappedCoordinates(ii)));\
+ const double eps(std::numeric_limits<double>::epsilon());\
+ ASSERT_ARRAY_NEAR(&(*ele->getNode(ii))[0], global.getCoords(), 3u, eps);\
+ }
+
+} //namespace
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimLineY)
+{
+ auto ele = TestLine2::createY();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(MeshLib::CoordinateSystemType::Y));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ double exp_R[3*3] = {0, -1, 0,
+ 1, 0, 0,
+ 0, 0, 1};
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimLineZ)
+{
+ auto ele = TestLine2::createZ();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele,
+ MeshLib::CoordinateSystem(MeshLib::CoordinateSystemType::Z));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ double exp_R[3*3] = {0, 0, -1, 0, 1, 0, 1, 0, 0};
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimLineXY)
+{
+ auto ele = TestLine2::createXY();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ double exp_R[3*3] = {0.70710678118654757, -0.70710678118654757, 0,
+ 0.70710678118654757, 0.70710678118654757, 0,
+ 0, 0, 1};
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimLineXYZ)
+{
+ auto ele = TestLine2::createXYZ();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ double exp_R[3*3] = {0.57735026918962584, -0.81649658092772626, 0,
+ 0.57735026918962584, 0.40824829046386313, -0.70710678118654757,
+ 0.57735026918962584, 0.40824829046386313, 0.70710678118654757};
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadXZ)
+{
+ auto ele = TestQuad4::createXZ();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ // results when using GeoLib::ComputeRotationMatrixToXY()
+ double exp_R[3*3] = { 1, 0, 0,
+ 0, 0, -1,
+ 0, 1, 0};
+
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadYZ)
+{
+ auto ele = TestQuad4::createYZ();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ // results when using GeoLib::ComputeRotationMatrixToXY()
+ double exp_R[3*3] = { 0, 0, 1,
+ 0, 1, 0,
+ -1, 0, 0};
+
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
+TEST(MeshLib, CoordinatesMappingLocalLowerDimQuadXYZ)
+{
+ auto ele = TestQuad4::createXYZ();
+ MeshLib::ElementCoordinatesMappingLocal mapping(*ele, MeshLib::CoordinateSystem(*ele));
+ auto matR(mapping.getRotationMatrixToGlobal());
+ //debugOutput(ele, mapping);
+
+ // results when using GeoLib::ComputeRotationMatrixToXY()
+ double exp_R[3*3] = { 1, 0, 0,
+ 0, 0.70710678118654757, -0.70710678118654757,
+ 0, 0.70710678118654757, 0.70710678118654757};
+
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 24ecfde7d9f863fbc773b5f52471efa79bdf9de1..834055a367cad392192ca7131fbb8b422cb91c3b 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -140,6 +140,10 @@ class NumLibFemIsoTest : public ::testing::Test, public T::TestFeType
std::vector<const MeshElementType*> vec_eles;
MeshElementType* mesh_element;
+ public:
+#ifdef OGS_USE_EIGEN
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW // required to use fixed size Eigen matrices
+#endif
}; // NumLibFemIsoTest
template <class T>
diff --git a/scripts/cmake/cmake/FindEigen3.cmake b/scripts/cmake/cmake/FindEigen3.cmake
index 2b9e6afae4dcc348c53d96b052efba5a3b161828..0f8c5d45c79ce9a97b224a4ae2f2949c3a58b30a 100644
--- a/scripts/cmake/cmake/FindEigen3.cmake
+++ b/scripts/cmake/cmake/FindEigen3.cmake
@@ -85,3 +85,6 @@ else ()
endif()
+if(EIGEN3_INCLUDE_DIR)
+ message(STATUS "Eigen version ${EIGEN3_VERSION}")
+endif()