diff --git a/fdaPDE/mesh/element.h b/fdaPDE/mesh/element.h index a1605a4c..c2dc06b5 100644 --- a/fdaPDE/mesh/element.h +++ b/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. diff --git a/fdaPDE/mesh/mesh.h b/fdaPDE/mesh/mesh.h index 383fe188..1270bfdb 100644 --- a/fdaPDE/mesh/mesh.h +++ b/fdaPDE/mesh/mesh.h @@ -51,7 +51,7 @@ template struct neighboring_structure { // access to domain's triangulation, M: tangent space dimension, N: embedding space dimension template class Mesh { private: - // coordinates of points costituting the vertices of mesh elements + // coordinates of points making the vertices of mesh elements DMatrix nodes_ {}; int n_nodes_ = 0; // identifiers of points (as row indexes in points_ matrix) composing each element, by row @@ -65,14 +65,13 @@ template class Mesh { std::vector edges_ {}; std::vector 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> cache_ {}; void fill_cache(); public: Mesh() = default; - Mesh(const DMatrix& nodes, const DMatrix& elements, - const typename neighboring_structure::type& neighbors, const DMatrix& boundary); + Mesh(const DMatrix& nodes, const DMatrix& elements, const DMatrix& boundary); // getters const Element& element(int ID) const { return cache_[ID]; } @@ -161,59 +160,103 @@ template 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 -Mesh::Mesh( - const DMatrix& nodes, const DMatrix& elements, - const typename neighboring_structure::type& neighbors, const DMatrix& 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::Mesh(const DMatrix& nodes, const DMatrix& elements, const DMatrix& 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, in a row major format auto edge_pattern = combinations(); std::unordered_map, int, std_array_hash> visited; - std::array 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::value) { + neighbors_ = DMatrix::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> 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> 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 diff --git a/test/src/mesh_test.cpp b/test/src/mesh_test.cpp index a30e942d..dc634626 100644 --- a/test/src/mesh_test.cpp +++ b/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 . -#include #include +#include #include // testing framework #include #include -#include #include +#include #include 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 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 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> mesh_edge_set; std::vector 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 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::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::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 diff --git a/test/src/utils/mesh_loader.h b/test/src/utils/mesh_loader.h index 1df2a767..f9ce2b93 100644 --- a/test/src/utils/mesh_loader.h +++ b/test/src/utils/mesh_loader.h @@ -81,18 +81,17 @@ template struct MeshLoader { points_ = double_reader.parse_file(points_file); // realign indexes to 0, if requested elements_ = (int_reader.parse_file(elements_file).array() - 1).matrix(); - edges_ = (int_reader.parse_file(edges_file).array() - 1).matrix(); + edges_ = (int_reader.parse_file(edges_file).array() > 0) + .select(int_reader.parse_file(edges_file).array()-1, -1).matrix(); boundary_ = int_reader.parse_file(boundary_file); if constexpr (!is_linear_network::value) - neighbors_ = (int_reader.parse_file(neighbors_file).array() - 1).matrix(); + neighbors_ = (int_reader.parse_file(neighbors_file).array() > 0) + .select(int_reader.parse_file(neighbors_file).array()-1, -1).matrix(); else { neighbors_ = int_reader.parse_file(neighbors_file); - for (int k = 0; k < neighbors_.outerSize(); ++k) { - for (SpMatrix::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)) {};