Skip to content

Commit

Permalink
space-time PDEs cannot be instantiated without time domain, type eras…
Browse files Browse the repository at this point in the history 
…ed PDEs expose also time_domain() and initial_condition()
  • Loading branch information
@AlePalu
AlePalu committed Jan 13, 2024
1 parent 189ee75 commit 182556d
Show file tree
Hide file tree
Showing 3 changed files with 60 additions and 49 deletions.
2 changes: 1 addition & 1 deletion fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ class FEMLinearParabolicSolver : public FEMSolverBase<D, E, F, Ts...> {
if (!this->is_init) throw std::runtime_error("solver must be initialized first!");

Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
this->set_deltaT((pde.time()[1] - pde.time()[0]));
this->set_deltaT((pde.time_domain()[1] - pde.time_domain()[0]));
std::size_t n = this->n_dofs(); // degrees of freedom in space
std::size_t m = pde.forcing_data().cols(); // number of iterations for time loop

Expand Down
7 changes: 5 additions & 2 deletions fdaPDE/pde/differential_operators.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,13 +40,16 @@ FDAPDE_DEFINE_DIFFERENTIAL_OPERATOR(BiLaplacian, bilaplacian);
template <typename E> struct is_symmetric {
static constexpr bool value = std::decay<E>::type::is_symmetric;
};

// trait to detect if the differential operator denotes a parabolic problem.
template <typename E_> struct is_parabolic {
typedef typename std::decay<E_>::type E;
// returns true if the time derivative operator dT() is detected in the expression
static constexpr bool value = has_instance_of<dT, decltype(std::declval<E>().get_operator_type())>::value;
};
// detects if the problem is stationary (does not involve time)
template <typename E_> struct is_stationary {
typedef typename std::decay<E_>::type E;
static constexpr bool value = !has_instance_of<dT, decltype(std::declval<E>().get_operator_type())>::value;
};

} // namespace core
} // namespace fdapde
Expand Down
100 changes: 54 additions & 46 deletions fdaPDE/pde/pde.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,48 +40,53 @@ template <typename D, // problem domain
typename... Ts> // parameters forwarded to S
class PDE {
public:
typedef D DomainType; // triangulated domain
static constexpr int M = DomainType::local_dimension;
static constexpr int N = DomainType::embedding_dimension;
typedef E OperatorType; // differential operator in its strong-formulation
using SpaceDomainType = D; // triangulated spatial domain
using TimeDomainType = DVector<double>; // time-interval [0,T] (for time-dependent PDEs)
static constexpr int M = SpaceDomainType::local_dimension;
static constexpr int N = SpaceDomainType::embedding_dimension;
using OperatorType = E; // differential operator in its strong-formulation
static_assert(
std::is_base_of<DifferentialExpr<OperatorType>, OperatorType>::value, "E is not a valid differential operator");
typedef F ForcingType; // type of forcing object (either a matrix or a callable object)
std::is_base_of<DifferentialExpr<OperatorType>, OperatorType>::value);
using ForcingType = F; // type of forcing object (either a matrix or a callable object)
static_assert(
std::is_same<DMatrix<double>, F>::value || std::is_base_of<ScalarExpr<D::embedding_dimension, F>, F>::value,
"forcing is not a matrix or a scalar expression || N != F::base");
typedef typename pde_solver_selector<S, D, E, F, Ts...>::type SolverType;
typedef typename SolverType::FunctionalBasis FunctionalBasis; // function space approximating the solution space
typedef typename SolverType::Quadrature Quadrature; // quadrature for numerical integral approximations
std::is_same<DMatrix<double>, ForcingType>::value ||
std::is_base_of<ScalarExpr<SpaceDomainType::embedding_dimension, ForcingType>, ForcingType>::value);
using SolverType = typename pde_solver_selector<S, SpaceDomainType, OperatorType, ForcingType, Ts...>::type;
using FunctionalBasis = typename SolverType::FunctionalBasis; // function space approximating the solution space
using Quadrature = typename SolverType::Quadrature; // quadrature for numerical integral approximations

// minimal constructor, use below setters to complete the construction of a PDE object
PDE(const D& domain) : domain_(domain), solver_(domain) { }
PDE(const D& domain, const DVector<double>& time) : domain_(domain), time_(time), solver_(domain) { }
PDE(const D& domain, E diff_op) : domain_(domain), diff_op_(diff_op), solver_(domain) { }
fdapde_enable_constructor_if(is_parabolic, E) PDE(const D& domain, const DVector<double>& time, E diff_op) :
domain_(domain), time_(time), diff_op_(diff_op), solver_(domain) { }
void set_forcing(const F& forcing_data) { forcing_data_ = forcing_data; }
void set_differential_operator(E diff_op) { diff_op_ = diff_op; }
// full constructors
PDE(const D& domain, E diff_op, const F& forcing_data) :
// space-only constructors
fdapde_enable_constructor_if(is_stationary, OperatorType) PDE(const D& domain) :
domain_(domain), solver_(domain) { }
fdapde_enable_constructor_if(is_stationary, OperatorType) PDE(const D& domain, E diff_op) :
domain_(domain), diff_op_(diff_op), solver_(domain) { }
fdapde_enable_constructor_if(is_stationary, OperatorType)
PDE(const SpaceDomainType& domain, OperatorType diff_op, const ForcingType& forcing_data) :
domain_(domain), diff_op_(diff_op), forcing_data_(forcing_data), solver_(domain) { }
fdapde_enable_constructor_if(is_parabolic, E)
PDE(const D& domain, const DVector<double>& time, E diff_op, const F& forcing_data) :
domain_(domain), time_(time), diff_op_(diff_op), forcing_data_(forcing_data), solver_(domain) { }
// space-time constructors
fdapde_enable_constructor_if(is_parabolic, OperatorType) PDE(const D& domain, const TimeDomainType& t) :
domain_(domain), time_domain_(t), solver_(domain) { }
fdapde_enable_constructor_if(is_parabolic, OperatorType) PDE(const D& domain, const TimeDomainType& t, E diff_op) :
domain_(domain), time_domain_(t), diff_op_(diff_op), solver_(domain) { }
fdapde_enable_constructor_if(is_parabolic, OperatorType) PDE(
const SpaceDomainType& domain, const TimeDomainType& t, OperatorType diff_op, const ForcingType& forcing_data) :
domain_(domain), time_domain_(t), diff_op_(diff_op), forcing_data_(forcing_data), solver_(domain) { }

// setters
void set_forcing(const ForcingType& forcing_data) { forcing_data_ = forcing_data; }
void set_differential_operator(OperatorType diff_op) { diff_op_ = diff_op; }
void set_dirichlet_bc(const DMatrix<double>& data) { boundary_data_ = data; }
void set_initial_condition(const DVector<double>& data) { initial_condition_ = data; };

// getters
const DomainType& domain() const { return domain_; }
const DVector<double>& time() const {return time_;}
const SpaceDomainType& domain() const { return domain_; }
const DVector<double>& time_domain() const { return time_domain_; }
OperatorType differential_operator() const { return diff_op_; }
const ForcingType& forcing_data() const { return forcing_data_; }
const DVector<double>& initial_condition() const { return initial_condition_; }
const DMatrix<double>& boundary_data() const { return boundary_data_; };
const Quadrature& integrator() const { return solver_.integrator(); }
const FunctionalBasis& basis() const { return solver_.basis(); }

// evaluates the functional basis defined over the pyhisical domain on a given set of locations
template <template <typename> typename EvaluationPolicy>
std::pair<SpMatrix<double>, DVector<double>> eval_functional_basis(const DMatrix<double>& locs) const {
return solver_.basis().template eval<EvaluationPolicy>(locs);
Expand All @@ -99,52 +104,55 @@ class PDE {
solver_.solve(*this);
}
private:
const DomainType& domain_; // triangulated problem domain
const DVector<double> time_;
const SpaceDomainType& domain_; // triangulated spatial domain
const TimeDomainType time_domain_; // time interval [0, T], for space-time PDEs
OperatorType diff_op_; // differential operator in its strong formulation
ForcingType forcing_data_; // forcing data
DVector<double> initial_condition_ {}; // initial condition, (for space-time problems only)
DVector<double> initial_condition_ {}; // initial condition, for space-time PDEs
SolverType solver_ {}; // problem solver
DMatrix<double> boundary_data_; // boundary conditions
};

// type-erasure wrapper (forcing type must be convertible to a DMatrix<double>)
// type-erasure wrapper (we require ForcingType_ to be convertible to a DMatrix<double>)
struct I_PDE {
void init() { invoke<void, 0>(*this); }
void solve() { invoke<void, 1>(*this); }
// getters
decltype(auto) solution() const { return invoke<const DMatrix<double>& , 2>(*this); }
decltype(auto) force() const { return invoke<const DMatrix<double>& , 3>(*this); }
decltype(auto) stiff() const { return invoke<const SpMatrix<double>&, 4>(*this); }
decltype(auto) mass() const { return invoke<const SpMatrix<double>&, 5>(*this); }
decltype(auto) quadrature_nodes() const { return invoke<DMatrix<double> , 6>(*this); }
decltype(auto) n_dofs() const { return invoke<std::size_t , 7>(*this); }
decltype(auto) dof_coords() const { return invoke<DMatrix<double> , 8>(*this); }
decltype(auto) forcing_data() const { return invoke<const DMatrix<double>& , 9>(*this); }
decltype(auto) solution() const { return invoke<const DMatrix<double>& , 2>(*this); }
decltype(auto) force() const { return invoke<const DMatrix<double>& , 3>(*this); }
decltype(auto) stiff() const { return invoke<const SpMatrix<double>&, 4>(*this); }
decltype(auto) mass() const { return invoke<const SpMatrix<double>&, 5>(*this); }
decltype(auto) quadrature_nodes() const { return invoke<DMatrix<double> , 6>(*this); }
decltype(auto) n_dofs() const { return invoke<std::size_t , 7>(*this); }
decltype(auto) dof_coords() const { return invoke<DMatrix<double> , 8>(*this); }
decltype(auto) forcing_data() const { return invoke<const DMatrix<double>& , 9>(*this); }
decltype(auto) time_domain() const { return invoke<const DVector<double>&, 10>(*this); }
decltype(auto) initial_condition() const { return invoke<const DVector<double>&, 11>(*this); }

struct eval_basis_ret_type { SpMatrix<double> Psi; DVector<double> D; };
std::optional<eval_basis_ret_type> eval_basis(eval e, const DMatrix<double>& locs) const {
using RetType = eval_basis_ret_type;
switch (e) { // run-time switch based on sampling strategy
case eval::pointwise:
return std::optional<RetType>{invoke<RetType, 10>(*this, locs)};
return std::optional<RetType>{invoke<RetType, 12>(*this, locs)};
case eval::areal:
return std::optional<RetType>{invoke<RetType, 11>(*this, locs)};
return std::optional<RetType>{invoke<RetType, 13>(*this, locs)};
}
return std::nullopt;
}
// setters
template <typename F> void set_forcing(const F& data) { fdapde::invoke<void, 12>(*this, DMatrix<double>(data)); }
void set_dirichlet_bc(const DMatrix<double>& data) { fdapde::invoke<void, 13>(*this, data); }
void set_initial_condition(const DVector<double>& data) { fdapde::invoke<void, 14>(*this, data); }
template <typename E> void set_differential_operator(E diff_op) { fdapde::invoke<void, 15>(*this, diff_op); }
template <typename F> void set_forcing(const F& data) { fdapde::invoke<void, 14>(*this, DMatrix<double>(data)); }
void set_dirichlet_bc(const DMatrix<double>& data) { fdapde::invoke<void, 15>(*this, data); }
void set_initial_condition(const DVector<double>& data) { fdapde::invoke<void, 16>(*this, data); }
template <typename E> void set_differential_operator(E diff_op) { fdapde::invoke<void, 17>(*this, diff_op); }

// function pointers forwardings
template <typename T>
using fn_ptrs = fdapde::mem_fn_ptrs<
&T::init, &T::solve, // initialization and pde solution
// getters
&T::solution, &T::force, &T::stiff, &T::mass, &T::quadrature_nodes, &T::n_dofs, &T::dof_coords, &T::forcing_data,
&T::time_domain, &T::initial_condition,
&T::template eval_functional_basis<pointwise_evaluation>, &T::template eval_functional_basis<areal_evaluation>,
// setters
&T::set_forcing, &T::set_dirichlet_bc, &T::set_initial_condition, &T::set_differential_operator>;
Expand Down

0 comments on commit 182556d

Please sign in to comment.