From c4d4ad6726ee6c60fa410ab33c0971451ddfcfde Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Thu, 23 Apr 2015 12:57:28 +0200
Subject: [PATCH] add ElementCoordinatesMappingLocal class
---
MeshLib/ElementCoordinatesMappingLocal.cpp | 97 +++++++
MeshLib/ElementCoordinatesMappingLocal.h | 79 +++++
Tests/MeshLib/TestCoordinatesMappingLocal.cpp | 270 ++++++++++++++++++
3 files changed, 446 insertions(+)
create mode 100644 MeshLib/ElementCoordinatesMappingLocal.cpp
create mode 100644 MeshLib/ElementCoordinatesMappingLocal.h
create mode 100644 Tests/MeshLib/TestCoordinatesMappingLocal.cpp
diff --git a/MeshLib/ElementCoordinatesMappingLocal.cpp b/MeshLib/ElementCoordinatesMappingLocal.cpp
new file mode 100644
index 00000000000..c6a5c513d76
--- /dev/null
+++ b/MeshLib/ElementCoordinatesMappingLocal.cpp
@@ -0,0 +1,97 @@
+/**
+ * \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)
+{
+ assert(e.getDimension() <= global_coords.getDimension());
+ for(size_t i = 0; i < e.getNNodes(); i++)
+ _point_vec.push_back(MeshLib::Node(*(e.getNode(i))));
+
+ getRotationMatrixToGlobal(e, global_coords, _point_vec, _matR2global);
+ rotateToLocal(e, global_coords, _point_vec, _matR2global.transpose(), _point_vec);
+}
+
+void ElementCoordinatesMappingLocal::rotateToLocal(
+ const Element &ele,
+ const CoordinateSystem &global_coords,
+ const std::vector<MeshLib::Node> &vec_pt,
+ const RotationMatrix &matR2local,
+ std::vector<MeshLib::Node> &local_pt) const
+{
+ // rotate the point coordinates
+ const std::size_t global_dim = global_coords.getDimension();
+ Eigen::VectorXd dx = Eigen::VectorXd::Zero(global_dim);
+ Eigen::VectorXd x_new = Eigen::VectorXd::Zero(3);
+ for(std::size_t i = 0; i < ele.getNNodes(); i++)
+ {
+ for (std::size_t j=0; j<global_dim; j++)
+ dx[j] = vec_pt[i].getCoords()[j];
+
+ x_new.head(global_dim) = matR2local * dx;
+ local_pt[i] = MeshLib::Node(x_new.data());
+ }
+};
+
+void ElementCoordinatesMappingLocal::getRotationMatrixToGlobal(
+ const Element &e,
+ const CoordinateSystem &global_coords,
+ const std::vector<MeshLib::Node> &vec_pt,
+ RotationMatrix &matR) const
+{
+ const std::size_t global_dim = global_coords.getDimension();
+
+ // compute R in x=R*x' where x is original coordinates and x' is local coordinates
+ matR = RotationMatrix::Zero(global_dim, global_dim);
+ if (global_dim == e.getDimension()) {
+ matR = RotationMatrix::Identity(global_dim, global_dim);
+ } else if (e.getDimension() == 1) {
+ MathLib::Vector3 xx(vec_pt[0], vec_pt[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) {
+ std::vector<MathLib::Point3d*> pnts;
+ for (auto& p : vec_pt)
+ pnts.push_back(const_cast<MeshLib::Node*>(&p));
+
+ // get plane normal
+ MathLib::Vector3 plane_normal;
+ double d;
+ GeoLib::getNewellPlane (pnts, plane_normal, d);
+ //std::cout << "pn=" << plane_normal << std::endl;
+ // compute a rotation matrix to XY
+ MathLib::DenseMatrix<double> matToXY(3,3,0.0);
+ GeoLib::computeRotationMatrixToXY(plane_normal, matToXY);
+ // set a transposed matrix
+ for (size_t i=0; i<global_dim; ++i)
+ for (size_t j=0; j<global_dim; ++j)
+ matR(i, j) = matToXY(j,i);
+ }
+
+}
+
+} // MeshLib
diff --git a/MeshLib/ElementCoordinatesMappingLocal.h b/MeshLib/ElementCoordinatesMappingLocal.h
new file mode 100644
index 00000000000..974ce68c009
--- /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 "MeshLib/Elements/Element.h"
+#include "MeshLib/CoordinateSystem.h"
+
+namespace MeshLib
+{
+
+#ifdef OGS_USE_EIGEN
+typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> 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 &_point_vec[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 Element &e, const CoordinateSystem &coordinate_system,
+ const std::vector<MeshLib::Node> &vec_pt, const RotationMatrix &matR2local,
+ std::vector<MeshLib::Node> &local_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> _point_vec;
+ RotationMatrix _matR2global;
+};
+
+}
+
+#endif // ELEMENTCOORDINATESMAPPINGLOCAL_H_
diff --git a/Tests/MeshLib/TestCoordinatesMappingLocal.cpp b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
new file mode 100644
index 00000000000..30e223f2b95
--- /dev/null
+++ b/Tests/MeshLib/TestCoordinatesMappingLocal.cpp
@@ -0,0 +1,270 @@
+/**
+ * \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 dim = ElementType::dimension;
+ 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 dim = 2; //ElementType::dimension;
+ 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 40 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);
+ }
+
+};
+
+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;
+ }
+}
+
+// check if using the rotation matrix results in the original coordinates
+#define CHECK_COORDS(ele, mapping)\
+ for (unsigned ii=0; ii<ele->getNNodes(); ii++) {\
+ Eigen::Map<Eigen::Vector3d> local(const_cast<double*>(&(*mapping.getMappedCoordinates(ii))[0]));\
+ Eigen::Vector3d global(matR*local);\
+ const double eps(std::numeric_limits<double>::epsilon());\
+ ASSERT_ARRAY_NEAR(&(*ele->getNode(ii))[0], global.data(), 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[2*2] = {0, -1,
+ 1, 0}; //row major
+ 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}; //row major
+ 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[2*2] = {0.70710678118654757, -0.70710678118654757,
+ 0.70710678118654757, 0.70710678118654757}; //row major
+ 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}; //row major
+ 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}; //row major
+// // results when using GeoLib::ComputeRotationMatrixToXY2()
+// double exp_R[3*3] = { -1, 0, 0,
+// 0, 0, -1,
+// 0, -1, 0}; //row major
+
+ 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}; //row major
+// // results when using GeoLib::ComputeRotationMatrixToXY2()
+// double exp_R[3*3] = { 0, 0, 1,
+// -1, 0, 0,
+// 0, -1, 0}; //row major
+
+ 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}; //row major
+// // results when using GeoLib::ComputeRotationMatrixToXY2()
+// double exp_R[3*3] = { -1, 0, 0,
+// 0, -0.70710678118654757, -0.70710678118654757,
+// 0, -0.70710678118654757, 0.70710678118654757}; //row major
+
+ const double eps(std::numeric_limits<double>::epsilon());
+ ASSERT_ARRAY_NEAR(exp_R, matR.data(), matR.size(), eps);
+ CHECK_COORDS(ele,mapping);
+}
+
--
GitLab