diff --git a/MeshLib/MeshEditing/ConvertToLinearMesh.cpp b/MeshLib/MeshEditing/ConvertToLinearMesh.cpp
index 06cc122804a5646865ec1f301d4e15f460eba604..d06d423e26520aafd0946d158802b3f7c01344c9 100644
--- a/MeshLib/MeshEditing/ConvertToLinearMesh.cpp
+++ b/MeshLib/MeshEditing/ConvertToLinearMesh.cpp
@@ -10,9 +10,9 @@
#include "ConvertToLinearMesh.h"
#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Elements/Hex.h"
#include "MeshLib/Elements/Line.h"
#include "MeshLib/Elements/Quad.h"
-#include "MeshLib/Elements/Hex.h"
#include "MeshLib/Elements/Tet.h"
#include "MeshLib/Elements/Tri.h"
#include "MeshLib/Elements/Utils.h"
@@ -22,33 +22,43 @@
#include "MeshLib/Properties.h"
#include "MeshLib/PropertyVector.h"
-
namespace MeshLib
{
-
namespace
{
-
template <typename T_ELEMENT>
T_ELEMENT* createLinearElement(MeshLib::Element const* e,
- std::vector<MeshLib::Node*> const& vec_new_nodes)
+ std::vector<MeshLib::Node*> const& vec_new_nodes)
{
auto const n_base_nodes = T_ELEMENT::n_base_nodes;
auto** nodes = new MeshLib::Node*[n_base_nodes];
for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
{
- nodes[i] =
- const_cast<MeshLib::Node*>(vec_new_nodes[e->getNode(i)->getID()]);
+ auto const it = find_if(
+ begin(vec_new_nodes), end(vec_new_nodes),
+ [node_i = e->getNode(i)](Node* const new_node) {
+ return *node_i ==
+ *new_node; // coordinate comparison up to epsilon
+ });
+ if (it == end(vec_new_nodes))
+ {
+ OGS_FATAL(
+ "A base node %d (with original global node id %d) not found in "
+ "the list for element %d.",
+ i, e->getNode(i)->getID(), e->getID());
+ }
+ nodes[i] = const_cast<MeshLib::Node*>(*it);
}
return new T_ELEMENT(nodes);
}
-} // unnamed namespace
+} // unnamed namespace
-
-std::unique_ptr<MeshLib::Mesh> convertToLinearMesh(MeshLib::Mesh const& org_mesh, std::string const& new_mesh_name)
+std::unique_ptr<MeshLib::Mesh> convertToLinearMesh(
+ MeshLib::Mesh const& org_mesh, std::string const& new_mesh_name)
{
- std::vector<MeshLib::Node*> vec_new_nodes = MeshLib::copyNodeVector(MeshLib::getBaseNodes(org_mesh.getElements()));
+ std::vector<MeshLib::Node*> vec_new_nodes =
+ MeshLib::copyNodeVector(MeshLib::getBaseNodes(org_mesh.getElements()));
// create new elements with the quadratic nodes
std::vector<MeshLib::Element*> vec_new_eles;
@@ -56,32 +66,33 @@ std::unique_ptr<MeshLib::Mesh> convertToLinearMesh(MeshLib::Mesh const& org_mesh
{
if (e->getCellType() == MeshLib::CellType::LINE3)
{
- vec_new_eles.push_back(createLinearElement<MeshLib::Line>(
- e, vec_new_nodes));
+ vec_new_eles.push_back(
+ createLinearElement<MeshLib::Line>(e, vec_new_nodes));
}
else if (e->getCellType() == MeshLib::CellType::QUAD8)
{
- vec_new_eles.push_back(createLinearElement<MeshLib::Quad>(
- e, vec_new_nodes));
+ vec_new_eles.push_back(
+ createLinearElement<MeshLib::Quad>(e, vec_new_nodes));
}
else if (e->getCellType() == MeshLib::CellType::TRI6)
{
- vec_new_eles.push_back(createLinearElement<MeshLib::Tri>(
- e, vec_new_nodes));
+ vec_new_eles.push_back(
+ createLinearElement<MeshLib::Tri>(e, vec_new_nodes));
}
else if (e->getCellType() == MeshLib::CellType::HEX20)
{
- vec_new_eles.push_back(createLinearElement<MeshLib::Hex>(
- e, vec_new_nodes));
+ vec_new_eles.push_back(
+ createLinearElement<MeshLib::Hex>(e, vec_new_nodes));
}
else if (e->getCellType() == MeshLib::CellType::TET10)
{
- vec_new_eles.push_back(createLinearElement<MeshLib::Tet>(
- e, vec_new_nodes));
+ vec_new_eles.push_back(
+ createLinearElement<MeshLib::Tet>(e, vec_new_nodes));
}
else
{
- OGS_FATAL("Mesh element type %s is not supported", MeshLib::CellType2String(e->getCellType()).c_str());
+ OGS_FATAL("Mesh element type %s is not supported",
+ MeshLib::CellType2String(e->getCellType()).c_str());
}
}
@@ -112,7 +123,7 @@ std::unique_ptr<MeshLib::Mesh> convertToLinearMesh(MeshLib::Mesh const& org_mesh
new_prop->resize(new_mesh->getNumberOfNodes() * n_src_comp);
// copy only base node values
- for (unsigned i=0; i<org_mesh.getNumberOfBaseNodes(); i++)
+ for (unsigned i = 0; i < org_mesh.getNumberOfBaseNodes(); i++)
{
for (int j = 0; j < n_src_comp; j++)
{
@@ -125,5 +136,4 @@ std::unique_ptr<MeshLib::Mesh> convertToLinearMesh(MeshLib::Mesh const& org_mesh
return new_mesh;
}
-} // end namespace MeshLib
-
+} // end namespace MeshLib
diff --git a/MeshLib/MeshEditing/ConvertToLinearMesh.h b/MeshLib/MeshEditing/ConvertToLinearMesh.h
index 9093004616ca2c4e4fb46c1066d59cda2f4c6b28..92253690b5bfa3937af5398881f9f3926637cb60 100644
--- a/MeshLib/MeshEditing/ConvertToLinearMesh.h
+++ b/MeshLib/MeshEditing/ConvertToLinearMesh.h
@@ -16,10 +16,9 @@
namespace MeshLib
{
-
/// Converts a non-linear mesh to a linear mesh. All the mesh properties will
/// be copied except for entries for non-linear nodes.
std::unique_ptr<MeshLib::Mesh> convertToLinearMesh(
const MeshLib::Mesh& mesh, const std::string& new_mesh_name);
-} // end namespace MeshLib
+} // end namespace MeshLib
diff --git a/Tests/MeshLib/ConvertToLinearMesh.cpp b/Tests/MeshLib/ConvertToLinearMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..068bea3f6f0024b8e30cfab2b7cae4d6c5633fb2
--- /dev/null
+++ b/Tests/MeshLib/ConvertToLinearMesh.cpp
@@ -0,0 +1,324 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <algorithm>
+#include <boost/variant.hpp>
+#include <numeric>
+#include <random>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshEditing/ConvertToLinearMesh.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/MeshGenerators/QuadraticMeshGenerator.h"
+#include "MeshLib/Node.h"
+
+using namespace MeshLib;
+using namespace std::string_literals;
+
+class ConvertToLinearMesh : public ::testing::Test
+{
+public:
+ ConvertToLinearMesh()
+ {
+ linear_mesh.reset(
+ MeshGenerator::generateRegularHexMesh(1.0, mesh_size));
+ quadratic_mesh = createQuadraticOrderMesh(*linear_mesh);
+ }
+
+ static constexpr int mesh_size = 5;
+
+ std::unique_ptr<Mesh> linear_mesh;
+ std::unique_ptr<Mesh> quadratic_mesh;
+};
+
+// Returns the (inverse) permutation vector of b_nodes to a_nodes, or an error
+// message.
+boost::variant<std::vector<std::size_t>, std::string> compareNodes(
+ std::vector<Node*> const& a_nodes, std::vector<Node*> const& b_nodes)
+{
+ // For each 'a' mesh node find a corresponding 'b' mesh node, not checking
+ // the ids, but store the id's in a map.
+ std::vector<std::size_t> b_node_id(a_nodes.size());
+
+ std::size_t const nnodes = a_nodes.size();
+ for (std::size_t i = 0; i < nnodes; ++i)
+ {
+ Node* a_node = a_nodes[i];
+ auto const b_it =
+ std::find_if(begin(b_nodes), end(b_nodes), [&](Node* const b_node) {
+ return *a_node ==
+ *b_node; // coordinate comparison up to epsilon
+ });
+ if (b_it == end(b_nodes))
+ {
+ return "Could not find a matching node for node " +
+ std::to_string(a_node->getID()) + ".";
+ }
+
+ b_node_id[i] = distance(begin(b_nodes), b_it);
+ }
+
+ return b_node_id;
+}
+
+// Two elements are equal if their corresponding node coordinates are equal up
+// to a (shift) permutation.
+// Returns a string in case of error containing a descriptive message.
+boost::optional<std::string> equal(Element const& a, Element const& b)
+{
+ if (typeid(a) != typeid(b)) // => (a.nodes.size == b.nodes.size)
+ {
+ return "Elements have different types."s;
+ }
+
+ auto const a_nodes = a.getNodes();
+ auto const b_nodes = b.getNodes();
+ auto const nnodes = a.getNumberOfNodes();
+
+ std::size_t first_b_node = 0;
+ while (first_b_node < nnodes)
+ {
+ if (*a_nodes[0] == *b_nodes[first_b_node])
+ break;
+ ++first_b_node;
+ }
+ if (first_b_node == nnodes)
+ {
+ return "Did not find first node of the element " +
+ std::to_string(a.getID()) + " in the nodes of the element " +
+ std::to_string(b.getID()) + ".";
+ }
+
+ for (std::size_t i = 1; i < nnodes; ++i)
+ {
+ if (!(*a_nodes[i] == *b_nodes[(first_b_node + i) % nnodes]))
+ {
+ return "node " + std::to_string(i) + " is not equal to the node " +
+ std::to_string((first_b_node + i) % nnodes) + ".";
+ }
+ }
+
+ return {}; // All nodes are equal.
+}
+
+// Returns a vector filled with numbers [0, size) randomly permuted.
+std::vector<std::size_t> generateRandomPermutationVector(std::size_t const size)
+{
+ std::random_device rd;
+ std::mt19937 random_generator(rd());
+
+ std::vector<std::size_t> permutation(size);
+ std::iota(begin(permutation), end(permutation), 0);
+ std::shuffle(begin(permutation), end(permutation), random_generator);
+ return permutation;
+}
+
+// Test the inverse permutation for being the inverse: p^-1(p(.)) = Id.
+// In case of failure an error message is returned.
+boost::optional<std::string> inversePermutationIdentityTest(
+ std::vector<std::size_t> const& permutation,
+ std::vector<std::size_t> const& inverse_permutation)
+{
+ if (permutation.size() != inverse_permutation.size())
+ {
+ return "Inverse permutation size differs from the permutation size: " +
+ std::to_string(inverse_permutation.size()) + " vs. " +
+ std::to_string(permutation.size()) + ".";
+ }
+
+ for (std::size_t i = 0; i < permutation.size(); ++i)
+ {
+ auto const p = permutation[i];
+ if (inverse_permutation[p] == i)
+ {
+ continue;
+ }
+ return "Inverse permutation element " + std::to_string(p) +
+ " does not point to " + std::to_string(i) +
+ "-th element but to " + std::to_string(inverse_permutation[p]) +
+ ".";
+ }
+ return {};
+}
+
+// Create new nodes in permuted order with new ids.
+std::vector<Node*> permuteNodes(std::vector<std::size_t> const& permutation,
+ std::vector<Node*> const& nodes)
+{
+ std::vector<Node*> permuted_nodes;
+ for (std::size_t i = 0; i < permutation.size(); ++i)
+ {
+ auto const node_p_i = nodes[permutation[i]];
+ permuted_nodes.push_back(new Node{node_p_i->getCoords(), i});
+ }
+ return permuted_nodes;
+}
+
+// Tests the nodes of two meshes and returns a string in case of error with
+// message, otherwise a permutation vector as in compareNodes().
+boost::variant<std::vector<std::size_t>, std::string> compareNodeVectors(
+ std::vector<Node*> const& a, std::vector<Node*> const& b)
+{
+ if (a.size() != b.size())
+ {
+ return "There are different number of nodes: " +
+ std::to_string(a.size()) + " and " + std::to_string(b.size()) +
+ ".";
+ }
+
+ auto const compare_nodes_result = compareNodes(a, b);
+ if (compare_nodes_result.type() == typeid(std::string))
+ {
+ return "Node comparison failed: " +
+ boost::get<std::string>(compare_nodes_result);
+ }
+
+ return boost::get<std::vector<std::size_t>>(compare_nodes_result);
+}
+
+TEST_F(ConvertToLinearMesh, GeneratedHexMeshRandomizedNodes)
+{
+ ASSERT_TRUE(linear_mesh != nullptr);
+ ASSERT_TRUE(quadratic_mesh != nullptr);
+
+ auto const permutation =
+ generateRandomPermutationVector(quadratic_mesh->getNumberOfNodes());
+
+ // Create new nodes in permuted order with new ids.
+ auto const permuted_quadratic_nodes =
+ permuteNodes(permutation, quadratic_mesh->getNodes());
+
+ //
+ // Test the permuted nodes and compute the inverse permutation map.
+ //
+ auto const compare_nodes_result = compareNodeVectors(
+ quadratic_mesh->getNodes(), permuted_quadratic_nodes);
+ ASSERT_FALSE(compare_nodes_result.type() == typeid(std::string))
+ << "Quadratic mesh nodes comparison with permuted quadratic nodes "
+ "failed.\n"
+ << boost::get<std::string>(compare_nodes_result);
+
+ auto const& inverse_permutation =
+ boost::get<std::vector<std::size_t>>(compare_nodes_result);
+ {
+ auto const result =
+ inversePermutationIdentityTest(permutation, inverse_permutation);
+ ASSERT_TRUE(result == boost::none)
+ << "Quadratic mesh nodes permutation test failed: " << *result;
+ }
+
+ auto clone_element_using_permuted_nodes = [&](Element* const e) {
+ Node** nodes = new Node*[e->getNumberOfNodes()];
+ for (std::size_t i = 0; i < e->getNumberOfNodes(); ++i)
+ {
+ auto const q_i = e->getNode(i)->getID();
+ auto const r = inverse_permutation[q_i];
+ nodes[i] = permuted_quadratic_nodes[r];
+ }
+ return e->clone(nodes, e->getID());
+ };
+
+ // Create new elements in the same order, but using new nodes.
+ std::vector<Element*> new_quadratic_elements;
+ std::transform(begin(quadratic_mesh->getElements()),
+ end(quadratic_mesh->getElements()),
+ back_inserter(new_quadratic_elements),
+ clone_element_using_permuted_nodes);
+
+ Mesh permuted_nodes_quadratic_mesh{"permuted_quadratic_mesh",
+ permuted_quadratic_nodes,
+ new_quadratic_elements};
+ //
+ // Test the quadratic meshes first.
+ //
+ {
+ ASSERT_EQ(quadratic_mesh->getNumberOfElements(),
+ permuted_nodes_quadratic_mesh.getNumberOfElements());
+
+ auto const nelements = quadratic_mesh->getNumberOfElements();
+ for (std::size_t i = 0; i < nelements; ++i)
+ {
+ auto const& element_a = *quadratic_mesh->getElement(i);
+ auto const& element_b =
+ *permuted_nodes_quadratic_mesh.getElement(i);
+ auto const elements_are_equal = equal(element_a, element_b);
+ ASSERT_TRUE(elements_are_equal == boost::none)
+ << *elements_are_equal << " For the element " << i << ".";
+ }
+ }
+
+ {
+ auto const converted_mesh = convertToLinearMesh(
+ permuted_nodes_quadratic_mesh, "back_to_linear_mesh");
+
+ //
+ // Two meshes are equal if all of their nodes have exactly same
+ // coordinates.
+ //
+ auto const compare_nodes_result = compareNodeVectors(
+ linear_mesh->getNodes(), converted_mesh->getNodes());
+ ASSERT_FALSE(compare_nodes_result.type() == typeid(std::string))
+ << "Linear mesh nodes comparison with converted linear mesh nodes "
+ "failed.\n"
+ << boost::get<std::string>(compare_nodes_result);
+ // TODO (naumov) check inverse permutation, but it requires the reduced
+ // to base nodes forward permutation vector first.
+
+ // Two meshes are equal if their elements share same nodes.
+ ASSERT_EQ(linear_mesh->getNumberOfElements(),
+ converted_mesh->getNumberOfElements());
+
+ auto const linear_mesh_nelements = linear_mesh->getNumberOfElements();
+ for (std::size_t i = 0; i < linear_mesh_nelements; ++i)
+ {
+ auto const& element_a = *linear_mesh->getElement(i);
+ auto const& element_b = *converted_mesh->getElement(i);
+ auto const elements_are_equal = equal(element_a, element_b);
+ ASSERT_TRUE(elements_are_equal == boost::none)
+ << *elements_are_equal << " For the element " << i << ".";
+ }
+ }
+}
+
+TEST_F(ConvertToLinearMesh, GeneratedHexMeshBackToLinear)
+{
+ ASSERT_TRUE(linear_mesh != nullptr);
+ ASSERT_TRUE(quadratic_mesh != nullptr);
+
+ auto const converted_mesh =
+ convertToLinearMesh(*quadratic_mesh, "back_to_linear_mesh");
+
+ //
+ // Two meshes are equal if all of their nodes have exactly same coordinates.
+ //
+ auto const compare_nodes_result =
+ compareNodeVectors(linear_mesh->getNodes(), converted_mesh->getNodes());
+ ASSERT_FALSE(compare_nodes_result.type() == typeid(std::string))
+ << "Linear mesh nodes comparison with converted linear mesh nodes "
+ "failed.\n"
+ << boost::get<std::string>(compare_nodes_result);
+
+ // Two meshes are equal if their elements share same nodes.
+ ASSERT_EQ(linear_mesh->getNumberOfElements(),
+ converted_mesh->getNumberOfElements());
+
+ auto const linear_mesh_nelements = linear_mesh->getNumberOfElements();
+ for (std::size_t i = 0; i < linear_mesh_nelements; ++i)
+ {
+ auto const& linear_mesh_element = *linear_mesh->getElement(i);
+ auto const& converted_mesh_element = *converted_mesh->getElement(i);
+ auto const elements_are_equal =
+ equal(linear_mesh_element, converted_mesh_element);
+ ASSERT_TRUE(elements_are_equal == boost::none)
+ << *elements_are_equal << " For the element " << i << ".";
+ }
+}