From d48c5bfd8b6dc1d1b29bd8b329cbacb67044c00b Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Mon, 20 Nov 2023 11:43:30 +0000 Subject: [PATCH 1/3] PDE sets temporal mesh. --- .../finite_elements/solvers/fem_linear_parabolic_solver.h | 1 + fdaPDE/pde/pde.h | 4 ++++ test/src/fem_pde_test.cpp | 6 +++--- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h b/fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h index f17120bb..03889117 100644 --- a/fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h +++ b/fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h @@ -40,6 +40,7 @@ class FEMLinearParabolicSolver : public FEMSolverBase { if (!this->is_init) throw std::runtime_error("solver must be initialized first!"); // define eigen system solver, use SparseLU decomposition. Eigen::SparseLU, Eigen::COLAMDOrdering> solver; + this->set_deltaT((pde.time()[1]-pde.time()[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 diff --git a/fdaPDE/pde/pde.h b/fdaPDE/pde/pde.h index 139bb0be..b90b6658 100644 --- a/fdaPDE/pde/pde.h +++ b/fdaPDE/pde/pde.h @@ -73,7 +73,9 @@ class PDE : public PDEBase { // minimal constructor, use below setters to complete the construction of a PDE object PDE(const D& domain) : domain_(domain) { } + PDE(const D& domain, const DVector& time) : domain_(domain), time_(time) { }; PDE(const D& domain, E diff_op) : domain_(domain), diff_op_(diff_op) { }; + PDE(const D& domain, E diff_op, const DVector& time) : domain_(domain), diff_op_(diff_op), time_(time) { }; void set_forcing(const F& forcing_data) { forcing_data_ = forcing_data; } void set_differential_operator(E diff_op) { diff_op_ = diff_op; } // full constructors @@ -86,6 +88,7 @@ class PDE : public PDEBase { // getters const DomainType& domain() const { return domain_; } + const DVector& time() const {return time_;} OperatorType differential_operator() const { return diff_op_; } const ForcingType& forcing_data() const { return forcing_data_; } const DVector& initial_condition() const { return initial_condition_; } @@ -110,6 +113,7 @@ class PDE : public PDEBase { private: const DomainType& domain_; // triangulated problem domain + const DVector time_; OperatorType diff_op_; // differential operator in its strong formulation ForcingType forcing_data_; // forcing data DVector initial_condition_ {}; // initial condition, (for space-time problems only) diff --git a/test/src/fem_pde_test.cpp b/test/src/fem_pde_test.cpp index a434b4a4..d50efbda 100644 --- a/test/src/fem_pde_test.cpp +++ b/test/src/fem_pde_test.cpp @@ -244,7 +244,7 @@ TEST(fem_pde_test, parabolic_isotropic_order2) { MeshLoader unit_square("unit_square"); auto L = dt() -laplacian(); - PDE, FEM, fem_order<2>> pde_(unit_square.mesh); + PDE, FEM, fem_order<2>> pde_(unit_square.mesh, times); pde_.set_differential_operator(L); // compute boundary condition and exact solution @@ -303,7 +303,7 @@ TEST(fem_pde_test, parabolic_isotropic_order2) { TEST(fem_pde_test, parabolic_isotropic_order1_convergence) { // exact solution constexpr double pi = 3.14159265358979323846; - int M = 101; + int M = 31; DMatrix times(M,1); double time_max = 1.; for(int j = 0; j < M; ++j){ @@ -323,7 +323,7 @@ TEST(fem_pde_test, parabolic_isotropic_order1_convergence) { std::string domain_name = "unit_square_" + std::to_string(N(n)); MeshLoader unit_square(domain_name); auto L = dt() -laplacian(); - PDE, FEM, fem_order<1>> pde_(unit_square.mesh); + PDE, FEM, fem_order<1>> pde_(unit_square.mesh, times); pde_.set_differential_operator(L); // compute boundary condition and exact solution From 106f7110c3f95595a351683a665dc473f5fdedc8 Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Mon, 20 Nov 2023 14:30:18 +0000 Subject: [PATCH 2/3] Added full constructor for Parabolic PDEs. --- fdaPDE/pde/pde.h | 5 +++-- fdaPDE/utils/type_erasure.h | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/fdaPDE/pde/pde.h b/fdaPDE/pde/pde.h index b90b6658..579d9859 100644 --- a/fdaPDE/pde/pde.h +++ b/fdaPDE/pde/pde.h @@ -75,13 +75,14 @@ class PDE : public PDEBase { PDE(const D& domain) : domain_(domain) { } PDE(const D& domain, const DVector& time) : domain_(domain), time_(time) { }; PDE(const D& domain, E diff_op) : domain_(domain), diff_op_(diff_op) { }; - PDE(const D& domain, E diff_op, const DVector& time) : domain_(domain), diff_op_(diff_op), time_(time) { }; + PDE(const D& domain, const DVector& time, E diff_op) : domain_(domain), time_(time), diff_op_(diff_op) { }; 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) : domain_(domain), diff_op_(diff_op), forcing_data_(forcing_data) { } - + PDE(const D& domain, const DVector& time, E diff_op, const F& forcing_data) : + domain_(domain), time_(time), diff_op_(diff_op), forcing_data_(forcing_data) { } // setters virtual void set_dirichlet_bc(const DMatrix& data) { boundary_data_ = data; } virtual void set_initial_condition(const DVector& data) { initial_condition_ = data; }; diff --git a/fdaPDE/utils/type_erasure.h b/fdaPDE/utils/type_erasure.h index 6438a52f..62160833 100644 --- a/fdaPDE/utils/type_erasure.h +++ b/fdaPDE/utils/type_erasure.h @@ -16,7 +16,7 @@ #ifndef __TYPE_ERASURE_H__ #define __TYPE_ERASURE_H__ - +#include #include "assert.h" #include "traits.h" From 7605b1c8f58ef02d109bd41424d9a2184ebdc0cb Mon Sep 17 00:00:00 2001 From: AlePalu Date: Mon, 20 Nov 2023 16:43:11 +0100 Subject: [PATCH 3/3] added macro for SFINAE selection of constructor in overload resolution set, time-dependent pde constructor enabled only for time-dependent PDEs --- fdaPDE/pde/pde.h | 7 +++++-- fdaPDE/utils/traits.h | 4 ++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/fdaPDE/pde/pde.h b/fdaPDE/pde/pde.h index 579d9859..30c432a1 100644 --- a/fdaPDE/pde/pde.h +++ b/fdaPDE/pde/pde.h @@ -24,6 +24,7 @@ #include "../utils/symbols.h" #include "../utils/integration/integrator.h" #include "differential_expressions.h" +#include "differential_operators.h" namespace fdapde { namespace core { @@ -75,13 +76,15 @@ class PDE : public PDEBase { PDE(const D& domain) : domain_(domain) { } PDE(const D& domain, const DVector& time) : domain_(domain), time_(time) { }; PDE(const D& domain, E diff_op) : domain_(domain), diff_op_(diff_op) { }; - PDE(const D& domain, const DVector& time, E diff_op) : domain_(domain), time_(time), diff_op_(diff_op) { }; + fdapde_enable_constructor_if(is_parabolic, E) PDE(const D& domain, const DVector& time, E diff_op) : + domain_(domain), time_(time), diff_op_(diff_op) {}; 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) : domain_(domain), diff_op_(diff_op), forcing_data_(forcing_data) { } - PDE(const D& domain, const DVector& time, E diff_op, const F& forcing_data) : + fdapde_enable_constructor_if(is_parabolic, E) + PDE(const D& domain, const DVector& time, E diff_op, const F& forcing_data) : domain_(domain), time_(time), diff_op_(diff_op), forcing_data_(forcing_data) { } // setters virtual void set_dirichlet_bc(const DMatrix& data) { boundary_data_ = data; } diff --git a/fdaPDE/utils/traits.h b/fdaPDE/utils/traits.h index 131dfc20..f364f0a8 100644 --- a/fdaPDE/utils/traits.h +++ b/fdaPDE/utils/traits.h @@ -177,6 +177,10 @@ struct fn_ptr_traits_impl : public fn_ptr_traits_base struct fn_ptr_traits : public fn_ptr_traits_impl { }; + +// macro for SFINAE based selection of constructor in overload resolution +#define fdapde_enable_constructor_if(TRAIT_CONDITION, T) \ + template ::value, int>::type = 0> } // namespace fdapde