Skip to content

Commit

Permalink
warning fixes, small changes in field module and BlockFrame, pde_ptr …
Browse files Browse the repository at this point in the history 
…removed (core exposes just erasable wrapper)
  • Loading branch information
@AlePalu
AlePalu committed Mar 8, 2024
1 parent cd138ee commit e3f061d
Show file tree
Hide file tree
Showing 15 changed files with 265 additions and 350 deletions.
8 changes: 4 additions & 4 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,9 @@ Documentation can be found on our [documentation site](https://fdapde.github.io/

## Dependencies
fdaPDE-core is an header-only library, therefore it does not require any installation. Just make sure to have it in your include path. Neverthless, to compile code including this library you need:
* A C++17 compliant compiler. Supported versions are:
* Linux: `gcc` 11 (or higher), `clang` 13 (or higher)
* macOS: `apple-clang` (the XCode version of `clang`).
* The [Eigen](https://eigen.tuxfamily.org/index.php?title=Main_Page) linear algebra library, version 3.3.9.
* A C++20 compliant compiler. Supported versions are:
* Linux: `gcc` 11 (or higher), `clang` 15 (or higher)
* macOS: `apple-clang` (the XCode version of `clang`, AppleClang 15 or higher).
* The [Eigen](https://eigen.tuxfamily.org/index.php?title=Main_Page) linear algebra library, version 3.4.0.

If you wish to run the test suite contained in the `test/` folder, be sure to have [Google Test](http://google.github.io/googletest/) installed.
2 changes: 1 addition & 1 deletion fdaPDE/fields/dot_product.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ template <int N, typename T1, typename T2> class DotProduct : public ScalarExpr<
}
return result;
}
template <typename T> const DotProduct<N, T1, T2>& forward(T i) {
template <typename T> const DotProduct<N, T1, T2>& forward(T i) const {
op1_.forward(i);
op2_.forward(i);
return *this;
Expand Down
6 changes: 3 additions & 3 deletions fdaPDE/fields/field_ptrs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ template <typename E> class ScalarPtr : public ScalarExpr<E::static_inner_size,
ScalarPtr(E* ptr) : ptr_(ptr) { fdapde_assert(bool(ptr_) == true); };
// delegate to pointed memory location
double operator()(const SVector<E::static_inner_size>& p) const { return ptr_->operator()(p); }
template <typename T> void forward(T i) { ptr_->forward(i); }
template <typename T> void forward(T i) const { ptr_->forward(i); }
// access to pointed element
E* operator->() { return ptr_; }
typedef E PtrType; // expose wrapped type
Expand All @@ -50,7 +50,7 @@ template <typename E> class VectorPtr : public VectorExpr<E::static_inner_size,
VectorPtr(E* ptr) : ptr_(ptr) { fdapde_assert(bool(ptr_) == true); };
// delegate to pointed memory location
auto operator[](std::size_t i) const { return ptr_->operator[](i); }
template <typename T> void forward(T i) { ptr_->forward(i); }
template <typename T> void forward(T i) const { ptr_->forward(i); }
// access to pointed element
E* operator->() { return ptr_; }
typedef E PtrType; // expose wrapped type
Expand All @@ -66,7 +66,7 @@ template <typename E> class MatrixPtr : public MatrixExpr<E::static_inner_size,
MatrixPtr(E* ptr) : ptr_(ptr) { fdapde_assert(bool(ptr_) == true); };
// delegate to pointed memory location
auto coeff(std::size_t i, std::size_t j) const { return ptr_->coeff(i, j); }
template <typename T> void forward(T i) { ptr_->forward(i); }
template <typename T> void forward(T i) const { ptr_->forward(i); }
// access to pointed element
E* operator->() { return ptr_; }
typedef E PtrType; // expose wrapped type
Expand Down
51 changes: 23 additions & 28 deletions fdaPDE/fields/matrix_expressions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ class MatrixVectorProduct : public VectorExpr<N, M, MatrixVectorProduct<N, M, K,
inline DotProduct<N, decltype(op1_.row(int())), T2> operator[](int i) const {
return DotProduct<N, decltype(op1_.row(int())), T2>(op1_.row(i), op2_);
}
template <typename T> const MatrixVectorProduct<N, M, K, T1, T2>& forward(T i) {
template <typename T> const MatrixVectorProduct<N, M, K, T1, T2>& forward(T i) const {
op1_.forward(i);
op2_.forward(i);
return *this;
Expand All @@ -69,7 +69,7 @@ class MatrixMatrixProduct : public MatrixExpr<N, M, K, MatrixMatrixProduct<N, M,
return DotProduct<N, decltype(op1_.row(int())), decltype(op2_.col(int()))>(
op1_.row(i), op2_.col(j));
}
template <typename T> const MatrixMatrixProduct<N, M, K, T1, T2>& forward(T i) {
template <typename T> const MatrixMatrixProduct<N, M, K, T1, T2>& forward(T i) const {
op1_.forward(i);
op2_.forward(i);
return *this;
Expand All @@ -85,7 +85,7 @@ template <int N, int M, int K, typename E> class MatrixRow : public VectorExpr<N
MatrixRow(const E& expr, int row) : expr_(expr), row_(row) {};
// subscripting the i-th element of a row triggers evaluation of the expression
auto operator[](int i) const { return expr_.coeff(row_, i); }
template <typename T> const MatrixRow<N, M, K, E>& forward(T i) {
template <typename T> const MatrixRow<N, M, K, E>& forward(T i) const {
expr_.forward(i);
return *this;
}
Expand All @@ -102,7 +102,7 @@ template <int N, int M, int K, typename E> class MatrixCol : public VectorExpr<N
MatrixCol(const E& expr, int col) : expr_(expr), col_(col) {};
// subscripting the i-th element of a column triggers evaluation of the expression
auto operator[](int i) const { return expr_.coeff(i, col_); }
template <typename T> const MatrixCol<N, M, K, E>& forward(T i) {
template <typename T> const MatrixCol<N, M, K, E>& forward(T i) const {
expr_.forward(i);
return *this;
}
Expand All @@ -119,19 +119,15 @@ template <int N, int M, int K, typename E> class MatrixCol : public VectorExpr<N
const MatrixExpr<N, M, K, E1>& op1, const MatrixExpr<N, M, K, E2>& op2) { \
return MatrixBinOp<N, M, K, E1, E2, FUNCTOR>(op1.get(), op2.get(), FUNCTOR()); \
} \
\
template <int N, int M, int K, typename E> \
MatrixBinOp<N, M, K, MatrixConst<N, M, K>, E, FUNCTOR> OPERATOR( \
MatrixBinOp<N, M, K, Matrix<N, M, K>, E, FUNCTOR> OPERATOR( \
SMatrix<M, K> op1, const MatrixExpr<N, M, K, E>& op2) { \
return MatrixBinOp<N, M, K, MatrixConst<N, M, K>, E, FUNCTOR>( \
MatrixConst<N, M, K>(op1), op2.get(), FUNCTOR()); \
return MatrixBinOp<N, M, K, Matrix<N, M, K>, E, FUNCTOR>(Matrix<N, M, K>(op1), op2.get(), FUNCTOR()); \
} \
\
template <int N, int M, int K, typename E> \
MatrixBinOp<N, M, K, E, MatrixConst<N, M, K>, FUNCTOR> OPERATOR( \
MatrixBinOp<N, M, K, E, Matrix<N, M, K>, FUNCTOR> OPERATOR( \
const MatrixExpr<N, M, K, E>& op1, SMatrix<M, K> op2) { \
return MatrixBinOp<N, M, K, E, MatrixConst<N, M, K>, FUNCTOR>( \
op1.get(), MatrixConst<N, M, K>(op2), FUNCTOR()); \
return MatrixBinOp<N, M, K, E, Matrix<N, M, K>, FUNCTOR>(op1.get(), Matrix<N, M, K>(op2), FUNCTOR()); \
}

// base class for matrix expressions
Expand All @@ -150,7 +146,7 @@ template <int M, int N, int K, typename E> class MatrixExpr : public MatrixBase
MatrixExpr() = default;
MatrixExpr(int inner_size, int outer_rows, int outer_cols) :
inner_size_(inner_size), outer_rows_(outer_rows), outer_cols_(outer_cols),
outer_size_(outer_rows * outer_cols) {};
outer_size_(outer_rows * outer_cols) { }

// access operator on (i,j)-th element on the base type E
auto coeff(int i, int j) const { return static_cast<const E&>(*this).coeff(i, j); }
Expand All @@ -176,7 +172,7 @@ template <int M, int N, int K, typename E> class MatrixExpr : public MatrixBase
MatrixRow<M, N, K, E> row(int i) const;
MatrixCol<M, N, K, E> col(int i) const;
template <typename F> MatrixVectorProduct<M, N, K, E, F> operator*(const VectorExpr<M, K, F>& op) const;
MatrixVectorProduct<M, N, K, E, VectorConst<M, K>> operator*(const SVector<K>& op) const;
MatrixVectorProduct<M, N, K, E, Vector<M, K>> operator*(const SVector<K>& op) const;
template <typename T> void forward([[maybe_unused]] T i) const { return; }
MatrixNegationOp<M, N, K, E> operator-() const { return MatrixNegationOp<M, N, K, E>(get()); }
};
Expand All @@ -192,16 +188,15 @@ template <int N, int M, int K, typename E> MatrixCol<N, M, K, E> MatrixExpr<N, M

// an expression node representing a constant matrix
template <int N, int M, int K>
class MatrixConst : public MatrixExpr<N, M, K, MatrixConst<N, M, K>> {
class Matrix : public MatrixExpr<N, M, K, Matrix<N, M, K>> {
private:
SMatrix<M, K> value_;
public:
// constructor
MatrixConst() = default;
MatrixConst(SMatrix<M, K> value) : value_(value) { }
Matrix() = default;
Matrix(SMatrix<M, K> value) : value_(value) { }
double coeff(int i, int j) const { return value_(i, j); }
// assignment operator
MatrixConst<N, M, K>& operator=(const SMatrix<M, K>& value) {
Matrix<N, M, K>& operator=(const SMatrix<M, K>& value) {
value_ = value;
return *this;
}
Expand All @@ -212,15 +207,15 @@ class MatrixConst : public MatrixExpr<N, M, K, MatrixConst<N, M, K>> {
template <int N, int M, int K>
class DiscretizedMatrixField : public MatrixExpr<N, M, K, DiscretizedMatrixField<N, M, K>> {
private:
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DataType;
using DataType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
DataType* data_;
Eigen::Map<SMatrix<M, K>> value_;
mutable Eigen::Map<SMatrix<M, K>> value_;
public:
// constructor
DiscretizedMatrixField() : value_(NULL) {};
DiscretizedMatrixField(DataType& data) : data_(&data), value_(NULL) {};
double coeff(int i, int j) const { return value_(i, j); }
void forward(int i) {
void forward(int i) const {
new (&value_) Eigen::Map<SMatrix<M, K>>(data_->data() + (i * M * K)); // construct map in place
}
};
Expand All @@ -238,7 +233,7 @@ class MatrixBinOp : public MatrixExpr<N, M, K, MatrixBinOp<N, M, K, OP1, OP2, Bi
// access operator. Apply the functor to each accessed element. This returns a callable object
auto coeff(int i, int j) const { return f_(op1_.coeff(i, j), op2_.coeff(i, j)); }
// call parameter evaluation on operands
template <typename T> const MatrixBinOp<N, M, K, OP1, OP2, BinaryOperation>& forward(T i) {
template <typename T> const MatrixBinOp<N, M, K, OP1, OP2, BinaryOperation>& forward(T i) const {
op1_.forward(i);
op2_.forward(i);
return *this;
Expand All @@ -255,14 +250,14 @@ MatrixVectorProduct<N, M, K, E, F> MatrixExpr<N, M, K, E>::operator*(const Vecto
}
// SMatrix<M,K> times VectorExpr<N,K>
template <int N, int M, int K, typename F>
MatrixVectorProduct<N, M, K, MatrixConst<N, M, K>, F>
MatrixVectorProduct<N, M, K, Matrix<N, M, K>, F>
operator*(const Eigen::Matrix<double, M, K>& op1, const VectorExpr<N, K, F>& op2) {
return MatrixVectorProduct<N, M, K, MatrixConst<N, M, K>, F>(MatrixConst<N, M, K>(op1), op2.get());
return MatrixVectorProduct<N, M, K, Matrix<N, M, K>, F>(Matrix<N, M, K>(op1), op2.get());
}
// MatrixExpr<N,M,K> times SVector<K>
template <int N, int M, int K, typename E>
MatrixVectorProduct<N, M, K, E, VectorConst<N, K>> MatrixExpr<N, M, K, E>::operator*(const SVector<K>& op) const {
return MatrixVectorProduct<N, M, K, E, VectorConst<N, K>>(this->get(), VectorConst<N, K>(op));
MatrixVectorProduct<N, M, K, E, Vector<N, K>> MatrixExpr<N, M, K, E>::operator*(const SVector<K>& op) const {
return MatrixVectorProduct<N, M, K, E, Vector<N, K>>(this->get(), Vector<N, K>(op));
}

// element wise multiplication of MatrixExpr by scalar
Expand Down Expand Up @@ -303,7 +298,7 @@ class MatrixNegationOp : public MatrixExpr<M, N, K, MatrixNegationOp<M, N, K, E>
MatrixNegationOp(const E& op) : op_(op) {};
auto coeff(int i, int j) const { return -(op_.coeff(i, j)); }
// call parameter evaluation on operands
template <typename T> const MatrixNegationOp<N, M, K, E>& forward(T i) {
template <typename T> const MatrixNegationOp<N, M, K, E>& forward(T i) const {
op_.forward(i);
return *this;
}
Expand Down
30 changes: 14 additions & 16 deletions fdaPDE/fields/matrix_field.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,28 +46,26 @@ class MatrixField : public MatrixExpr<M, N, K, MatrixField<M, N, K, F>> {
std::is_same<typename std::invoke_result<F, InnerVectorType>::type, double>::value);
// static constructors
MatrixField() requires(N != Dynamic && K != Dynamic) { field_.resize(N * K); }
MatrixField(const std::vector<FieldType>& v) requires(N != Dynamic && K != Dynamic){
fdapde_assert(int(v.size()) == outer_size());
field_.reserve(v.size());
for (std::size_t i = 0; i < v.size(); ++i) { field_.emplace_back(v[i]); }
MatrixField(const std::vector<FieldType>& v) requires(N != Dynamic && K != Dynamic) {
fdapde_assert(int(v.size()) == outer_size());
field_.reserve(v.size());
for (std::size_t i = 0; i < v.size(); ++i) { field_.emplace_back(v[i]); }
}
// fully dynamic constructor
MatrixField(int m, int n, int k) requires (N == Dynamic || K == Dynamic): Base(m, n, k) {
field_.resize(n * k, ScalarField<M, FieldType>(m));
}
// call operator (evaluate field at point)
OuterMatrixType operator()(const InnerVectorType& point) const {
OuterMatrixType result;
if constexpr(N == Dynamic || K == Dynamic) result.resize(outer_rows(), outer_cols());
for (int i = 0; i < outer_rows(); ++i) {
for (int j = 0; j < outer_cols(); ++j) result(i, j) = field_[j + i * outer_cols()](point);
}
return result;
OuterMatrixType result;
if constexpr (N == Dynamic || K == Dynamic) result.resize(outer_rows(), outer_cols());
for (int i = 0; i < outer_rows(); ++i) {
for (int j = 0; j < outer_cols(); ++j) result(i, j) = field_[j + i * outer_cols()](point);
}
return result;
}
// const and non-const access operations
const ScalarField<M, FieldType>& operator()(int i, int j) const {
return field_[j + i * outer_cols()];
};
const ScalarField<M, FieldType>& operator()(int i, int j) const { return field_[j + i * outer_cols()]; };
ScalarField<M, FieldType>& operator()(int i, int j) { return field_[j + i * outer_cols()]; };
const ScalarField<M, FieldType>& coeff(int i, int j) const { return operator()(i, j); };
protected:
Expand All @@ -77,10 +75,10 @@ class MatrixField : public MatrixExpr<M, N, K, MatrixField<M, N, K, F>> {
// out of class definitions of MatrixField arithmetic
// rhs multiplication by SVector
template <int M, int N, int K, typename F>
MatrixVectorProduct<M, N, K, MatrixField<M, N, K, F>, VectorConst<M, K>>
MatrixVectorProduct<M, N, K, MatrixField<M, N, K, F>, Vector<M, K>>
operator*(const MatrixField<M, N, K, F>& op1, const static_dynamic_vector_selector_t<K>& op2) {
return MatrixVectorProduct<M, N, K, MatrixField<M, N, K, F>, VectorConst<M, K>>(
op1, VectorConst<M, K>(op2, op1.inner_size(), op1.outer_cols()));
return MatrixVectorProduct<M, N, K, MatrixField<M, N, K, F>, Vector<M, K>>(
op1, Vector<M, K>(op2, op1.inner_size(), op1.outer_cols()));
}
// rhs multiplication by VectorField
template <int M, int N, int K, typename F1, typename F2>
Expand Down
5 changes: 1 addition & 4 deletions fdaPDE/fields/scalar_expressions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,10 +64,8 @@ template <int N, typename E> class ScalarExpr : public ScalarBase {
static constexpr int cols = 1;
static constexpr int static_inner_size = N; // dimensionality of base space (can be Dynamic)
static constexpr int NestAsRef = 0; // whether to store the node by reference of by copy

ScalarExpr() = default;
ScalarExpr(int dynamic_inner_size) : dynamic_inner_size_(dynamic_inner_size) { }

// call operator() on the base type E
inline double operator()(const VectorType& p) const { return static_cast<const E&>(*this)(p); }
const E& get() const { return static_cast<const E&>(*this); }
Expand All @@ -77,8 +75,7 @@ template <int N, typename E> class ScalarExpr : public ScalarBase {
ScalarNegationOp<N, E> operator-() const { return ScalarNegationOp<N, E>(get(), dynamic_inner_size_); }
inline constexpr int inner_size() const { return (N == Dynamic) ? dynamic_inner_size_ : static_inner_size; }
// dynamic resizing of base space dimension only allowed for Dynamic expressions
template <int N_ = N> typename std::enable_if<N_ == Dynamic, void>::type resize(int n) { dynamic_inner_size_ = n; }

void resize(int n) requires(N == Dynamic) { dynamic_inner_size_ = n; }
void set_step(double h) { h_ = h; } // set step size in derivative approximation
double step() const { return h_; }
ScalarExprGradient<N, E> derive() const { return ScalarExprGradient<N, E>(get(), h_); }
Expand Down
26 changes: 11 additions & 15 deletions fdaPDE/fields/scalar_field.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,28 +56,24 @@ class ScalarField : public ScalarExpr<N, ScalarField<N, F>> {
explicit ScalarField() requires(N != Dynamic) { }
explicit ScalarField(int n) requires (N == Dynamic) : Base(n) { }
explicit ScalarField(const FieldType& f) : f_(f) {};

// assignement and constructor from a ScalarExpr requires the base type F to be a std::function<>
template <
typename E, typename U = FieldType,
typename std::enable_if<std::is_same<U, std::function<double(VectorType)>>::value, int>::type = 0>
ScalarField(const ScalarExpr<N, E>& f) {
template <typename E>
ScalarField(const ScalarExpr<N, E>& f)
requires(std::is_same<FieldType, std::function<double(VectorType)>>::value) {
E op = f.get();
std::function<double(VectorType)> field_expr = [op](SVector<N> x) -> double { return op(x); };
f_ = field_expr;
f_ = [op](SVector<N> x) -> double { return op(x); };
}
template <typename E, typename U = FieldType>
typename std::enable_if<std::is_same<U, std::function<double(VectorType)>>::value, ScalarField<N>&>::type
operator=(const ScalarExpr<N, E>& f) {
template <typename E>
ScalarField& operator=(const ScalarExpr<N, E>& f)
requires(std::is_same<FieldType, std::function<double(VectorType)>>::value) {
E op = f.get();
std::function<double(VectorType)> field_expr = [op](VectorType x) -> double { return op(x); };
f_ = field_expr;
f_ = [op](VectorType x) -> double { return op(x); };
return *this;
}
// assignment from lambda expression
template <typename L, typename U = FieldType>
typename std::enable_if<std::is_same<U, std::function<double(VectorType)>>::value, ScalarField<N>&>::type
operator=(const L& lambda) {
template <typename L>
ScalarField& operator=(const L& lambda)
requires(std::is_same<FieldType, std::function<double(VectorType)>>::value) {
f_ = lambda;
return *this;
}
Expand Down
Loading

0 comments on commit e3f061d

Please sign in to comment.