mesh independency from neighboring input information (computed at con…

…struction time from element matrix). Tested on 1.5D, 2D, 2.5D, 3D geometries

fdaPDE/mesh/element.h

@@ -43,6 +43,8 @@ constexpr int ct_nnodes(const int M, const int R) {
 constexpr int ct_nvertices(const int M) { return M + 1; }
 // number of edges of an M-dimensional element
 constexpr int ct_nedges(const int M) { return (M * (M + 1)) / 2; }
+// number of faces of an M-dimensional element
+constexpr int ct_nneighbors(const int M) { return (M == 1) ? fdapde::Dynamic : (M + 1); }
 
 // A single mesh element. This object represents the main **geometrical** abstraction of a physical element.
 // No functional information is carried by instances of this object.

fdaPDE/mesh/mesh.h

@@ -51,7 +51,7 @@ template <int M, int N> struct neighboring_structure {
 // access to domain's triangulation, M: tangent space dimension, N: embedding space dimension
 template <int M, int N> class Mesh {
    private:
-    // coordinates of points costituting the vertices of mesh elements
+    // coordinates of points making the vertices of mesh elements
     DMatrix<double> nodes_ {};
     int n_nodes_ = 0;
     // identifiers of points (as row indexes in points_ matrix) composing each element, by row
@@ -65,14 +65,13 @@ template <int M, int N> class Mesh {
     std::vector<int> edges_ {};
     std::vector<int> edge_map_ {};   // the i-th row refers to the elements' identifiers insisting on the i-th edge
     int n_edges_ = 0;
-  
+
     // precomputed set of elements (cached for fast access)
     std::vector<Element<M, N>> cache_ {};
     void fill_cache();
    public:
     Mesh() = default;
-    Mesh(const DMatrix<double>& nodes, const DMatrix<int>& elements,
-	 const typename neighboring_structure<M, N>::type& neighbors, const DMatrix<int>& boundary);
+    Mesh(const DMatrix<double>& nodes, const DMatrix<int>& elements, const DMatrix<int>& boundary);
 
     // getters
     const Element<M, N>& element(int ID) const { return cache_[ID]; }
@@ -161,59 +160,103 @@ template <int M, int N> class Mesh {
         embedding_dimension = N,
         n_vertices = ct_nvertices(M),
         n_edges_per_element = ct_nedges(M),
-        n_neighbors = ct_nedges(M),
-        n_vertices_per_edge = M, // generalize wrt other dimensionalities
-	n_elements_per_edge = 2  // same here...
+        n_neighbors_per_element = ct_nneighbors(M),
+        n_vertices_per_edge = M,
+	n_elements_per_edge = 2  // is n_elements_per_face, change
     };
 };
 
 // implementative details
 
-// construct from raw matrices (NB: matrix enumeration is assumed to start from 0)
 template <int M, int N>
-Mesh<M, N>::Mesh(
-  const DMatrix<double>& nodes, const DMatrix<int>& elements,
-  const typename neighboring_structure<M, N>::type& neighbors, const DMatrix<int>& boundary) :
-    nodes_(nodes), neighbors_(neighbors), elements_(elements), boundary_(boundary) {
-    // store number of nodes and number of elements
-    n_nodes_ = nodes_.rows();
-    n_elements_ = elements_.rows();
+Mesh<M, N>::Mesh(const DMatrix<double>& nodes, const DMatrix<int>& elements, const DMatrix<int>& boundary) :
+  nodes_(nodes), elements_(elements), boundary_(boundary) {
+  // store number of nodes and number of elements
+  n_nodes_ = nodes_.rows();
+  n_elements_ = elements_.rows();
 
-    // compute mesh limits
-    for (size_t dim = 0; dim < N; ++dim) {
-        range_[dim].first = nodes_.col(dim).minCoeff();
-        range_[dim].second = nodes_.col(dim).maxCoeff();
-    }
-    // scan the whole mesh and precompute elements informations for fast access
-    fill_cache();
+  // compute mesh limits
+  for (size_t dim = 0; dim < N; ++dim) {
+    range_[dim].first = nodes_.col(dim).minCoeff();
+    range_[dim].second = nodes_.col(dim).maxCoeff();
+  }
 
     // compute edges informations
-    // edges_ are contigously stored in memory as a std::vector<int>, in a row major format
     auto edge_pattern = combinations<n_vertices_per_edge, n_vertices>();
     std::unordered_map<std::array<int, n_vertices_per_edge>, int, std_array_hash<int, n_vertices_per_edge>> visited;
-
     std::array<int, n_vertices_per_edge> edge;
+    // cycle over all elements
     for (int i = 0; i < n_elements_; ++i) {
-        for (int j = 0; j < edge_pattern.rows(); ++j) {
-            // construct edge
-            for (int k = 0; k < n_vertices_per_edge; ++k) { edge[k] = elements_(i, edge_pattern(j, k)); }
-            // check if edge already processed
-            std::sort(edge.begin(), edge.end());
-            auto it = visited.find(edge);
-            if (it != visited.end()) {
-                // free memory (only two elements share the same edge) and update edge to element structure
-                *(edge_map_.begin() + 2 * (it->second) + 1) = i;
-                visited.erase(it);
-            } else {
-                // store edge and update edge to element information
-                for (int k = 0; k < n_vertices_per_edge; ++k) { edges_.emplace_back(edge[k]); }
-                visited.insert({edge, n_edges_});
-                n_edges_++;
-                edge_map_.insert(edge_map_.end(), {i, -1});   // -1 flags no incident element at edge
-            }
-        }
+      for (int j = 0; j < edge_pattern.rows(); ++j) {
+	// construct edge
+	for (int k = 0; k < n_vertices_per_edge; ++k) { edge[k] = elements_(i, edge_pattern(j, k)); }
+	std::sort(edge.begin(), edge.end());    // normalize wrt node ordering
+	auto it = visited.find(edge);
+	if (it != visited.end()) {
+	  // free memory (two, and only two, elements share the same face) and update face to element bounding
+	  *(edge_map_.begin() + n_elements_per_edge * (it->second) + 1) = i;
+	  visited.erase(it);
+	} else {
+	  // store edge and update face to element bounding
+	  for (int k = 0; k < n_vertices_per_edge; ++k) { edges_.emplace_back(edge[k]); }
+	  visited.insert({edge, n_edges_});
+	  edge_map_.insert(edge_map_.end(), {i, -1});   // -1 flags no incident element at edge
+	  n_edges_++;
+	}
+      }
+    }
+
+
+  if constexpr (!is_linear_network<M, N>::value) {  
+    neighbors_ = DMatrix<int>::Constant(n_elements_, n_neighbors_per_element, -1);
+
+    // compute neighboring informations, such that the neighbor indexed i lies opposite to vertex i
+    for (auto it = edge_begin(); it != edge_end(); ++it) {
+      for (int i = 0; i < n_elements_per_edge; ++i) {   // each face is shared by two, and only two, adjacent elements
+	int element_id = (*it).adjacent_elements()[i];
+	if (element_id >= 0) {   // not a boundary face
+	  // search point opposite to this face (the j-th node which is not a node of this face)
+	  int j = 0;
+	  for (; j < n_vertices; ++j) {
+	    bool found = false;
+	    for (int k = 0; k < n_vertices_per_edge; ++k) {
+	      if ((*it).node_ids()[k] == elements_(element_id, j)) { found = true; }
+	    }
+	    if (!found) break;
+	  }
+	  neighbors_(element_id, j) = (*it).adjacent_elements()[(i + 1) % 2];
+	}
+      }
+    }
+
+    // for 3D need to compute faces
+
+  } else { // linear network specialization
+    std::map<int, std::vector<int>> adjacent; // for each node, the elements insisting on it
+
+    for (std::size_t i = 0; i < n_elements_; ++i) {
+      adjacent[elements_(i,0)].push_back(i);
+      adjacent[elements_(i,1)].push_back(i);
+    }
+
+    // recover adjacency matrix
+    std::vector<Eigen::Triplet<int>> triplet_list;
+    for(const auto& node : adjacent) {
+      for(std::size_t i = 0; i < node.second.size(); ++i) {
+	for(std::size_t j = i+1; j < node.second.size(); ++j) {
+	  triplet_list.emplace_back(node.second[i], node.second[j], 1);
+	  triplet_list.emplace_back(node.second[j], node.second[i], 1);  // adjacency relation is symmetric
+	}
+      }
     }
-    return;
+    neighbors_.resize(n_elements_, n_elements_);
+    neighbors_.setFromTriplets(triplet_list.begin(), triplet_list.end());
+    neighbors_.makeCompressed();
+  }
+
+  // precompute elements informations for fast access
+  fill_cache();
+  return;
 }
 
 // fill the cache_ data structure with pointers to element objects

test/src/mesh_test.cpp

@@ -14,23 +14,23 @@
 // You should have received a copy of the GNU General Public License
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-#include <fdaPDE/utils.h>
 #include <fdaPDE/mesh.h>
+#include <fdaPDE/utils.h>
 #include <gtest/gtest.h>   // testing framework
 
 #include <memory>
 #include <random>
-#include <unordered_set>
 #include <set>
+#include <unordered_set>
 #include <vector>
 using fdapde::core::Element;
 using fdapde::core::Mesh;
 
 #include "utils/mesh_loader.h"
 #include "utils/utils.h"
+using fdapde::testing::almost_equal;
 using fdapde::testing::MESH_TYPE_LIST;
 using fdapde::testing::MeshLoader;
-using fdapde::testing::almost_equal;
 
 // test suite for testing both non-manifold meshes (2D/3D) and manifold mesh (2.5D/1.5D)
 template <typename E> struct mesh_test : public ::testing::Test {
@@ -48,15 +48,13 @@ TYPED_TEST(mesh_test, elements_construction) {
         auto e = this->mesh_loader.mesh.element(i);
         // check coordinates stored in element built from Mesh object match raw informations
         int j = 0;
-	auto raw_elements = this->mesh_loader.elements_.row(i);
-	auto raw_points = this->mesh_loader.points_;
+        auto raw_elements = this->mesh_loader.elements_.row(i);
+        auto raw_points = this->mesh_loader.points_;
         for (int k = 0; k < raw_elements.size(); ++k) {
             int nodeID = raw_elements[k];
-	    auto p = raw_points.row(nodeID);
+            auto p = raw_points.row(nodeID);
             SVector<TestFixture::N> ePoint = e.coords()[j];
-            for (std::size_t idx = 0; idx < TestFixture::N; ++idx) {
-                EXPECT_TRUE(almost_equal(p[idx], ePoint[idx]));
-            }
+            for (std::size_t idx = 0; idx < TestFixture::N; ++idx) { EXPECT_TRUE(almost_equal(p[idx], ePoint[idx])); }
             j++;
         }
     }
@@ -73,26 +71,53 @@ TYPED_TEST(mesh_test, edges_construction) {
         std::sort(e.begin(), e.end());   // normalize wrt ordering of edge's nodes
         expected_edge_set.push_back(e);
     }
-    // load mesh edges and compute set difference
+    // load mesh edges and compute mask
     std::vector<std::vector<int>> mesh_edge_set;
     std::vector<bool> edge_mask(expected_edge_set.size(), false);
     for (auto it = this->mesh_loader.mesh.edge_begin(); it != this->mesh_loader.mesh.edge_end(); ++it) {
         std::vector<int> e {};
-	for (std::size_t j = 0; j < K; ++j) { e.push_back((*it).node_ids()[j]); }
+        for (std::size_t j = 0; j < K; ++j) { e.push_back((*it).node_ids()[j]); }
         std::sort(e.begin(), e.end());   // normalize wrt ordering of edge's node
 
-	// find this edge in expected set
-	auto search_it = std::find(expected_edge_set.begin(), expected_edge_set.end(), e);
-	if(search_it != expected_edge_set.end()) {
-	  edge_mask[std::distance(expected_edge_set.begin(), search_it)] = true;
-	}
-    }; 
+        // find this edge in expected set
+        auto search_it = std::find(expected_edge_set.begin(), expected_edge_set.end(), e);
+        if (search_it != expected_edge_set.end()) {
+            edge_mask[std::distance(expected_edge_set.begin(), search_it)] = true;
+        }
+    };
     // check all expected edges are indeed computed
     bool result = true;
-    for(bool b : edge_mask) { result &= b; }
+    for (bool b : edge_mask) { result &= b; }
     EXPECT_TRUE(result == true);
 }
 
+// check neighbors informations are computed correctly (up to a permutation of the elements)
+TYPED_TEST(mesh_test, neighbors_construction) {
+    // same number of elements
+    EXPECT_TRUE(this->mesh_loader.neighbors_.size() == this->mesh_loader.mesh.neighbors().size());
+
+    if constexpr (!fdapde::core::is_linear_network<TestFixture::M, TestFixture::N>::value) {
+        bool result = true;
+	int n_row = this->mesh_loader.neighbors_.rows();
+	int n_col = this->mesh_loader.neighbors_.cols();
+	for(std::size_t i = 0; i < n_row; ++i) {
+	  for(std::size_t j = 0; j < n_col; ++j){
+	    if(this->mesh_loader.neighbors_(i,j) != this->mesh_loader.mesh.neighbors()(i,j)) { result = false; }
+	  }
+	}
+        EXPECT_TRUE(result == true);
+    } else {
+      // for linear networks, neighbors are stored as a sparse adjacency matrix
+      bool result = true;
+      for (int k = 0; k < this->mesh_loader.neighbors_.outerSize(); ++k) {
+	for (SpMatrix<int>::InnerIterator it(this->mesh_loader.neighbors_,k); it; ++it) {
+	  if(it.value() != this->mesh_loader.mesh.neighbors().coeff(it.row(), it.col())) { result = false; }
+	}
+      }
+      EXPECT_TRUE(result == true);
+    }
+}
+
 // performs some checks on the mesh topology, e.g. checks that stated neighbors shares exactly M points
 TYPED_TEST(mesh_test, boundary_checks) {
     // cycle over all mesh elements

test/src/utils/mesh_loader.h

@@ -81,18 +81,17 @@ template <typename E> struct MeshLoader {
 	points_ = double_reader.parse_file<DenseStorage>(points_file);
         // realign indexes to 0, if requested
         elements_ = (int_reader.parse_file<DenseStorage>(elements_file).array() - 1).matrix();
-        edges_ = (int_reader.parse_file<DenseStorage>(edges_file).array() - 1).matrix();
+        edges_ = (int_reader.parse_file<DenseStorage>(edges_file).array() > 0)
+	    .select(int_reader.parse_file<DenseStorage>(edges_file).array()-1, -1).matrix();
         boundary_ = int_reader.parse_file<DenseStorage>(boundary_file);
 	if constexpr (!is_linear_network<M, N>::value)
-            neighbors_ = (int_reader.parse_file<DenseStorage>(neighbors_file).array() - 1).matrix();
+	    neighbors_ = (int_reader.parse_file<DenseStorage>(neighbors_file).array() > 0)
+	      .select(int_reader.parse_file<DenseStorage>(neighbors_file).array()-1, -1).matrix();
         else {
             neighbors_ = int_reader.parse_file<SparseStorage>(neighbors_file);
-            for (int k = 0; k < neighbors_.outerSize(); ++k) {
-                for (SpMatrix<int>::InnerIterator it(neighbors_, k); it; ++it) { it.valueRef() = it.value() - 1; }
-            }
         }
         // initialize mesh
-        mesh = E(points_, elements_, neighbors_, boundary_);
+	mesh = E(points_, elements_, boundary_);
     }
     // load default mesh according to dimensionality
     MeshLoader() : MeshLoader(standard_mesh_selector(E::local_dimension, E::embedding_dimension)) {};