Skip to content

Commit 7956032

Browse files
AlePaluAlePalu
committed
expressions node which might return a temporary always casts result to Scalar, adjusting triangular matrix type system, minor changes
1 parent f9e1cc1 commit 7956032

File tree

20 files changed

+579
-431
lines changed

20 files changed

+579
-431
lines changed

‎fdaPDE/linear_algebra.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,17 @@ template <typename XprType> struct MatrixCoeffWiseExpr;
3939
[[maybe_unused]] static constexpr int LhsMode = 0;
4040
[[maybe_unused]] static constexpr int RhsMode = 1;
4141

42+
namespace internals {
43+
44+
// tag types to enable/disable costly data integrity checks
45+
struct checked_t { };
46+
struct unchecked_t { };
47+
48+
} // namespace internals
49+
50+
[[maybe_unused]] inline constexpr internals::checked_t checked {};
51+
[[maybe_unused]] inline constexpr internals::unchecked_t unchecked {};
52+
4253
} // namespace fdapde
4354

4455
#include "src/linear_algebra/traits.h"

‎fdaPDE/src/geometry/polygon.h‎

Lines changed: 28 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -29,13 +29,13 @@ template <int LocalDim, int EmbedDim> class Polygon {
2929

3030
// constructors
3131
Polygon() noexcept = default;
32-
Polygon(const Eigen::Matrix<double, Dynamic, Dynamic>& nodes) noexcept : triangulation_() {
32+
Polygon(const Matrix<double, Dynamic, Dynamic>& nodes) noexcept : triangulation_() {
3333
fdapde_assert(nodes.rows() > 0 && nodes.cols() == embed_dim);
3434
if (internals::are_2d_counterclockwise_sorted(nodes)) {
3535
triangulate_(nodes);
3636
} else { // nodes are in clocwise order, reverse node ordering
3737
int n_nodes = nodes.rows();
38-
Eigen::Matrix<double, Dynamic, Dynamic> reversed_nodes(n_nodes, embed_dim);
38+
Matrix<double, Dynamic, Dynamic> reversed_nodes(n_nodes, embed_dim);
3939
for (int i = 0; i < n_nodes; ++i) { reversed_nodes.row(i) = nodes.row(n_nodes - 1 - i); }
4040
triangulate_(reversed_nodes);
4141
}
@@ -44,11 +44,9 @@ template <int LocalDim, int EmbedDim> class Polygon {
4444
Polygon(const Polygon&) noexcept = default;
4545
Polygon(Polygon&&) noexcept = default;
4646
// observers
47-
const Eigen::Matrix<double, Dynamic, Dynamic>& nodes() const { return triangulation_.nodes(); }
48-
const Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>& cells() const { return triangulation_.cells(); }
49-
Eigen::Map<const Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>> edges() const {
50-
return triangulation_.edges();
51-
}
47+
const Matrix<double, Dynamic, Dynamic>& nodes() const { return triangulation_.nodes(); }
48+
const Matrix<int, Dynamic, Dynamic>& cells() const { return triangulation_.cells(); }
49+
MatrixView<const int, Dynamic, Dynamic> edges() const { return triangulation_.edges(); }
5250
const Triangulation<local_dim, embed_dim>& triangulation() const { return triangulation_; }
5351
double measure() const { return triangulation_.measure(); }
5452
int n_nodes() const { return triangulation_.n_nodes(); }
@@ -57,35 +55,37 @@ template <int LocalDim, int EmbedDim> class Polygon {
5755
// test whether point p is contained in polygon
5856
template <int Rows, int Cols>
5957
requires((Rows == embed_dim && Cols == 1) || (Cols == embed_dim && Rows == 1))
60-
bool contains(const Eigen::Matrix<double, Rows, Cols>& p) const {
58+
bool contains(const Matrix<double, Rows, Cols>& p) const {
6159
return triangulation_.locate(p) != -1;
6260
}
6361
// random sample points in polygon
64-
Eigen::Matrix<double, Dynamic, Dynamic> sample(int n_samples, int seed = random_seed) const {
62+
Matrix<double, Dynamic, Dynamic> sample(int n_samples, int seed = random_seed) const {
6563
return triangulation_.sample(n_samples, seed);
6664
}
6765
private:
6866
// perform polygon triangulation
69-
void triangulate_(const Eigen::Matrix<double, Dynamic, Dynamic>& nodes) {
67+
void triangulate_(const Matrix<double, Dynamic, Dynamic>& nodes) {
7068
std::vector<int> cells;
7169
// perform monotone partitioning
7270
std::vector<std::vector<int>> poly_partition = monotone_partition_(nodes);
7371
// triangulate each monotone polygon
7472
for (const std::vector<int>& poly : poly_partition) {
75-
std::vector<int> local_cells = triangulate_monotone_(nodes(poly, Eigen::all));
73+
Matrix<double, Dynamic, Dynamic> local_nodes(poly.size(), nodes.cols());
74+
for (int i = 0, n = poly.size(); i < n; ++i) { local_nodes.row(i) = nodes.row(poly[i]); }
75+
std::vector<int> local_cells = triangulate_monotone_(local_nodes);
7676
// move local node numbering to global node numbering
7777
for (std::size_t i = 0; i < local_cells.size(); ++i) { local_cells[i] = poly[local_cells[i]]; }
7878
cells.insert(cells.end(), local_cells.begin(), local_cells.end());
7979
}
8080
// set-up face-based data structure
8181
triangulation_ = Triangulation<LocalDim, EmbedDim>(
82-
nodes, Eigen::Map<Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>>(cells.data(), cells.size() / 3, 3),
83-
Eigen::Matrix<int, Dynamic, 1>::Ones(nodes.rows()));
82+
nodes, MatrixView<int, Dynamic, Dynamic>(cells.data(), cells.size() / 3, 3),
83+
Matrix<int, Dynamic, 1>::Ones(nodes.rows()));
8484
}
8585

8686
// partition an arbitrary polygon P into a set of monotone polygons (plane sweep approach, section 3.2 of De Berg,
8787
// M. (2000). Computational geometry: algorithms and applications. Springer Science & Business Media.)
88-
std::vector<std::vector<int>> monotone_partition_(const Eigen::Matrix<double, Dynamic, Dynamic>& coords) {
88+
std::vector<std::vector<int>> monotone_partition_(const Matrix<double, Dynamic, Dynamic>& coords) {
8989
using poly_t = DCEL<local_dim, embed_dim>;
9090
using halfedge_t = typename poly_t::halfedge_t;
9191
using halfedge_ptr_t = std::add_pointer_t<halfedge_t>;
@@ -253,7 +253,7 @@ template <int LocalDim, int EmbedDim> class Polygon {
253253
return monotone_partition;
254254
}
255255
// triangulate monotone polygon (returns a RowMajor ordered matrix of cells).
256-
std::vector<int> triangulate_monotone_(const Eigen::Matrix<double, Dynamic, Dynamic>& nodes) {
256+
std::vector<int> triangulate_monotone_(const Matrix<double, Dynamic, Dynamic>& nodes) {
257257
// every triangulation of a polygon of n points has n - 2 triangles (lemma 1.2.2 of (1))
258258
std::vector<int> cells;
259259
int n_nodes = nodes.rows();
@@ -361,15 +361,15 @@ template <int LocalDim, int EmbedDim> class MultiPolygon {
361361

362362
MultiPolygon() : n_polygons_(0) { }
363363
// rings is a vector of matrix of coordinates, where, each inner matrix defines a closed non self-intersecting loop
364-
MultiPolygon(const std::vector<Eigen::Matrix<double, Dynamic, Dynamic>>& rings) : n_polygons_(0) {
364+
MultiPolygon(const std::vector<Matrix<double, Dynamic, Dynamic>>& rings) : n_polygons_(0) {
365365
// ESRI shapefile format specification: loops in clockwise order define the outer border of a polygon,
366366
// loops in counterclockwise order defines a hole inside the last found outer polygonal ring
367367
int n_rings = rings.size();
368368
if (n_rings == 1) {
369369
triangulation_ = Polygon<local_dim, embed_dim>(rings[0]).triangulation();
370370
} else {
371371
// detect polygons with holes
372-
using iterator = typename std::vector<Eigen::Matrix<double, Dynamic, Dynamic>>::const_iterator;
372+
using iterator = typename std::vector<Matrix<double, Dynamic, Dynamic>>::const_iterator;
373373
constexpr int n_nodes_per_cell = 3;
374374
std::vector<std::pair<iterator, iterator>> polygons;
375375
iterator begin = rings.begin();
@@ -385,10 +385,10 @@ template <int LocalDim, int EmbedDim> class MultiPolygon {
385385
polygons.emplace_back(begin, end);
386386
begin = end;
387387
}
388-
Eigen::Matrix<double, Dynamic, Dynamic> nodes(n_nodes, embed_dim);
388+
Matrix<double, Dynamic, Dynamic> nodes(n_nodes, embed_dim);
389389
std::vector<int> cells;
390390
int nodes_off_ = 0;
391-
for (const Eigen::Matrix<double, Dynamic, Dynamic>& coords : rings) {
391+
for (const Matrix<double, Dynamic, Dynamic>& coords : rings) {
392392
// triangulate polygon
393393
Triangulation<local_dim, embed_dim> poly_tri = Polygon<local_dim, embed_dim>(coords).triangulation();
394394
nodes.middleRows(nodes_off_, poly_tri.n_nodes()) = poly_tri.nodes();
@@ -402,22 +402,18 @@ template <int LocalDim, int EmbedDim> class MultiPolygon {
402402
}
403403
// set up face-based storage
404404
triangulation_ = Triangulation<local_dim, embed_dim>(
405-
nodes,
406-
Eigen::Map<Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>>(
407-
cells.data(), cells.size() / n_nodes_per_cell, n_nodes_per_cell),
408-
Eigen::Matrix<int, Dynamic, 1>::Ones(nodes.rows()));
405+
nodes, MatrixView<int, Dynamic, Dynamic>(cells.data(), cells.size() / n_nodes_per_cell, n_nodes_per_cell),
406+
Matrix<int, Dynamic, 1>::Ones(nodes.rows()));
409407
}
410408
}
411409
// observers
412-
const Eigen::Matrix<double, Dynamic, Dynamic>& nodes() const { return triangulation_.nodes(); }
413-
const Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>& cells() const { return triangulation_.cells(); }
414-
Eigen::Map<const Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>> edges() const {
415-
return triangulation_.edges();
416-
}
410+
const Matrix<double, Dynamic, Dynamic>& nodes() const { return triangulation_.nodes(); }
411+
const Matrix<int, Dynamic, Dynamic>& cells() const { return triangulation_.cells(); }
412+
MatrixView<const int, Dynamic, Dynamic> edges() const { return triangulation_.edges(); }
417413
// computes only the boundary edges of the multipoligon (discards triangulation's edges)
418-
Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor> boundary_edges() const {
414+
Matrix<int, Dynamic, Dynamic> boundary_edges() const {
419415
int n_edges = triangulation_.n_boundary_edges();
420-
Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor> m(n_edges, 2);
416+
Matrix<int, Dynamic, Dynamic> m(n_edges, 2);
421417
int j = 0;
422418
for (int i = 0, k = triangulation_.n_edges(); i < k; ++i) {
423419
if (triangulation_.is_edge_on_boundary(i)) { m.row(j++) = triangulation_.edges().row(i); }
@@ -434,11 +430,11 @@ template <int LocalDim, int EmbedDim> class MultiPolygon {
434430
// test whether point p is contained in polygon
435431
template <int Rows, int Cols>
436432
requires((Rows == embed_dim && Cols == 1) || (Cols == embed_dim && Rows == 1))
437-
bool contains(const Eigen::Matrix<double, Rows, Cols>& p) const {
433+
bool contains(const Matrix<double, Rows, Cols>& p) const {
438434
return triangulation_.locate(p) != -1;
439435
}
440436
// random sample points in polygon
441-
Eigen::Matrix<double, Dynamic, Dynamic> sample(int n_samples, int seed = random_seed) const {
437+
Matrix<double, Dynamic, Dynamic> sample(int n_samples, int seed = random_seed) const {
442438
return triangulation_.sample(n_samples, seed);
443439
}
444440
private:

‎fdaPDE/src/geometry/r_tree.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -379,8 +379,8 @@ class RTree {
379379
item_t() = default;
380380
template <typename BBoxT>
381381
requires(std::is_convertible_v<BBoxT, bbox_t>)
382-
explicit item_t(int index, const BBoxT& bbox) : data_(index), type_(type_t::DATA), bbox_(bbox) { }
383-
explicit item_t(node_t* node) : node_(node), type_(type_t::NODE), bbox_(node_->bbox()) { }
382+
explicit item_t(int index, const BBoxT& bbox) : type_(type_t::DATA), data_(index), bbox_(bbox) { }
383+
explicit item_t(node_t* node) : type_(type_t::NODE), node_(node), bbox_(node_->bbox()) { }
384384
// copy/move semantic
385385
item_t(const item_t& other) = default;
386386
item_t(item_t&& other) = default;

‎fdaPDE/src/geometry/simplex.h‎

Lines changed: 10 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -53,8 +53,8 @@ template <int Order_, int EmbedDim_> class Simplex {
5353
double measure() const { return measure_; }
5454
// simplex minimal enclosing rectangle
5555
std::array<double, 2 * embed_dim> bbox() const {
56-
NodeType ll = coords_.rowwise().minCoeff();
57-
NodeType ur = coords_.rowwise().maxCoeff();
56+
NodeType ll = coords_.rowwise().min();
57+
NodeType ur = coords_.rowwise().max();
5858
std::array<double, 2 * embed_dim> bbox_;
5959
for (int i = 0; i < embed_dim; ++i) {
6060
bbox_[i] = ll[i];
@@ -69,8 +69,8 @@ template <int Order_, int EmbedDim_> class Simplex {
6969
// writes the point p in the barycentric coordinate system of this simplex
7070
Matrix<double, local_dim + 1, 1> barycentric_coords(const NodeType& p) const {
7171
Matrix<double, local_dim + 1, 1> z;
72-
z.bottomRows(local_dim) = invJ_ * (p - coords_.col(0));
73-
z[0] = 1 - z.bottomRows(local_dim).sum();
72+
z.bottom_rows(local_dim) = invJ_ * (p - coords_.col(0));
73+
z[0] = 1 - z.bottom_rows(local_dim).sum();
7474
return z;
7575
}
7676

@@ -130,10 +130,10 @@ template <int Order_, int EmbedDim_> class Simplex {
130130
}
131131
// move x to barycentric coordinates
132132
Matrix<double, local_dim + 1, 1> z;
133-
z.bottomRows(local_dim) = invJ_ * (x - coords_.col(0));
134-
z[0] = 1 - z.bottomRows(local_dim).sum();
135-
if ((z.array() < -machine_epsilon).any()) return ContainsReturnType::OUTSIDE;
136-
int nonzeros = (z.array() > machine_epsilon).count();
133+
z.bottom_rows(local_dim) = invJ_ * (x - coords_.col(0));
134+
z[0] = 1 - z.bottom_rows(local_dim).sum();
135+
if ((z.cwise() < -machine_epsilon).any()) return ContainsReturnType::OUTSIDE;
136+
int nonzeros = (z.cwise() > machine_epsilon).count();
137137
if (nonzeros == 1) return ContainsReturnType::ON_VERTEX;
138138
if (nonzeros == n_nodes_per_face) return ContainsReturnType::ON_FACE;
139139
return ContainsReturnType::INSIDE;
@@ -170,11 +170,9 @@ template <int Order_, int EmbedDim_> class Simplex {
170170
Matrix<double, local_dim + 1, 1> q = barycentric_coords(p);
171171
// check if point inside simplex
172172
if constexpr (local_dim != embed_dim) {
173-
if (
174-
(q.array() > -machine_epsilon).all() && supporting_plane().distance(p) < machine_epsilon)
175-
return p;
173+
if ((q.cwise() > -machine_epsilon).all() && supporting_plane().distance(p) < machine_epsilon) return p;
176174
} else {
177-
if ((q.array() > -machine_epsilon).all()) return p;
175+
if ((q.cwise() > -machine_epsilon).all()) return p;
178176
}
179177
if constexpr (Order_ == 1) { // end of recursion
180178
if (q[0] < 0) return coords_.col(1);

‎fdaPDE/src/geometry/triangulation.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -595,7 +595,7 @@ template <int N> class Triangulation<2, N> : public TriangulationBase<2, N, Tria
595595
for (int i = 0, n = pts.rows(); i < n; ++i) {
596596
std::vector<int> candidates = spatial_index_->locate_query(pts.row(i));
597597
for (int j : candidates) { // perform exact containment test
598-
if (cell(j).contains(pts.row(i))) cells[i] = j;
598+
if (cell(j).contains(pts.row(i))) { cells[i] = j; }
599599
}
600600
}
601601
return cells;

‎fdaPDE/src/geometry/walk.h‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -33,17 +33,17 @@ template <typename MeshType> class WalkSearch {
3333
static_assert(MeshType::local_dim == MeshType::embed_dim);
3434
};
3535
// finds element containing p, returns nullptr if element not found
36-
int locate(const Eigen::Matrix<double, embed_dim, 1>& p) const {
36+
int locate(const Matrix<double, embed_dim, 1>& p) const {
3737
// start search from random element
3838
std::random_device rng {};
3939
std::uniform_int_distribution<std::size_t> uniform_int(0, mesh_->n_cells() - 1);
40-
std::size_t next = uniform_int(rng);
40+
std::size_t next = uniform_int(rng);
4141

4242
std::unordered_set<std::size_t> visited_;
4343
while (!mesh_->cell(next).contains(p) || visited_.find(next) != visited_.end()) {
4444
visited_.insert(next);
4545
// compute barycantric coordinates
46-
Eigen::Matrix<double, embed_dim + 1, 1> bary_coord = mesh_->cell(next).barycentric_coords(p);
46+
Matrix<double, embed_dim + 1, 1> bary_coord = mesh_->cell(next).barycentric_coords(p);
4747
// find minimum baricentric coordinate and move to element insisting of opposite face
4848
std::size_t min_bary_coord_index;
4949
bary_coord.minCoeff(&min_bary_coord_index);
@@ -54,10 +54,10 @@ template <typename MeshType> class WalkSearch {
5454
}
5555
template <typename CoordsMatrix>
5656
requires(internals::is_eigen_dense_xpr_v<CoordsMatrix>)
57-
Eigen::Matrix<int, Dynamic, 1> locate(const CoordsMatrix& locs) const {
57+
Matrix<int, Dynamic, 1> locate(const CoordsMatrix& locs) const {
5858
fdapde_assert(locs.cols() == embed_dim);
59-
Eigen::Matrix<int, Dynamic, 1> ids(locs.rows());
60-
for (int i = 0; i < locs.rows(); ++i) { ids[i] = locate(Eigen::Matrix<double, embed_dim, 1>(locs.row(i))); }
59+
Matrix<int, Dynamic, 1> ids(locs.rows());
60+
for (int i = 0; i < locs.rows(); ++i) { ids[i] = locate(Matrix<double, embed_dim, 1>(locs.row(i))); }
6161
return ids;
6262
}
6363
};

0 commit comments

Comments
 (0)