From eae3943d29c1dcf1032bafec6b0c4fb36f3e26f8 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri, 14 May 2021 14:52:01 +0200
Subject: [PATCH] [MeshLib] Added a utility of getElementRotationMatrices
---
MeshLib/Utils/GetElementRotationMatrices.cpp | 48 ++++++++++++++++++++
MeshLib/Utils/GetElementRotationMatrices.h | 24 ++++++++++
MeshLib/Utils/GetSpaceDimension.cpp | 42 +++++++++++++++++
MeshLib/Utils/GetSpaceDimension.h | 45 ++++++++++++++++++
MeshLib/Utils/SetMeshSpaceDimension.cpp | 29 +-----------
5 files changed, 161 insertions(+), 27 deletions(-)
create mode 100644 MeshLib/Utils/GetElementRotationMatrices.cpp
create mode 100644 MeshLib/Utils/GetElementRotationMatrices.h
create mode 100644 MeshLib/Utils/GetSpaceDimension.cpp
create mode 100644 MeshLib/Utils/GetSpaceDimension.h
diff --git a/MeshLib/Utils/GetElementRotationMatrices.cpp b/MeshLib/Utils/GetElementRotationMatrices.cpp
new file mode 100644
index 00000000000..8ec79cc2380
--- /dev/null
+++ b/MeshLib/Utils/GetElementRotationMatrices.cpp
@@ -0,0 +1,48 @@
+/**
+ * \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 14, 2021, 2:38 PM
+ */
+
+#include "GetElementRotationMatrices.h"
+
+#include "GetSpaceDimension.h"
+#include "MeshLib/ElementCoordinatesMappingLocal.h"
+#include "MeshLib/Elements/Elements.h"
+#include "MeshLib/Mesh.h"
+
+namespace MeshLib
+{
+std::vector<Eigen::MatrixXd> getElementRotationMatrices(
+ int const space_dimension, int const mesh_dimension,
+ std::vector<Element*> const& elements)
+{
+ std::vector<Eigen::MatrixXd> element_rotation_matrices;
+ element_rotation_matrices.reserve(elements.size());
+ for (auto const& element : elements)
+ {
+ int const element_dimension = static_cast<int>(element->getDimension());
+ if (element_dimension == space_dimension)
+ {
+ element_rotation_matrices.emplace_back(Eigen::MatrixXd::Identity(
+ element_dimension, element_dimension));
+ }
+ else
+ {
+ MeshLib::ElementCoordinatesMappingLocal coordinates_mapping(
+ *element, mesh_dimension);
+
+ element_rotation_matrices.emplace_back(
+ coordinates_mapping.getRotationMatrixToGlobal().topLeftCorner(
+ space_dimension, element_dimension));
+ }
+ }
+
+ return element_rotation_matrices;
+}
+} // namespace MeshLib
diff --git a/MeshLib/Utils/GetElementRotationMatrices.h b/MeshLib/Utils/GetElementRotationMatrices.h
new file mode 100644
index 00000000000..e09be8c695b
--- /dev/null
+++ b/MeshLib/Utils/GetElementRotationMatrices.h
@@ -0,0 +1,24 @@
+/**
+ * \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 14, 2021, 2:38 PM
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace MeshLib
+{
+class Element;
+
+std::vector<Eigen::MatrixXd> getElementRotationMatrices(
+ int const space_dimension, int const mesh_dimension,
+ std::vector<Element*> const& elements);
+} // namespace MeshLib
diff --git a/MeshLib/Utils/GetSpaceDimension.cpp b/MeshLib/Utils/GetSpaceDimension.cpp
new file mode 100644
index 00000000000..6ebddd2fa41
--- /dev/null
+++ b/MeshLib/Utils/GetSpaceDimension.cpp
@@ -0,0 +1,42 @@
+/**
+ * \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 June 8, 2021, 12:16 PM
+ */
+
+#include "GetSpaceDimension.h"
+
+#include <algorithm>
+#include <array>
+#include <limits>
+
+#include "BaseLib/Error.h"
+#include "MeshLib/Node.h"
+
+namespace MeshLib
+{
+int getSpaceDimension(std::vector<Node*> const& nodes)
+{
+ std::array x_magnitude = {0.0, 0.0, 0.0};
+
+ double const* const x_ref = nodes[0]->getCoords();
+ for (auto const& node : nodes)
+ {
+ auto const x = node->getCoords();
+ for (int i = 0; i < 3; i++)
+ {
+ x_magnitude[i] += std::fabs(x[i] - x_ref[i]);
+ }
+ }
+
+ return static_cast<int>(std::count_if(
+ x_magnitude.begin(), x_magnitude.end(),
+ [](const double x_i_magnitude)
+ { return x_i_magnitude > std::numeric_limits<double>::epsilon(); }));
+}
+}; // namespace MeshLib
diff --git a/MeshLib/Utils/GetSpaceDimension.h b/MeshLib/Utils/GetSpaceDimension.h
new file mode 100644
index 00000000000..abd3fa76234
--- /dev/null
+++ b/MeshLib/Utils/GetSpaceDimension.h
@@ -0,0 +1,45 @@
+/**
+ * \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 June 8, 2021, 12:16 PM
+ */
+
+#pragma once
+
+#include <vector>
+
+namespace MeshLib
+{
+class Node;
+
+/**
+ * \brief Computes dimension of the embedding space containing the set of given
+ * points.
+ *
+ * The space dimension is computed by accounting the non-zero norms of
+ * \f$\mathbf x\f$, \f$\mathbf y\f$, \f$\mathbf z\f$, which are
+ * the coordinates of all nodes of the mesh. With this concept,
+ * the space dimension of a mesh is:
+ * - 1, if the mesh is 1D and and the mesh is parallel either to
+ * \f$x\f$, \f$y\f$ or to \f$z\f$ axis.
+ * - 2, if the mesh is 1D but it contains inclined elements on the origin
+ * coordinate plane of \f$x-y,\, y-z,\, \text{or }\, x-z\f$.
+ * That means the coordinates of all nodes are
+ * \f$(x, y, 0)\f$, \f$(x, 0, z)\f$ or \f$(0, y, z)\f$.
+ * - 3, if the mesh is 1D and but it is not on any origin
+ * coordinate plane of \f$x-y,\, y-z,\, \text{or }\, x-z\f$.
+ * - 2, if the mesh is 2D and it is on the origin
+ * coordinate plane of \f$x-y,\, y-z,\, \text{or }\, x-z\f$.
+ * - 2, if the mesh is 2D and it is on vertical or horizontal plane
+ * that is parallel to the the origin coordinate
+ * plane of \f$x-y,\, y-z,\, \text{or }\, x-z\f$ but with an offset.
+ * - 3, if the mesh contains inclined 2D elements.
+ * - 3, if the mesh contains 3D elements.
+ */
+int getSpaceDimension(std::vector<Node*> const& nodes);
+}; // namespace MeshLib
diff --git a/MeshLib/Utils/SetMeshSpaceDimension.cpp b/MeshLib/Utils/SetMeshSpaceDimension.cpp
index 35c908e2f66..ed28a84aad8 100644
--- a/MeshLib/Utils/SetMeshSpaceDimension.cpp
+++ b/MeshLib/Utils/SetMeshSpaceDimension.cpp
@@ -11,41 +11,16 @@
#include "SetMeshSpaceDimension.h"
-#include <algorithm>
-#include <array>
-#include <limits>
-
-#include "BaseLib/Error.h"
+#include "GetSpaceDimension.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Mesh.h"
-#include "MeshLib/Node.h"
namespace MeshLib
{
-unsigned getSpaceDimension(Mesh const& mesh)
-{
- std::array x_magnitude = {0.0, 0.0, 0.0};
-
- auto const nodes = mesh.getNodes();
- for (auto const& node : nodes)
- {
- auto const x = node->getCoords();
- for (int i = 0; i < 3; i++)
- {
- x_magnitude[i] += std::fabs(x[i]);
- }
- }
-
- return static_cast<unsigned>(std::count_if(
- x_magnitude.begin(), x_magnitude.end(), [](const double x_i_magnitude) {
- return x_i_magnitude > std::numeric_limits<double>::epsilon();
- }));
-}
-
void setMeshSpaceDimension(std::vector<std::unique_ptr<Mesh>> const& meshes)
{
// Get the space dimension from the bulk mesh:
- auto const space_dimension = getSpaceDimension(*meshes[0]);
+ auto const space_dimension = getSpaceDimension(meshes[0]->getNodes());
for (auto& mesh : meshes)
{
auto elements = mesh->getElements();
--
GitLab