Skip to content

Commit 9e56bed

Browse files
AlePaluAlePalu
committed
coeffwise operations, matrix-scalar multiplication made explicit, vector_assignment_executor (can assign vector expressions of different shapes), iterator support for blocks, assignment from std::initializer_lists, MatrixExpr conversion to Matrix<bool>, various bug fixing. Started Eigen removal from geometry module
1 parent afe2349 commit 9e56bed

34 files changed

+1044
-641
lines changed

‎fdaPDE/geometry.h‎

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -20,29 +20,29 @@
2020
// clang-format off
2121

2222
// include required modules
23-
#include "linear_algebra.h" // pull Eigen first
23+
#include "linear_algebra.h"
2424
#include "utility.h"
25-
#include "io.h"
25+
/* #include "io.h" */
2626

2727
// minimal geometric entites
2828
#include "src/geometry/utility.h"
29-
#include "src/geometry/primitives.h"
29+
/* #include "src/geometry/primitives.h" */
3030
#include "src/geometry/hyperplane.h"
3131
#include "src/geometry/simplex.h"
3232
#include "src/geometry/segment.h"
3333
#include "src/geometry/triangle.h"
3434
#include "src/geometry/tetrahedron.h"
3535
// algorithms
36-
#include "src/geometry/kd_tree.h"
37-
#include "src/geometry/tree_search.h"
38-
#include "src/geometry/walk_search.h"
39-
#include "src/geometry/projection.h"
36+
/* #include "src/geometry/kd_tree.h" */
37+
/* #include "src/geometry/tree_search.h" */
38+
/* #include "src/geometry/walk_search.h" */
39+
/* #include "src/geometry/projection.h" */
4040
// data structures
4141
#include "src/geometry/triangulation.h"
42-
#include "src/geometry/interval.h"
43-
#include "src/geometry/linear_network.h"
44-
#include "src/geometry/dcel.h"
45-
#include "src/geometry/polygon.h"
42+
/* #include "src/geometry/interval.h" */
43+
/* #include "src/geometry/linear_network.h" */
44+
/* #include "src/geometry/dcel.h" */
45+
/* #include "src/geometry/polygon.h" */
4646

4747
// clang-format on
4848

‎fdaPDE/linear_algebra.h‎

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,7 @@ template <int Rows, int Cols, typename XprType> struct MatrixExpr;
4040

4141
namespace internals {
4242

43-
// detects whether XprType represents a static sized or dynamic sized expression
43+
// sizing traits
4444
template <typename XprType> struct is_dynamic_sized {
4545
static constexpr bool value = std::decay_t<XprType>::Rows == Dynamic || std::decay_t<XprType>::Cols == Dynamic;
4646
};
@@ -49,6 +49,18 @@ template <typename XprType> struct is_static_sized {
4949
static constexpr bool value = !is_dynamic_sized_v<XprType>;
5050
};
5151
template <typename XprType> static constexpr bool is_static_sized_v = is_static_sized<XprType>::value;
52+
template <typename LhsXprType, typename RhsXprType> struct same_static_size {
53+
private:
54+
using LhsXprTypeClean = std::decay_t<LhsXprType>;
55+
using RhsXprTypeClean = std::decay_t<RhsXprType>;
56+
public:
57+
static constexpr bool value =
58+
!is_dynamic_sized_v<LhsXprTypeClean> && !is_dynamic_sized_v<RhsXprTypeClean> &&
59+
(LhsXprTypeClean::Rows * LhsXprTypeClean::Cols == RhsXprTypeClean::Rows * RhsXprTypeClean::Cols);
60+
};
61+
template <typename LhsXprType, typename RhsXprType>
62+
static constexpr bool same_static_size_v = same_static_size<LhsXprType, RhsXprType>::value;
63+
// shaping traits
5264
template <typename LhsXprType, typename RhsXprType> struct same_static_shape {
5365
private:
5466
using LhsXprTypeClean = std::decay_t<LhsXprType>;
@@ -60,11 +72,17 @@ template <typename LhsXprType, typename RhsXprType> struct same_static_shape {
6072
};
6173
template <typename LhsXprType, typename RhsXprType>
6274
static constexpr bool same_static_shape_v = same_static_shape<LhsXprType, RhsXprType>::value;
63-
75+
template <typename XprType> struct is_vector_shaped {
76+
using XprTypeClean = std::decay_t<XprType>;
77+
static constexpr bool value = (XprTypeClean::Rows == 1 || XprTypeClean::Cols == 1);
78+
};
79+
template <typename XprType> static constexpr bool is_vector_shaped_v = is_vector_shaped<XprType>::value;
80+
6481
} // namespace internals
6582
} // namespace fdapde
6683

6784
#include "src/linear_algebra/matrix.h"
85+
#include "src/linear_algebra/coeffwise_op.h"
6886
#include "src/linear_algebra/binary_op.h"
6987
#include "src/linear_algebra/block.h"
7088
#include "src/linear_algebra/unary_op.h"
@@ -77,15 +95,18 @@ static constexpr bool same_static_shape_v = same_static_shape<LhsXprType, RhsXpr
7795
#include "src/linear_algebra/bool.h"
7896

7997
// #include "src/linear_algebra/skew.h"
80-
// #include "src/linear_algebra/permutation.h"
98+
#include "src/linear_algebra/permutation.h"
8199
// #include "src/linear_algebra/spd.h"
82100

101+
102+
#include "src/linear_algebra/partial_piv_lu.h"
103+
83104
// expression template system
84105
#include "src/linear_algebra/xpr.h"
85106

86107
// algorithms
87108
#include "src/linear_algebra/evd.h"
88-
// #include "src/linear_algebra/partial_piv_lu.h"
109+
89110
// #include "src/linear_algebra/qr.h"
90111

91112
// clang-format on

‎fdaPDE/src/geometry/dcel.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ template <int LocalDim, int EmbedDim> class DCEL {
3434
// internal data structures
3535
struct node_t {
3636
private:
37-
using coords_t = Eigen::Matrix<double, embed_dim, 1>;
37+
using coords_t = Matrix<double, embed_dim, 1>;
3838
int id_; // global node index
3939
halfedge_t* halfedge_; // any edge having this node as its origin
4040
bool boundary_; // asserted true if node is on boundary
@@ -65,7 +65,7 @@ template <int LocalDim, int EmbedDim> class DCEL {
6565
node_t(int id, bool boundary, CoordsType&&... coords) :
6666
node_t(id, nullptr, boundary, coords...) { }
6767
// observers
68-
const Eigen::Matrix<double, embed_dim, 1>& coords() const { return coords_; }
68+
const Matrix<double, embed_dim, 1>& coords() const { return coords_; }
6969
halfedge_t* halfedge() const { return halfedge_; }
7070
void set_halfedge(halfedge_t* halfedge) { halfedge_ = halfedge; }
7171
int id() const { return id_; }
@@ -158,7 +158,7 @@ template <int LocalDim, int EmbedDim> class DCEL {
158158
// constructors
159159
DCEL() : nodes_(), halfedges_(), n_nodes_(0), n_halfedges_(0), n_cells_(0) { }
160160
// constructs a closed loop structure linking nodes one after the other
161-
static DCEL<local_dim, embed_dim> make_polygon(const Eigen::Matrix<double, Dynamic, Dynamic>& nodes) {
161+
static DCEL<local_dim, embed_dim> make_polygon(const Matrix<double, Dynamic, Dynamic>& nodes) {
162162
fdapde_assert(nodes.cols() == embed_dim);
163163
int n_nodes = nodes.rows();
164164
DCEL<local_dim, embed_dim> dcel;
@@ -239,8 +239,8 @@ template <int LocalDim, int EmbedDim> class DCEL {
239239
}
240240

241241
// observers
242-
Eigen::Matrix<double, Dynamic, Dynamic> nodes() const { // matrix of nodes coordinates
243-
Eigen::Matrix<double, Dynamic, Dynamic> coords(n_nodes_, embed_dim);
242+
Matrix<double, Dynamic, Dynamic> nodes() const { // matrix of nodes coordinates
243+
Matrix<double, Dynamic, Dynamic> coords(n_nodes_, embed_dim);
244244
for (auto it = nodes_.begin(); it != nodes_.end(); ++it) { coords.row(it->id()) = it->coords(); }
245245
return coords;
246246
}

‎fdaPDE/src/geometry/hyperplane.h‎

Lines changed: 32 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,9 @@ template <int LocalDim, int EmbedDim> class HyperPlane {
2929
static constexpr int local_dim = LocalDim;
3030
HyperPlane() = default;
3131
// constructs a hyperplane passing through 2 points, e.g. a line
32-
HyperPlane(const Eigen::Matrix<double, embed_dim, 1>& x1, const Eigen::Matrix<double, embed_dim, 1>& x2) : p_(x1) {
32+
HyperPlane(const Matrix<double, embed_dim, 1>& x1, const Matrix<double, embed_dim, 1>& x2) : p_(x1) {
3333
fdapde_static_assert(local_dim == 1, THIS_METHOD_IS_ONLY_FOR_LINES);
34-
Eigen::Matrix<double, embed_dim, 1> tmp = x2 - x1;
34+
Matrix<double, embed_dim, 1> tmp = x2 - x1;
3535
basis_.col(0) = tmp.normalized();
3636
if constexpr (embed_dim == 2) normal_ << -tmp[0], tmp[1];
3737
if constexpr (embed_dim == 3) normal_ << -tmp[0], tmp[1], 0; // one of the (infintely-many) normals to 3D line
@@ -40,28 +40,31 @@ template <int LocalDim, int EmbedDim> class HyperPlane {
4040
}
4141
// constructs an hyperplane passing through 3 (non-collinear) points, e.g. a plane
4242
HyperPlane(
43-
const Eigen::Matrix<double, embed_dim, 1>& x1, const Eigen::Matrix<double, embed_dim, 1>& x2,
44-
const Eigen::Matrix<double, embed_dim, 1>& x3) :
43+
const Matrix<double, embed_dim, 1>& x1, const Matrix<double, embed_dim, 1>& x2,
44+
const Matrix<double, embed_dim, 1>& x3) :
4545
p_(x1) {
4646
fdapde_static_assert(local_dim == 2, THIS_METHOD_IS_ONLY_FOR_PLANES);
4747
basis_.col(0) = (x2 - x1).normalized();
48-
basis_.col(1) = ((x3 - x1) - orthogonal_project(
49-
Eigen::Matrix<double, embed_dim, 1>(x3 - x1),
50-
Eigen::Matrix<double, embed_dim, 1>(basis_.col(0))))
51-
.normalized();
48+
basis_.col(1) =
49+
((x3 - x1) -
50+
orthogonal_project(Matrix<double, embed_dim, 1>(x3 - x1), Matrix<double, embed_dim, 1>(basis_.col(0))))
51+
.normalized();
5252
normal_ = ((x2 - x1).cross(x3 - x1)).normalized();
5353
offset_ = -x1.dot(normal_);
5454
}
5555
// constructors from matrix coordinates
56-
HyperPlane(const Eigen::Matrix<double, embed_dim, 2>& coords) requires(local_dim == 1) :
57-
HyperPlane(coords.col(0), coords.col(1)) { }
58-
HyperPlane(const Eigen::Matrix<double, embed_dim, 3>& coords) requires (local_dim == 2) :
59-
HyperPlane(coords.col(0), coords.col(1), coords.col(2)) { }
56+
HyperPlane(const Matrix<double, embed_dim, 2>& coords)
57+
requires(local_dim == 1)
58+
: HyperPlane(coords.col(0), coords.col(1)) { }
59+
HyperPlane(const Matrix<double, embed_dim, 3>& coords)
60+
requires(local_dim == 2)
61+
: HyperPlane(coords.col(0), coords.col(1), coords.col(2)) { }
6062
// general hyperplane constructor
61-
HyperPlane(const Eigen::Matrix<double, embed_dim, local_dim + 1>& coords) requires(local_dim > 2) :
62-
p_(coords.col(0)) {
63+
HyperPlane(const Matrix<double, embed_dim, local_dim + 1>& coords)
64+
requires(local_dim > 2)
65+
: p_(coords.col(0)) {
6366
basis_ = coords.rightCols(local_dim).colwise() - coords.col(0);
64-
// basis orthonormalization via modified Gram-Schmidt method
67+
// basis orthonormalization via modified Gram-Schmidt method
6568
basis_.col(0) /= basis_.col(0).norm();
6669
for (int i = 1; i < local_dim; ++i) {
6770
for (int j = 0; j < i; ++j) {
@@ -73,46 +76,43 @@ template <int LocalDim, int EmbedDim> class HyperPlane {
7376
offset_ = -coords.col(0).dot(normal_);
7477
}
7578
// projection
76-
Eigen::Matrix<double, local_dim, 1> project_onto(const Eigen::Matrix<double, embed_dim, 1>& x) {
79+
Matrix<double, local_dim, 1> project_onto(const Matrix<double, embed_dim, 1>& x) {
7780
if constexpr (local_dim == embed_dim) {
7881
return x;
7982
} else {
8083
// build the projection onto the space spanned by basis_
81-
Eigen::Matrix<double, local_dim, 1> proj;
84+
Matrix<double, local_dim, 1> proj;
8285
for (int i = 0; i < local_dim; ++i) { proj[i] = (x - p_).dot(basis_.col(i)); }
8386
return proj;
8487
}
8588
}
86-
Eigen::Matrix<double, embed_dim, 1> project(const Eigen::Matrix<double, embed_dim, 1>& x) const {
89+
Matrix<double, embed_dim, 1> project(const Matrix<double, embed_dim, 1>& x) const {
8790
if constexpr (local_dim == embed_dim) {
8891
return x;
8992
} else {
9093
return basis_ * basis_.transpose() * (x - p_) + p_;
9194
}
9295
}
93-
double distance(const Eigen::Matrix<double, embed_dim, 1>& x) {
94-
return (x - project(x)).norm();
95-
} // point-plane distance
96-
double eval(const Eigen::Matrix<double, embed_dim, 1>& coeffs) const { return normal_.dot(coeffs) + offset_; }
97-
Eigen::Matrix<double, embed_dim, 1> operator()(const Eigen::Matrix<double, local_dim, 1>& coeffs) const {
98-
Eigen::Matrix<double, embed_dim, 1> res = p_;
96+
double distance(const Matrix<double, embed_dim, 1>& x) { return (x - project(x)).norm(); } // point-plane distance
97+
double eval(const Matrix<double, embed_dim, 1>& coeffs) const { return normal_.dot(coeffs) + offset_; }
98+
Matrix<double, embed_dim, 1> operator()(const Matrix<double, local_dim, 1>& coeffs) const {
99+
Matrix<double, embed_dim, 1> res = p_;
99100
for (int i = 0; i < local_dim; ++i) { res += coeffs[i] * basis_.col(i); }
100101
return res;
101102
}
102-
const Eigen::Matrix<double, embed_dim, 1>& normal() const { return normal_; } // plane normal direction
103-
const Eigen::Matrix<double, embed_dim, 1>& point() const { return p_; } // a point belonging to the plane
104-
const Eigen::Matrix<double, embed_dim, local_dim>& basis() const { return basis_; }
103+
const Matrix<double, embed_dim, 1>& normal() const { return normal_; } // plane normal direction
104+
const Matrix<double, embed_dim, 1>& point() const { return p_; } // a point belonging to the plane
105+
const Matrix<double, embed_dim, local_dim>& basis() const { return basis_; }
105106
private:
106107
// let x_1, x_2, \ldots, x_{N+1} be a set of N+1 points through which the plane passes
107-
Eigen::Matrix<double, embed_dim, local_dim> basis_; // matrix [x2 - x1, x3 - x1, \ldots, x_{N+1} - x1]
108-
Eigen::Matrix<double, embed_dim, 1> normal_; // normal vector to the hyperplane
109-
Eigen::Matrix<double, embed_dim, 1> p_; // point through which the plane is guaranteed to pass
108+
Matrix<double, embed_dim, local_dim> basis_; // matrix [x2 - x1, x3 - x1, \ldots, x_{N+1} - x1]
109+
Matrix<double, embed_dim, 1> normal_; // normal vector to the hyperplane
110+
Matrix<double, embed_dim, 1> p_; // point through which the plane is guaranteed to pass
110111
double offset_; // hyperplane's intercept (the d in the equation ax + by + cz + d = 0)
111112

112113
// orthogonal projection of vector v over u
113114
template <int N>
114-
constexpr Eigen::Matrix<double, N, 1>
115-
orthogonal_project(const Eigen::Matrix<double, N, 1>& v, const Eigen::Matrix<double, N, 1>& u) {
115+
constexpr Matrix<double, N, 1> orthogonal_project(const Matrix<double, N, 1>& v, const Matrix<double, N, 1>& u) {
116116
return (v.dot(u) / u.squaredNorm() * u.array()).matrix();
117117
}
118118
};

‎fdaPDE/src/geometry/interval.h‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula
3333
using Base::nodes_; // physical coordinates of nodes
3434

3535
Triangulation() = default;
36-
Triangulation(const Eigen::Matrix<double, Dynamic, 1>& nodes, int flags = 0) : Base() {
36+
Triangulation(const Matrix<double, Dynamic, 1>& nodes, int flags = 0) : Base() {
3737
fdapde_assert(nodes.rows() > 1 && nodes.cols() == 1);
3838
Base::flags_ = flags;
3939
nodes_ = nodes;
@@ -49,7 +49,7 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula
4949
cells_(i, 0) = i;
5050
cells_(i, 1) = i + 1;
5151
}
52-
neighbors_ = Eigen::Matrix<int, Dynamic, Dynamic>::Constant(n_cells_, n_neighbors_per_cell, -1);
52+
neighbors_ = Matrix<int, Dynamic, Dynamic>::Constant(n_cells_, n_neighbors_per_cell, -1);
5353
neighbors_(0, 1) = 1;
5454
for (int i = 1; i < n_cells_ - 1; ++i) {
5555
neighbors_(i, 0) = i - 1;
@@ -62,7 +62,7 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula
6262
Base::boundary_markers_.set(n_nodes_ - 1);
6363
};
6464
// construct from interval's bounds [a, b] and the number of equidistant nodes n used to split [a, b]
65-
Triangulation(double a, double b, int n) : Triangulation(Eigen::Matrix<double, Dynamic, 1>::LinSpaced(n, a, b)) { }
65+
Triangulation(double a, double b, int n) : Triangulation(Matrix<double, Dynamic, 1>::LinSpaced(n, a, b)) { }
6666
Triangulation(const std::string& nodes, bool header, bool index_col, int flags = 0) :
6767
Triangulation(read_table<double>(nodes, header, index_col).as_matrix(), flags) { }
6868

@@ -71,7 +71,7 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula
7171
static Triangulation<1, 1> UnitInterval(int n_nodes) { return Triangulation<1, 1>::Interval(0.0, 1.0, n_nodes); }
7272

7373
// getters
74-
const Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor>& neighbors() const { return neighbors_; }
74+
const Matrix<int, Dynamic, Dynamic>& neighbors() const { return neighbors_; }
7575
const typename Base::CellType& cell(int id) const {
7676
if (Base::flags_ & cache_cells) { // cell caching enabled
7777
return cell_cache_[id];
@@ -97,16 +97,16 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula
9797
}
9898
// point location
9999
template <int Rows, int Cols>
100-
std::conditional_t<Rows == Dynamic || Cols == Dynamic, Eigen::Matrix<int, Dynamic, 1>, int>
101-
locate(const Eigen::Matrix<double, Rows, Cols>& p) const {
100+
std::conditional_t<Rows == Dynamic || Cols == Dynamic, Matrix<int, Dynamic, 1>, int>
101+
locate(const Matrix<double, Rows, Cols>& p) const {
102102
fdapde_static_assert(
103103
(Cols == 1 && Rows == 1) || (Cols == Dynamic && Rows == Dynamic),
104104
YOU_PASSED_A_MATRIX_OF_POINTS_TO_LOCATE_OF_WRONG_DIMENSIONS);
105105
if constexpr (Cols == 1) {
106106
return locate_(p[0]);
107107
} else {
108108
fdapde_assert(p.rows() > 0 && p.cols() == 1);
109-
Eigen::Matrix<int, Dynamic, 1> result;
109+
Matrix<int, Dynamic, 1> result;
110110
result.resize(p.rows());
111111
// start search
112112
for (int i = 0; i < p.rows(); ++i) { result[i] = locate_(p(i, 0)); }
@@ -146,7 +146,7 @@ template <> class Triangulation<1, 1> : public TriangulationBase<1, 1, Triangula
146146
}
147147
return -1;
148148
}
149-
Eigen::Matrix<int, Dynamic, Dynamic, Eigen::RowMajor> neighbors_ {}; // adjacent cells ids (-1: no adjacent cell)
149+
Matrix<int, Dynamic, Dynamic> neighbors_ {}; // adjacent cells ids (-1: no adjacent cell)
150150
// cell caching
151151
std::vector<typename Base::CellType> cell_cache_;
152152
mutable typename Base::CellType cell_; // used in case cell caching is off

‎fdaPDE/src/geometry/kd_tree.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ template <int K> class KDTree {
2727
using Container = BinaryTree<int>;
2828
using node_type = Container::node_type;
2929
using node_pointer = Container::node_pointer;
30-
using data_type = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
30+
using data_type = Matrix<double, Dynamic, Dynamic>;
3131
Container kdtree_; // the actual BinaryTree container
3232
data_type data_; // set of data indexed by this tree
3333
public:
@@ -73,7 +73,7 @@ template <int K> class KDTree {
7373
iterator end() { return kdtree_.end(); }
7474

7575
// returns an iterator to the nearest neighbor of p. Average O(log(n)) complexity (worst case is O(n))
76-
iterator nn_search(const Eigen::Matrix<double, Dynamic, 1>& p) const {
76+
iterator nn_search(const Matrix<double, Dynamic, 1>& p) const {
7777
fdapde_assert(p.size() == K);
7878
if (kdtree_.empty()) return kdtree_.cend(); // nothing to search
7979
const data_type& data = data_;
@@ -113,7 +113,7 @@ template <int K> class KDTree {
113113

114114
// solves a (rectangular) range query in a K-dimensional euclidean space
115115
struct RangeType {
116-
Eigen::Matrix<double, K, 1> ll, ur; // lower-left and upper-right corner
116+
Matrix<double, K, 1> ll, ur; // lower-left and upper-right corner
117117
};
118118
// returns a set of iterators to the nodes contained in the query
119119
std::unordered_set<int> range_search(const RangeType& query) const {

‎fdaPDE/src/geometry/segment.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,8 +38,8 @@ template <typename Triangulation> class Segment : public Simplex<Triangulation::
3838
}
3939
// getters
4040
int id() const { return id_; }
41-
Eigen::Matrix<int, Dynamic, 1> neighbors() const { return mesh_->neighbors().row(id_); }
42-
Eigen::Matrix<int, Dynamic, 1> node_ids() const { return mesh_->cells().row(id_); }
41+
Matrix<int, Dynamic, 1> neighbors() const { return mesh_->neighbors().row(id_); }
42+
Matrix<int, Dynamic, 1> node_ids() const { return mesh_->cells().row(id_); }
4343
bool on_boundary() const { return boundary_; }
4444
operator bool() const { return mesh_ != nullptr; }
4545
protected:

0 commit comments

Comments
 (0)