Merge pull request #2034 from endJunction/FixBoundaryElementeSearchHex20Bug

Fix boundary element search hex20 bug

MeshGeoToolsLib/MeshNodeSearcher.cpp

@@ -137,12 +137,12 @@ MeshNodesAlongSurface& MeshNodeSearcher::getMeshNodesAlongSurface(
     std::vector<MeshNodesAlongSurface*>::const_iterator it(_mesh_nodes_along_surfaces.begin());
     for (; it != _mesh_nodes_along_surfaces.end(); ++it) {
         if (&(*it)->getSurface() == &sfc) {
-            // we calculated mesh nodes for this polyline already
+            // we calculated mesh nodes on this surface already
             return *(*it);
         }
     }
 
-    // compute nodes (and supporting points) along polyline
+    // compute nodes (and supporting points) on surface
     _mesh_nodes_along_surfaces.push_back(
         new MeshNodesAlongSurface(_mesh,
                                   sfc,

MeshGeoToolsLib/MeshNodesAlongSurface.h

@@ -39,14 +39,15 @@ public:
      * GeoLib::Surface object within a given search radius.
      * @param mesh Mesh the search will be performed on.
      * @param sfc Along the GeoLib::Surface sfc the mesh nodes are searched.
-     * @param epsilon Euclidean distance tolerance value. Is the distance
+     * @param epsilon_radius Euclidean distance tolerance value. Is the distance
      * between a mesh node and the surface smaller than that value it is a mesh
      * node near the surface.
      * @param search_all_nodes switch between searching all mesh nodes and
      * searching the base nodes.
      */
     MeshNodesAlongSurface(MeshLib::Mesh const& mesh, GeoLib::Surface const& sfc,
-                          double epsilon, SearchAllNodes search_all_nodes);
+                          double epsilon_radius,
+                          SearchAllNodes search_all_nodes);
 
     /// return the mesh object
     MeshLib::Mesh const& getMesh() const;

MeshLib/Elements/HexRule20.cpp

@@ -18,15 +18,13 @@
 #include "Line.h"
 
 namespace MeshLib {
-
-const unsigned HexRule20::face_nodes[6][8] =
-{
-    {0, 3, 2, 1, 11, 10,  9,  8}, // Face 0
-    {0, 1, 5, 4,  8, 17, 12, 11}, // Face 1
-    {1, 2, 6, 5,  9, 18, 13, 17}, // Face 2
-    {2, 3, 7, 6, 10, 19, 14, 18}, // Face 3
-    {3, 0, 4, 7, 11, 16, 15, 19}, // Face 4
-    {4, 5, 6, 7, 12, 13, 14, 15}  // Face 5
+const unsigned HexRule20::face_nodes[6][8] = {
+    {0, 3, 2, 1, 11, 10, 9, 8},    // Face 0
+    {0, 1, 5, 4, 8, 17, 12, 16},   // Face 1
+    {1, 2, 6, 5, 9, 18, 13, 17},   // Face 2
+    {2, 3, 7, 6, 10, 19, 14, 18},  // Face 3
+    {3, 0, 4, 7, 11, 16, 15, 19},  // Face 4
+    {4, 5, 6, 7, 12, 13, 14, 15}   // Face 5
 };
 
 const unsigned HexRule20::edge_nodes[12][3] =

MeshLib/MeshSearch/ElementSearch.cpp

@@ -11,8 +11,9 @@
 
 #include <logog/include/logog.hpp>
 
-#include "MeshLib/Node.h"
+#include "BaseLib/makeVectorUnique.h"
 #include "MeshLib/Elements/Element.h"
+#include "MeshLib/Node.h"
 
 namespace MeshLib {
 
@@ -77,9 +78,8 @@ std::size_t ElementSearch::searchByNodeIDs(const std::vector<std::size_t> &nodes
             connected_elements.push_back(e->getID());
         }
     }
-    std::sort(connected_elements.begin(), connected_elements.end());
-    auto it = std::unique(connected_elements.begin(), connected_elements.end());
-    connected_elements.resize(std::distance(connected_elements.begin(),it));
+
+    BaseLib::makeVectorUnique(connected_elements);
 
     this->updateUnion(connected_elements);
     return connected_elements.size();

Tests/MeshLib/TestBoundaryElementSearch.cpp

@@ -14,14 +14,15 @@
 #include "GeoLib/Polyline.h"
 #include "GeoLib/Surface.h"
 
-#include "MeshLib/Mesh.h"
-#include "MeshLib/Node.h"
+#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+#include "MeshGeoToolsLib/HeuristicSearchLength.h"
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
 #include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/MeshGenerators/QuadraticMeshGenerator.h"
 #include "MeshLib/MeshSearch/NodeSearch.h"
-#include "MeshGeoToolsLib/MeshNodeSearcher.h"
-#include "MeshGeoToolsLib/HeuristicSearchLength.h"
-#include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
+#include "MeshLib/Node.h"
 
 using namespace MeshLib;
 
@@ -33,10 +34,17 @@ public:
         _quad_mesh(MeshGenerator::generateRegularQuadMesh(_geometric_size, _number_of_subdivisions_per_direction))
     {}
 
+    ~MeshLibBoundaryElementSearchInSimpleQuadMesh() override
+    {
+        for (auto p : _pnts)
+            delete p;
+    }
+
 protected:
     const double _geometric_size;
     const std::size_t _number_of_subdivisions_per_direction;
     std::unique_ptr<Mesh> _quad_mesh;
+    std::vector<GeoLib::Point*> _pnts;
 };
 
 class MeshLibBoundaryElementSearchInSimpleHexMesh : public testing::Test
@@ -45,17 +53,34 @@ public:
     MeshLibBoundaryElementSearchInSimpleHexMesh() :
         _geometric_size(10.0), _number_of_subdivisions_per_direction(10),
         _hex_mesh(MeshGenerator::generateRegularHexMesh(_geometric_size, _number_of_subdivisions_per_direction))
-    {}
+    {
+        // Points for the surfaces. Corners of a cube.
+        _pnts.push_back(new GeoLib::Point(0.0, 0.0, 0.0));
+        _pnts.push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
+        _pnts.push_back(
+            new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
+        _pnts.push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
+        _pnts.push_back(
+            new GeoLib::Point(_geometric_size, 0.0, _geometric_size));
+        _pnts.push_back(new GeoLib::Point(0.0, 0.0, _geometric_size));
+    }
 
     ~MeshLibBoundaryElementSearchInSimpleHexMesh() override
     {
-        delete _hex_mesh;
+        for (auto p : _pnts)
+            delete p;
     }
 
+    void TestHexSurfacesBoundaryElementSearcher(
+        Mesh const& mesh,
+        std::size_t const n_nodes_2d,
+        std::size_t const n_eles_2d) const;
+
 protected:
     const double _geometric_size;
     const std::size_t _number_of_subdivisions_per_direction;
-    Mesh* _hex_mesh;
+    std::unique_ptr<Mesh> _hex_mesh;
+    std::vector<GeoLib::Point*> _pnts;
 };
 
 TEST_F(MeshLibBoundaryElementSearchInSimpleQuadMesh, PolylineSearch)
@@ -65,12 +90,11 @@ TEST_F(MeshLibBoundaryElementSearchInSimpleQuadMesh, PolylineSearch)
     const unsigned n_nodes_per_dir = _number_of_subdivisions_per_direction + 1;
 
     // create a polyline (closed)
-    std::vector<GeoLib::Point*> pnts;
-    pnts.push_back(new GeoLib::Point(0.0, 0.0, 0.0));
-    pnts.push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
-    pnts.push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
-    pnts.push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
-    GeoLib::Polyline ply0(pnts);
+    _pnts.push_back(new GeoLib::Point(0.0, 0.0, 0.0));
+    _pnts.push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
+    _pnts.push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
+    _pnts.push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
+    GeoLib::Polyline ply0(_pnts);
     ply0.addPoint(0);
     ply0.addPoint(1);
     ply0.addPoint(2);
@@ -110,68 +134,82 @@ TEST_F(MeshLibBoundaryElementSearchInSimpleQuadMesh, PolylineSearch)
         ASSERT_EQ(n_nodes_per_dir*(n_nodes_per_dir-1)-n_nodes_per_dir*i, edge3->getNodeIndex(0));
         ASSERT_EQ(n_nodes_per_dir*(n_nodes_per_dir-1)-n_nodes_per_dir*(i+1), edge3->getNodeIndex(1));
     }
-
-    for (auto p : pnts)
-        delete p;
 }
 
-TEST_F(MeshLibBoundaryElementSearchInSimpleHexMesh, SurfaceSearch)
+template <typename ElementIterator>
+double computeAreaOfFaceElements(ElementIterator first, ElementIterator last)
 {
-    ASSERT_TRUE(_hex_mesh != nullptr);
-    const std::size_t n_nodes_1d = _number_of_subdivisions_per_direction + 1;
-    const std::size_t n_nodes_2d = n_nodes_1d * n_nodes_1d;
-    const std::size_t n_eles_2d = (n_nodes_1d-1) * (n_nodes_1d-1);
-
-    // create bottom and front surfaces of a cubic
-    std::vector<GeoLib::Point*> pnts;
-    pnts.push_back(new GeoLib::Point(0.0, 0.0, 0.0));
-    pnts.push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
-    pnts.push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
-    pnts.push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
-    pnts.push_back(new GeoLib::Point(_geometric_size, 0.0, _geometric_size));
-    pnts.push_back(new GeoLib::Point(0.0, 0.0, _geometric_size));
-
-    GeoLib::Surface sfc_bottom(pnts);
-    sfc_bottom.addTriangle(0, 1, 2);
-    sfc_bottom.addTriangle(0, 2, 3);
-
-    GeoLib::Surface sfc_front(pnts);
-    sfc_front.addTriangle(0, 1, 4);
-    sfc_front.addTriangle(0, 4, 5);
+    return std::accumulate(first, last, 0.0, [](double v, MeshLib::Element* e) {
+        return v + e->getContent();
+    });
+}
 
-    // perform search on the bottom surface
+void MeshLibBoundaryElementSearchInSimpleHexMesh::
+    TestHexSurfacesBoundaryElementSearcher(Mesh const& mesh,
+                                           std::size_t const n_nodes_2d,
+                                           std::size_t const n_eles_2d) const
+{
     MeshGeoToolsLib::MeshNodeSearcher mesh_node_searcher(
-        *_hex_mesh,
-        std::make_unique<MeshGeoToolsLib::HeuristicSearchLength>(*_hex_mesh),
+        mesh,
+        std::make_unique<MeshGeoToolsLib::SearchLength>(),
         MeshGeoToolsLib::SearchAllNodes::Yes);
     MeshGeoToolsLib::BoundaryElementsSearcher boundary_element_searcher(
-        *_hex_mesh, mesh_node_searcher);
+        mesh, mesh_node_searcher);
+
+    // perform search on the bottom surface
+    GeoLib::Surface sfc_bottom(_pnts);
+    sfc_bottom.addTriangle(0, 1, 2);
+    sfc_bottom.addTriangle(0, 2, 3);
+
     std::vector<MeshLib::Element*> const& found_faces_sfc_b(
         boundary_element_searcher.getBoundaryElements(sfc_bottom));
     ASSERT_EQ(n_eles_2d, found_faces_sfc_b.size());
-    double sum_area_b = std::accumulate(
-        found_faces_sfc_b.begin(), found_faces_sfc_b.end(), 0.0,
-        [](double v, MeshLib::Element* e) { return v + e->getContent(); });
-    ASSERT_EQ(_geometric_size*_geometric_size, sum_area_b);
+    ASSERT_EQ(_geometric_size * _geometric_size,
+              computeAreaOfFaceElements(found_faces_sfc_b.begin(),
+                                        found_faces_sfc_b.end()));
     auto connected_nodes_b = MeshLib::getUniqueNodes(found_faces_sfc_b);
     ASSERT_EQ(n_nodes_2d, connected_nodes_b.size());
     for (auto node : connected_nodes_b)
-        ASSERT_EQ(0.0, (*node)[2]); // check z coordinates
+        ASSERT_EQ(0.0, (*node)[2]);  // check z coordinates
 
     // perform search on the front surface
+    GeoLib::Surface sfc_front(_pnts);
+    sfc_front.addTriangle(0, 1, 4);
+    sfc_front.addTriangle(0, 4, 5);
+
     std::vector<MeshLib::Element*> const& found_faces_sfc_f(
         boundary_element_searcher.getBoundaryElements(sfc_front));
     ASSERT_EQ(n_eles_2d, found_faces_sfc_f.size());
-    double sum_area_f = std::accumulate(
-        found_faces_sfc_f.begin(), found_faces_sfc_f.end(), 0.0,
-        [](double v, MeshLib::Element* e) { return v + e->getContent(); });
-    ASSERT_EQ(_geometric_size*_geometric_size, sum_area_f);
+    ASSERT_EQ(_geometric_size * _geometric_size,
+              computeAreaOfFaceElements(found_faces_sfc_f.begin(),
+                                        found_faces_sfc_f.end()));
     auto connected_nodes_f = MeshLib::getUniqueNodes(found_faces_sfc_f);
     ASSERT_EQ(n_nodes_2d, connected_nodes_f.size());
     for (auto node : connected_nodes_f)
-        ASSERT_EQ(0.0, (*node)[1]); // check y coordinates
+        ASSERT_EQ(0.0, (*node)[1]);  // check y coordinates
+}
+
+TEST_F(MeshLibBoundaryElementSearchInSimpleHexMesh, SurfaceSearch)
+{
+    ASSERT_TRUE(_hex_mesh != nullptr);
+    const std::size_t& s = _number_of_subdivisions_per_direction;
+    const std::size_t n_nodes_2d = (s + 1) * (s + 1);
+    const std::size_t n_eles_2d = s * s;
 
-    for (auto p : pnts)
-        delete p;
+    TestHexSurfacesBoundaryElementSearcher(*_hex_mesh, n_nodes_2d, n_eles_2d);
 }
 
+// This is identical to the above
+// MeshLibBoundaryElementSearchInSimpleHexMesh.SurfaceSearch test but
+// creates a quadratic mesh from the original hex mesh.
+TEST_F(MeshLibBoundaryElementSearchInSimpleHexMesh, QuadElementsSurfaceSearch)
+{
+    ASSERT_TRUE(_hex_mesh != nullptr);
+    auto mesh = MeshLib::createQuadraticOrderMesh(*_hex_mesh);
+
+    const std::size_t& s = _number_of_subdivisions_per_direction;
+    const std::size_t n_nodes_2d = (s + 1) * (3 * s + 1);
+    const std::size_t n_eles_2d = s * s;
+
+    TestHexSurfacesBoundaryElementSearcher(*mesh, n_nodes_2d, n_eles_2d);
+}