Merge branch 'develop' of https://github.com/fdaPDE/fdaPDE-core into …

…develop

fdaPDE/finite_elements/solvers/fem_linear_parabolic_solver.h

@@ -40,6 +40,7 @@ class FEMLinearParabolicSolver : public FEMSolverBase<D, E, F, Ts...> {
         if (!this->is_init) throw std::runtime_error("solver must be initialized first!");
         // define eigen system solver, use SparseLU decomposition.
         Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> 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
 

fdaPDE/pde/pde.h

@@ -25,6 +25,7 @@
 #include "../utils/integration/integrator.h"
 #include "differential_expressions.h"
 #include "../utils/type_erasure.h"
+#include "differential_operators.h"
 
 namespace fdapde {
 namespace core {
@@ -57,19 +58,25 @@ class PDE {
 
     // 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<double>& time) : domain_(domain), time_(time) { };
     PDE(const D& domain, E diff_op) : domain_(domain), diff_op_(diff_op) { };
+    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) {};
     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) { }
-
+    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) { }
     // setters
     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_;}
     OperatorType differential_operator() const { return diff_op_; }
     const ForcingType& forcing_data() const { return forcing_data_; }
     const DVector<double>& initial_condition() const { return initial_condition_; }
@@ -95,6 +102,7 @@ class PDE {
     }
    private:
     const DomainType& domain_;               // triangulated problem domain
+    const DVector<double> time_;
     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)

fdaPDE/utils/traits.h

@@ -177,6 +177,10 @@ struct fn_ptr_traits_impl<R (T::*)(Args...) const> : public fn_ptr_traits_base<R
     using MemFnPtrType = R (T::*)(Args...) const;
 };
 template <auto FnPtr> struct fn_ptr_traits : public fn_ptr_traits_impl<decltype(FnPtr)> { };
+
+// macro for SFINAE based selection of constructor in overload resolution
+#define fdapde_enable_constructor_if(TRAIT_CONDITION, T)                                                               \
+    template <typename T_ = T, typename std::enable_if<TRAIT_CONDITION<T_>::value, int>::type = 0>
 
 }   // namespace fdapde
 

fdaPDE/utils/type_erasure.h

@@ -20,6 +20,7 @@
 #include <typeinfo>
 #include <typeindex>
 #include <unordered_map>
+#include <utility>
 #include "assert.h"
 #include "traits.h"
 

test/src/fem_pde_test.cpp

@@ -245,7 +245,7 @@ TEST(fem_pde_test, parabolic_isotropic_order2) {
 
     MeshLoader<Mesh2D> unit_square("unit_square");
     auto L = dt<FEM>() -laplacian<FEM>();
-    PDE<decltype(unit_square.mesh), decltype(L), DMatrix<double>, FEM, fem_order<2>> pde_(unit_square.mesh);
+    PDE<decltype(unit_square.mesh), decltype(L), DMatrix<double>, FEM, fem_order<2>> pde_(unit_square.mesh, times);
     pde_.set_differential_operator(L);
 
     // compute boundary condition and exact solution
@@ -304,7 +304,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<double> times(M,1);
     double time_max = 1.;
     for(int j = 0; j < M; ++j){
@@ -324,7 +324,7 @@ TEST(fem_pde_test, parabolic_isotropic_order1_convergence) {
       std::string domain_name = "unit_square_" + std::to_string(N(n));    
       MeshLoader<Mesh2D> unit_square(domain_name);
       auto L = dt<FEM>() -laplacian<FEM>();
-      PDE<decltype(unit_square.mesh), decltype(L), DMatrix<double>, FEM, fem_order<1>> pde_(unit_square.mesh);
+      PDE<decltype(unit_square.mesh), decltype(L), DMatrix<double>, FEM, fem_order<1>> pde_(unit_square.mesh, times);
       pde_.set_differential_operator(L);
 
       // compute boundary condition and exact solution