From d3e3e4fad1d53a319c4f8210e32017463c25784d Mon Sep 17 00:00:00 2001 From: AlePalu Date: Sat, 1 Jun 2024 17:37:36 +0200 Subject: [PATCH] better projection API, linear networks restored (unstable), node_patch returns std::vector --- fdaPDE/geometry/linear_network.h | 137 +++++++++++++++++++++++++++++++ fdaPDE/geometry/project.h | 32 ++++---- fdaPDE/geometry/segment.h | 2 +- fdaPDE/geometry/tree_search.h | 6 +- fdaPDE/geometry/triangulation.h | 6 +- fdaPDE/utils/symbols.h | 4 +- test/src/simplex_text.cpp | 119 --------------------------- 7 files changed, 164 insertions(+), 142 deletions(-) create mode 100644 fdaPDE/geometry/linear_network.h delete mode 100644 test/src/simplex_text.cpp diff --git a/fdaPDE/geometry/linear_network.h b/fdaPDE/geometry/linear_network.h new file mode 100644 index 00000000..47fb7664 --- /dev/null +++ b/fdaPDE/geometry/linear_network.h @@ -0,0 +1,137 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __LINEAR_NETWORK_H__ +#define __LINEAR_NETWORK_H__ + +#include "segment.h" +#include "../linear_algebra/binary_matrix.h" + +namespace fdapde { +namespace core { + +// template specialization for 1D meshes (bounded intervals) +template class Triangulation; +template <> class Triangulation<1, 2> { + public: + static constexpr int local_dim = 1; + static constexpr int embed_dim = 2; + static constexpr int n_nodes_per_cell = 2; + static constexpr int n_neighbors_per_cell = Dynamic; + static constexpr bool is_manifold = true; + using VertexType = SVector<2>; + using LocationPolicy = TreeSearch>; + + struct CellType : public Segment> { + private: + using Base = Segment>; + using Base::id_; + using Base::mesh_; + public: + CellType() = default; + CellType(int id, const Triangulation* mesh) : Segment>(id, mesh) { } + + DVector neighbors() const { + const auto& v1 = mesh_->node_to_cells().at(mesh_->cells()(id_, 0)); + const auto& v2 = mesh_->node_to_cells().at(mesh_->cells()(id_, 1)); + DVector result(v1.size() + v2.size()); + int i = 0; + for (; i < v1.size(); ++i) result[i] = v1[i]; + for (; i < v1.size() + v2.size(); ++i) result[i] = v2[i]; + return result; + } + }; + + Triangulation() = default; + Triangulation(const DMatrix& nodes, const DMatrix& cells, const DMatrix& boundary) : + nodes_(nodes), cells_(cells), nodes_markers_(boundary) { + // store number of nodes and number of elements + n_nodes_ = nodes_.rows(); + n_cells_ = cells_.rows(); + // compute mesh limits + range_.row(0) = nodes_.colwise().minCoeff(); + range_.row(1) = nodes_.colwise().maxCoeff(); + neighbors_.resize(n_cells_, n_cells_); + // compute node to cells boundings + for (int i = 0; i < n_cells_; ++i) { + node_to_cells_[cells_(i, 0)].push_back(i); + node_to_cells_[cells_(i, 1)].push_back(i); + } + // recover adjacency matrix + SpMatrix adjoint_neighbors; + std::vector> triplet_list; + for (const auto& [node, edges] : node_to_cells_) { + for (int i = 0; i < edges.size(); ++i) { + for (int j = i + 1; j < edges.size(); ++j) triplet_list.emplace_back(edges[j], edges[i], 1); + } + } + adjoint_neighbors.resize(n_cells_, n_cells_); + adjoint_neighbors.setFromTriplets(triplet_list.begin(), triplet_list.end()); + neighbors_ = adjoint_neighbors.selfadjointView(); // symmetrize neighboring relation + }; + + // getters + CellType cell(int id) const { return CellType(id, this); } + VertexType node(int id) const { return nodes_.row(id); } + bool is_node_on_boundary(int id) const { return nodes_markers_[id]; } + const DMatrix& nodes() const { return nodes_; } + const DMatrix& cells() const { return cells_; } + const SpMatrix& neighbors() const { return neighbors_; } + const BinaryVector& boundary() const { return nodes_markers_; } + int n_cells() const { return n_cells_; } + int n_nodes() const { return n_nodes_; } + SVector<2> range() const { return range_; } + + // iterators support + class cell_iterator : public index_based_iterator { + using Base = index_based_iterator; + using Base::index_; + friend Base; + const Triangulation* mesh_; + cell_iterator& operator()(int i) { + Base::val_ = mesh_->cell(i); + return *this; + } + public: + cell_iterator(int index, const Triangulation* mesh) : Base(index, 0, mesh->n_cells_), mesh_(mesh) { + if (index_ < mesh_->n_cells_) operator()(index_); + } + }; + cell_iterator cells_begin() const { return cell_iterator(0, this); } + cell_iterator cells_end() const { return cell_iterator(n_cells_, this); } + + // point location + DVector locate(const DMatrix& points) const { + if (!location_policy_.has_value()) location_policy_ = LocationPolicy(this); + return location_policy_->locate(points); + } + // the set of cells which have node id as vertex + std::vector node_patch(int id) const { return node_to_cells_.at(id); } + protected: + DMatrix nodes_; // physical coordinates of mesh's vertices + DMatrix cells_ {}; // nodes (as row indexes in nodes_ matrix) composing each cell + SpMatrix neighbors_ {}; // ids of faces adjacent to a given face (-1 if no adjacent face) + BinaryVector nodes_markers_ {}; // j-th element is 1 \iff node j is on boundary + std::unordered_map> node_to_cells_; // for each node, the ids of cells sharing it + SVector<2> range_ {}; // mesh bounding box (min and max coordinates) + int n_nodes_ = 0, n_cells_ = 0; + mutable std::optional location_policy_ {}; +}; + +} // namespace core +} // namespace fdapde + +#endif // __LINEAR_NETWORK_H__ diff --git a/fdaPDE/geometry/project.h b/fdaPDE/geometry/project.h index c9cdb3e6..5cbb5eca 100644 --- a/fdaPDE/geometry/project.h +++ b/fdaPDE/geometry/project.h @@ -23,15 +23,20 @@ namespace fdapde { namespace core { -template struct Project; +template class Projection { + private: + const TriangulationType* mesh_; + mutable std::optional> tree_; + public: + Projection() = default; + explicit Projection(const TriangulationType& mesh) : mesh_(&mesh) { } -template <> struct Project { - template static DMatrix compute(const DMatrix& points, const MeshType& mesh) { + DMatrix operator()(const DMatrix& points, tag_exact) const { DVector best = DVector::Constant(points.rows(), std::numeric_limits::max()); - DMatrix proj(points.rows(), MeshType::embed_dim); - for (typename MeshType::cell_iterator it = mesh.cells_begin(); it != mesh.cells_end(); ++it) { + DMatrix proj(points.rows(), TriangulationType::embed_dim); + for (typename TriangulationType::cell_iterator it = mesh_->cells_begin(); it != mesh_->cells_end(); ++it) { for (int i = 0; i < points.rows(); ++i) { - SVector proj_point = it->nearest(points.row(i)); + SVector proj_point = it->nearest(points.row(i)); double dist = (proj_point - points.row(i).transpose()).norm(); if (dist < best[i]) { best[i] = dist; @@ -41,20 +46,18 @@ template <> struct Project { } return proj; } -}; -template <> struct Project { - template static DMatrix compute(const DMatrix& points, const MeshType& mesh) { - DMatrix proj(points.rows(), MeshType::embed_dim); + DMatrix operator()(const DMatrix& points, tag_not_exact) const { + DMatrix proj(points.rows(), TriangulationType::embed_dim); // build kdtree of mesh nodes for fast nearest neighborhood searches - KDTree tree_(mesh.nodes()); + if (!tree_.has_value()) tree_ = KDTree(mesh_->nodes()); for (int i = 0; i < points.rows(); ++i) { // find nearest mesh node (in euclidean sense, approximation) - typename KDTree::iterator it = tree_.nn_search(points.row(i)); + typename KDTree::iterator it = tree_->nn_search(points.row(i)); // search nearest element in the node patch double best = std::numeric_limits::max(); - for (int j : mesh.node_patch(*it)) { - SVector proj_point = mesh.cell(j).nearest(points.row(i)); + for (int j : mesh_->node_patch(*it)) { + SVector proj_point = mesh_->cell(j).nearest(points.row(i)); double dist = (proj_point - points.row(i).transpose()).norm(); if (dist < best) { best = dist; @@ -64,6 +67,7 @@ template <> struct Project { } return proj; } + DMatrix operator()(const DMatrix& points) const { return operator()(points, fdapde::NotExact); } }; } // namespace core diff --git a/fdaPDE/geometry/segment.h b/fdaPDE/geometry/segment.h index 309525af..f4042748 100644 --- a/fdaPDE/geometry/segment.h +++ b/fdaPDE/geometry/segment.h @@ -44,7 +44,7 @@ template class Segment : public Simplex node_ids() const { return mesh_->cells().row(id_); } bool on_boundary() const { return boundary_; } operator bool() const { return mesh_ != nullptr; } - private: + protected: int id_ = 0; // segment ID in the physical mesh const Triangulation* mesh_ = nullptr; bool boundary_ = false; // true if the element has at least one vertex on the boundary diff --git a/fdaPDE/geometry/tree_search.h b/fdaPDE/geometry/tree_search.h index 9dae63c7..f3816494 100644 --- a/fdaPDE/geometry/tree_search.h +++ b/fdaPDE/geometry/tree_search.h @@ -59,11 +59,11 @@ template class TreeSearch { tree_ = KDTree<2 * embed_dim>(std::move(data)); // organize elements in a KD-tree structure } // finds all the elements containing p - std::unordered_set all_locate(const SVector& p) const { - std::unordered_set result; + std::vector all_locate(const SVector& p) const { + std::vector result; for (int id : tree_.range_search(query(p))) { typename MeshType::CellType c = mesh_->cell(id); - if (c.contains(p)) { result.insert(c.id()); } + if (c.contains(p)) { result.push_back(c.id()); } } return result; } diff --git a/fdaPDE/geometry/triangulation.h b/fdaPDE/geometry/triangulation.h index 36f0fd26..df25fc79 100644 --- a/fdaPDE/geometry/triangulation.h +++ b/fdaPDE/geometry/triangulation.h @@ -253,8 +253,8 @@ template class Triangulation<2, N> : public TriangulationBase<2, N, Tria if (!location_policy_.has_value()) location_policy_ = LocationPolicy(this); return location_policy_->locate(points); } - // computes the set of elements which have node id as vertex - std::unordered_set node_patch(int id) const { + // the set of cells which have node id as vertex + std::vector node_patch(int id) const { if (!location_policy_.has_value()) location_policy_ = LocationPolicy(this); return location_policy_->all_locate(Base::node(id)); } @@ -466,7 +466,7 @@ template <> class Triangulation<3, 3> : public TriangulationBase<3, 3, Triangula return location_policy_->locate(points); } // computes the set of elements which have node id as vertex - std::unordered_set node_patch(int id) const { + std::vector node_patch(int id) const { if (!location_policy_.has_value()) location_policy_ = LocationPolicy(this); return location_policy_->all_locate(Base::node(id)); } diff --git a/fdaPDE/utils/symbols.h b/fdaPDE/utils/symbols.h index 55e8125c..c40165ed 100644 --- a/fdaPDE/utils/symbols.h +++ b/fdaPDE/utils/symbols.h @@ -41,8 +41,8 @@ constexpr int Dynamic = -1; // used when the size of a vector or matrix is constexpr int random_seed = -1; // signals that a random seed is used somewhere // algorithm computation policies -struct Exact { }; -struct NotExact { }; +static struct tag_exact { } Exact; +static struct tag_not_exact { } NotExact; template struct static_dynamic_vector_selector { using type = typename std::conditional, SVector>::type; diff --git a/test/src/simplex_text.cpp b/test/src/simplex_text.cpp deleted file mode 100644 index 9089d025..00000000 --- a/test/src/simplex_text.cpp +++ /dev/null @@ -1,119 +0,0 @@ -// This file is part of fdaPDE, a C++ library for physics-informed -// spatial and functional data analysis. -// -// This program is free software: you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation, either version 3 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see . - -#include // testing framework -#include - -#include -#include -using fdapde::core::Simplex; - -#include "utils/utils.h" -using fdapde::testing::almost_equal; - -TEST(simplex_test, triangle_2d) { - SMatrix<2, 3> coords; - coords.col(0) = SVector<2>(0.0, 0.0); - coords.col(1) = SVector<2>(0.5, 0.0); - coords.col(2) = SVector<2>(0.0, 0.8); - Simplex<2, 2> t(coords); - - EXPECT_TRUE(almost_equal(t.measure(), (0.5 * 0.8) / 2)); - EXPECT_TRUE(almost_equal(t.circumcenter(), SVector<2>(0.25, 0.4))); - EXPECT_TRUE(almost_equal(t.circumradius(), std::sqrt(0.25 * 0.25 + 0.4 * 0.4))); - EXPECT_TRUE(almost_equal(t.barycenter(), SVector<2>(0.166666666666667, 0.2666666666666667))); - EXPECT_TRUE(t.bounding_box().first == SVector<2>(0, 0) && t.bounding_box().second == SVector<2>(0.5, 0.8)); - EXPECT_TRUE(t.contains(SVector<2>(0.25, 1.25)) == 0); // point falls outside - EXPECT_TRUE(t.contains(SVector<2>(0.25, 0.25)) == 1); // point falls inside - EXPECT_TRUE(t.contains(SVector<2>(0.25, 0.00)) == 2); // point falls on edge - EXPECT_TRUE(t.contains(SVector<2>(0.00, 0.00)) == 3); // point falls on vertex - - // compute perimeter of triangle by face iterator - EXPECT_TRUE(almost_equal( - std::accumulate( - t.face_begin(), t.face_end(), 0.0, [](double perimeter, auto& f) { return perimeter + f.measure(); }), - 0.5 + 0.8 + std::sqrt(0.5 * 0.5 + 0.8 * 0.8))); - // extract face - Simplex<2, 2>::FaceType f = t.face(2); - EXPECT_TRUE(f[0] == SVector<2>(0.5, 0) && f[1] == SVector<2>(0, 0.8)); - EXPECT_TRUE(almost_equal(f.measure(), std::sqrt(0.5 * 0.5 + 0.8 * 0.8))); - EXPECT_TRUE(almost_equal(f.circumcenter(), f.barycenter())); -} - -TEST(simplex_test, triangle_3d) { - SMatrix<3, 3> coords; - coords.col(0) = SVector<3>(0.0, 0.0, 0.0); - coords.col(1) = SVector<3>(0.5, 0.2, 0.0); - coords.col(2) = SVector<3>(0.0, 0.8, 0.6); - Simplex<2, 3> t(coords); - - EXPECT_TRUE(almost_equal(t.measure(), 0.25709920264364883)); - EXPECT_TRUE(almost_equal( - t.circumcenter(), t[0] + ((t[2] - t[0]).squaredNorm() * ((t[1] - t[0]).cross(t[2] - t[0])).cross(t[1] - t[0]) + - (t[1] - t[0]).squaredNorm() * ((t[2] - t[0]).cross(t[1] - t[0])).cross(t[2] - t[0])) / - (2 * ((t[1] - t[0]).cross(t[2] - t[0])).squaredNorm()))); - EXPECT_TRUE(almost_equal(t.circumradius(), t.circumcenter().norm())); // distance from zero - EXPECT_TRUE(almost_equal(t.barycenter(), SVector<3>(0.1666666666666667, 0.3333333333333333, 0.2000000000000000))); - EXPECT_TRUE(t.bounding_box().first == SVector<3>(0, 0, 0) && t.bounding_box().second == SVector<3>(0.5, 0.8, 0.6)); - EXPECT_TRUE(t.contains(SVector<3>(0.00, 1.00, 1.00)) == 0); // point falls outside - EXPECT_TRUE(t.contains(SVector<3>(0.25, 0.30, 0.15)) == 1); // point falls inside - EXPECT_TRUE(t.contains(SVector<3>(0.00, 0.40, 0.30)) == 2); // point falls on edge - EXPECT_TRUE(t.contains(SVector<3>(0.50, 0.20, 0.00)) == 3); // point falls on vertex - EXPECT_TRUE(almost_equal(t.normal(), ((t[1] - t[0]).cross(t[2] - t[0])).normalized())); - - // compute perimeter of triangle by face iterator - EXPECT_TRUE(almost_equal( - std::accumulate( - t.face_begin(), t.face_end(), 0.0, [](double perimeter, auto& f) { return perimeter + f.measure(); }), - 2.523402260893061)); - // extract face - Simplex<2, 3>::FaceType f = t.face(0); - EXPECT_TRUE(f[0] == SVector<3>(0.0, 0.0, 0.0) && f[1] == SVector<3>(0.5, 0.2, 0.0)); - EXPECT_TRUE(almost_equal(f.measure(), std::sqrt(0.5 * 0.5 + 0.2 * 0.2))); - EXPECT_TRUE(almost_equal(f.circumcenter(), f.barycenter())); -} - -TEST(simplex_test, tetrahedron_3d) { - SMatrix<3, 4> coords; - coords.col(0) = SVector<3>(0.0, 0.0, 0.0); - coords.col(1) = SVector<3>(0.4, 0.2, 0.0); - coords.col(2) = SVector<3>(0.0, 0.8, 0.6); - coords.col(3) = SVector<3>(0.4, 0.6, 0.8); - Simplex<3, 3> t(coords); - - EXPECT_TRUE(almost_equal(t.measure(), 0.0266666666666666)); - EXPECT_TRUE(almost_equal(t.circumcenter(), SVector<3>(0.11, 0.28, 0.46))); - EXPECT_TRUE(almost_equal(t.circumradius(), std::sqrt(0.11 * 0.11 + 0.28 * 0.28 + 0.46 * 0.46))); - EXPECT_TRUE(almost_equal(t.barycenter(), SVector<3>(0.2, 0.4, 0.35))); - EXPECT_TRUE(t.bounding_box().first == SVector<3>(0, 0, 0) && t.bounding_box().second == SVector<3>(0.4, 0.8, 0.8)); - EXPECT_TRUE(t.contains(SVector<3>(0.25, 1.25, 0.25)) == 0); // point falls outside - EXPECT_TRUE(t.contains(t.barycenter()) == 1); // point falls inside - EXPECT_TRUE(t.contains(t.face(0).barycenter()) == 2); // point falls on face - EXPECT_TRUE(t.contains(SVector<3>(0.40, 0.60, 0.80)) == 3); // point falls on vertex - - // compute surface area of tetrahedron by face iterator - EXPECT_TRUE(almost_equal( - std::accumulate( - t.face_begin(), t.face_end(), 0.0, [](double perimeter, auto& f) { return perimeter + f.measure(); }), - 0.86430301)); - - // extract face - Simplex<3, 3>::FaceType f = t.face(1); - EXPECT_TRUE( - f[0] == SVector<3>(0.0, 0.0, 0.0) && f[1] == SVector<3>(0.4, 0.2, 0.0) && f[2] == SVector<3>(0.4, 0.6, 0.8)); - EXPECT_TRUE(almost_equal(f.measure(), 0.19595918)); - EXPECT_TRUE(f.normal() == ((f[1] - f[0]).cross(f[2] - f[0])).normalized()); -}