Merge pull request #2519 from endJunction/ConvertToLinearMeshTests

Convert to linear mesh tests

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

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

Tests/MeshLib/ConvertToLinearMesh.cpp

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