From 37d27f636a4c7772e95b9339c86c20a6df16c95c Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 18 May 2021 17:10:16 +0200
Subject: [PATCH] [UnitTest] Added a test for getElementRotationMatrices with
3D mesh
---
.../TestGetElementRotationMatrices.cpp | 133 ++++++++++++++++++
1 file changed, 133 insertions(+)
create mode 100644 Tests/MeshLib/TestGetElementRotationMatrices.cpp
diff --git a/Tests/MeshLib/TestGetElementRotationMatrices.cpp b/Tests/MeshLib/TestGetElementRotationMatrices.cpp
new file mode 100644
index 00000000000..685f3b082f4
--- /dev/null
+++ b/Tests/MeshLib/TestGetElementRotationMatrices.cpp
@@ -0,0 +1,133 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on May 18, 2021, 12:31 PM
+ */
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Elements/Elements.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshEnums.h"
+#include "MeshLib/Node.h"
+#include "MeshLib/Utils/GetElementRotationMatrices.h"
+#include "MeshLib/Utils/GetSpaceDimension.h"
+
+TEST(MeshLib, GetElementRotationMatrices3DMesh)
+{
+ // The memory of the nodes are allocated by new operator, and it is released
+ // in the destructor of MeshLib::Mesh.
+
+ // Construct a 3D mesh, which contains two inclined 2D triangle elements
+ // and an inclined line element.
+ std::vector<MeshLib::Node*> nodes(10);
+ nodes[0] = new MeshLib::Node(0.0, 0.0, 0.0, 0);
+ nodes[1] = new MeshLib::Node(1.0, 0.0, 0.0, 1);
+ nodes[2] = new MeshLib::Node(1.0, 1.0, 0.0, 2);
+ nodes[3] = new MeshLib::Node(0.0, 1.0, 0.0, 3);
+
+ nodes[4] = new MeshLib::Node(0.0, 0.0, 1.0, 4);
+ nodes[5] = new MeshLib::Node(1.0, 0.0, 1.0, 5);
+ nodes[6] = new MeshLib::Node(1.0, 1.0, 1.0, 6);
+ nodes[7] = new MeshLib::Node(0.0, 1.0, 1.0, 7);
+
+ nodes[8] = new MeshLib::Node(0.0, 0.5, 1.5, 8);
+ nodes[9] = new MeshLib::Node(1.0, 0.5, 1.5, 9);
+
+ std::vector<MeshLib::Element*> elements;
+
+ // One hexahedral element:
+ std::array<MeshLib::Node*, 8> hex_element_nodes;
+ std::copy_n(nodes.begin(), 8, hex_element_nodes.begin());
+ elements.push_back(new MeshLib::Hex(hex_element_nodes));
+
+ // Two inclined triangle elements:
+ std::array<MeshLib::Node*, 3> tri_element_nodes{nodes[6], nodes[7],
+ nodes[8]};
+ elements.push_back(new MeshLib::Tri(tri_element_nodes));
+ tri_element_nodes[0] = nodes[6];
+ tri_element_nodes[1] = nodes[8];
+ tri_element_nodes[2] = nodes[9];
+ elements.push_back(new MeshLib::Tri(tri_element_nodes));
+
+ // One inclined line element:
+ std::array<MeshLib::Node*, 2> line_element_nodes{nodes[6], nodes[8]};
+ elements.push_back(new MeshLib::Line(line_element_nodes));
+
+ std::vector<std::unique_ptr<MeshLib::Mesh>> meshes;
+ meshes.push_back(
+ std::make_unique<MeshLib::Mesh>("a_mesh", nodes, elements));
+
+ int const space_dimension = MeshLib::getSpaceDimension(nodes);
+ auto const element_rotation_matrices = MeshLib::getElementRotationMatrices(
+ space_dimension, meshes[0]->getDimension(), meshes[0]->getElements());
+
+ // First element, the hexahedral element has identity matrix:
+ EXPECT_EQ(9, element_rotation_matrices[0].size());
+
+ // Second and third elements, the inclined triangle elements on the
+ // same plane:
+ auto const& rotation_matrix_tri1 = element_rotation_matrices[1];
+ auto const& rotation_matrix_tri2 = element_rotation_matrices[2];
+
+ double const diff = (rotation_matrix_tri1 - rotation_matrix_tri2).norm();
+ ASSERT_LE(diff, 1e-10);
+
+ Eigen::VectorXd b(3);
+ b[0] = 0.0;
+ b[1] = 0.0;
+ b[2] = 1.0;
+
+ // Projection test
+ Eigen::VectorXd const b_local = rotation_matrix_tri1.transpose() * b;
+ EXPECT_EQ(elements[1]->getDimension(), b_local.size());
+ EXPECT_EQ(0.0, b_local[0]);
+ double const expected_b_local1 = -std::sqrt(2.0) / 2.0;
+ EXPECT_LE(std::fabs(b_local[1] - expected_b_local1), 1.e-16);
+
+ // Forth element, the inclined line element:
+ Eigen::VectorXd const b_local_1D =
+ element_rotation_matrices[3].transpose() * b;
+ EXPECT_EQ(elements[3]->getDimension(), b_local_1D.size());
+
+ double const* const x_6 = nodes[6]->getCoords();
+ double const* const x_8 = nodes[8]->getCoords();
+ // b_local_1D = |b| (x_8-x_6) * b/(|(x_8-x_6)| |b|)
+ double dx[3];
+ for (int i = 0; i < 3; i++)
+ {
+ dx[i] = x_8[i] - x_6[i];
+ }
+
+ double const dx_norm =
+ std::sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+ double const expected_b_local_1D = dx[2] / dx_norm;
+ EXPECT_LE(std::fabs(expected_b_local_1D - b_local_1D[0]), 1.e-16);
+
+ // Test of rotation of a local vector to the global system
+ Eigen::VectorXd const local_1D_to_global =
+ element_rotation_matrices[3] * b_local_1D;
+ EXPECT_EQ(meshes[0]->getDimension(), local_1D_to_global.size());
+
+ double const local_1D_to_global_norm =
+ std::sqrt(local_1D_to_global[0] * local_1D_to_global[0] +
+ local_1D_to_global[1] * local_1D_to_global[1] +
+ local_1D_to_global[2] * local_1D_to_global[2]);
+ for (int i = 0; i < 3; i++)
+ {
+ EXPECT_LE(std::fabs(dx[i] / dx_norm -
+ local_1D_to_global[i] / local_1D_to_global_norm),
+ 1.e-15);
+ }
+}
--
GitLab