add ElementCoordinatesMappingLocal class

MeshLib/ElementCoordinatesMappingLocal.cpp

+/**
+ * \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

MeshLib/ElementCoordinatesMappingLocal.h

+/**
+ * \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_

Tests/MeshLib/TestCoordinatesMappingLocal.cpp

+/**
+ * \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);
+}
+