diff --git a/fdaPDE/assembly.h b/fdaPDE/assembly.h index d6b37b7..cc6ff81 100644 --- a/fdaPDE/assembly.h +++ b/fdaPDE/assembly.h @@ -24,6 +24,15 @@ namespace fdapde { [[maybe_unused]] static constexpr int CellMajor = 0; [[maybe_unused]] static constexpr int FaceMajor = 1; +// test function forward decl +template struct TestFunction; +template // deduction guide +TestFunction(FunctionSpace function_space) -> TestFunction; +// trial function forward decl +template struct TrialFunction; +template // deduction guide +TrialFunction(FunctionSpace function_space) -> TrialFunction; + namespace internals { // set of internal utilities to write weak form assembly loops @@ -34,7 +43,7 @@ template constexpr decltype(auto) trial_space(Xpr&& xpr) { decltype([]() { return requires { typename Xpr_::TrialSpace; }; }), std::decay_t>(); fdapde_static_assert(found, NO_TRIAL_SPACE_FOUND_IN_EXPRESSION); return meta::xpr_query< - decltype([](Xpr_&& xpr) -> auto& { return xpr.fe_space(); }), + decltype([](Xpr_&& xpr) -> auto& { return xpr.function_space(); }), decltype([]() { return requires { typename Xpr_::TrialSpace; }; })>(std::forward(xpr)); } template using trial_space_t = std::decay_t()))>; @@ -44,7 +53,7 @@ template constexpr decltype(auto) test_space(Xpr&& xpr) { decltype([]() { return requires { typename Xpr_::TestSpace; }; }), std::decay_t>(); fdapde_static_assert(found, NO_TEST_SPACE_FOUND_IN_EXPRESSION); return meta::xpr_query< - decltype([](Xpr_&& xpr) -> auto& { return xpr.fe_space(); }), + decltype([](Xpr_&& xpr) -> auto& { return xpr.function_space(); }), decltype([]() { return requires { typename Xpr_::TestSpace; }; })>(std::forward(xpr)); } template using test_space_t = std::decay_t()))>; diff --git a/fdaPDE/fields/spline.h b/fdaPDE/fields/spline.h index f3c776c..e769939 100644 --- a/fdaPDE/fields/spline.h +++ b/fdaPDE/fields/spline.h @@ -27,7 +27,7 @@ class Spline : public ScalarBase<1, Spline> { public: using Base = ScalarBase<1, Spline>; static constexpr int StaticInputSize = 1; - static constexpr int NestAsRef = 0; // avoid nesting as reference, .derive() generates temporaries + static constexpr int NestAsRef = 0; static constexpr int XprBits = 0; static constexpr int Order = Dynamic; using Scalar = double; @@ -124,12 +124,17 @@ class Spline : public ScalarBase<1, Spline> { requires(fdapde::is_subscriptable) constexpr Scalar operator()(const InputType_& p_) const { auto p = p_[0]; + std::cout << "order: " << order_ << std::endl; + for(auto xyz : knots_) std::cout << xyz << ", "; + std::cout << std::endl; int m = knots_.size() - 1; std::vector N(order_ + 1, 0.0); // special cases if ((i_ == 0 && p == knots_[0]) || (i_ == m - order_ - 1 && p == knots_[m])) return 1.0; // local property: return 0 if p is outside the range of this basis function - if (p < knots_[i_] || p >= knots_[i_ + order_ + 1]) return 0.0; + std::cout << "i_ " << i_ << "/" << knots_.size() << " : " << knots_[i_] << ", " << p << ", " << knots_[i_ + order_ + 1] << std::endl; + //if (p < knots_[i_] || p >= knots_[i_ + order_ + 1]) return 0.0; + std::cout << ":)" << std::endl; // initialize 0th degree basis functions for (int j = 0; j <= order_; ++j) { if (p >= knots_[i_ + j] && p < knots_[i_ + j + 1]) { N[j] = 1.0; } @@ -150,6 +155,7 @@ class Spline : public ScalarBase<1, Spline> { } } } + std::cout << N[0] << std::endl; return N[0]; } constexpr Derivative gradient(int n = 1) const { return Derivative(knots_, i_, order_, n); } diff --git a/fdaPDE/finite_elements.h b/fdaPDE/finite_elements.h index 7843751..09570b4 100644 --- a/fdaPDE/finite_elements.h +++ b/fdaPDE/finite_elements.h @@ -21,7 +21,6 @@ #include "finite_elements/dof_segment.h" #include "finite_elements/dof_tetrahedron.h" #include "finite_elements/dof_triangle.h" -#include "finite_elements/fe_assembler.h" #include "finite_elements/fe_bilinear_form_assembler.h" #include "finite_elements/fe_function.h" #include "finite_elements/fe_integration.h" diff --git a/fdaPDE/finite_elements/fe_assembler_base.h b/fdaPDE/finite_elements/fe_assembler_base.h index 85cacb5..2e65a1d 100644 --- a/fdaPDE/finite_elements/fe_assembler_base.h +++ b/fdaPDE/finite_elements/fe_assembler_base.h @@ -29,7 +29,7 @@ namespace fdapde { template struct FeMap; -enum fe_assembler_flags { +enum class fe_assembler_flags { compute_shape_values = 0x0001, compute_shape_grad = 0x0002, compute_second_derivatives = 0x0004, diff --git a/fdaPDE/finite_elements/fe_bilinear_form_assembler.h b/fdaPDE/finite_elements/fe_bilinear_form_assembler.h index 87e29c8..6407d09 100644 --- a/fdaPDE/finite_elements/fe_bilinear_form_assembler.h +++ b/fdaPDE/finite_elements/fe_bilinear_form_assembler.h @@ -114,7 +114,7 @@ class fe_bilinear_form_assembly_loop : cexpr::Matrix trial_divs; std::unordered_map> fe_map_buff; - if constexpr (Form::XprBits & fe_assembler_flags::compute_physical_quad_nodes) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_physical_quad_nodes)) { Base::distribute_quadrature_nodes( fe_map_buff, begin, end); // distribute quadrature nodes on physical mesh (if required) } @@ -124,14 +124,14 @@ class fe_bilinear_form_assembly_loop : for (iterator it = begin; it != end; ++it) { // update fe_packet content based on form requests fe_packet.cell_measure = it->measure(); - if constexpr (Form::XprBits & fe_assembler_flags::compute_cell_diameter) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_cell_diameter)) { fe_packet.cell_diameter = it->diameter(); } - if constexpr (Form::XprBits & fe_assembler_flags::compute_shape_grad) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_shape_grad)) { Base::eval_shape_grads_on_cell(it, test_shape_grads_, test_grads); if constexpr (is_petrov_galerkin) Base::eval_shape_grads_on_cell(it, trial_shape_grads_, trial_grads); } - if constexpr (Form::XprBits & fe_assembler_flags::compute_shape_div) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_shape_div)) { // divergence is defined only for vector elements, skeep computation for scalar element case if constexpr (n_test_components != 1) Base::eval_shape_div_on_cell(it, test_shape_grads_, test_divs); if constexpr (is_petrov_galerkin && n_trial_components != 1) @@ -145,24 +145,24 @@ class fe_bilinear_form_assembly_loop : for (int j = 0; j < n_test_basis; ++j) { // test function loop double value = 0; for (int q_k = 0; q_k < n_quadrature_nodes; ++q_k) { - if constexpr (Form::XprBits & fe_assembler_flags::compute_shape_values) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_shape_values)) { fe_packet.trial_value.assign_inplace_from(trial_shape_values_.template slice<0, 1>(i, q_k)); fe_packet.test_value .assign_inplace_from(test_shape_values_ .template slice<0, 1>(j, q_k)); } - if constexpr (Form::XprBits & fe_assembler_flags::compute_shape_grad) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_shape_grad)) { fe_packet.trial_grad.assign_inplace_from(is_galerkin ? test_grads.template slice<0, 1>(i, q_k) : trial_grads.template slice<0, 1>(i, q_k)); fe_packet.test_grad .assign_inplace_from(test_grads.template slice<0, 1>(j, q_k)); } - if constexpr (Form::XprBits & fe_assembler_flags::compute_shape_div) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_shape_div)) { if constexpr (n_trial_components != 1) { fe_packet.trial_div = (is_galerkin && n_test_components != 1) ? test_divs(i, q_k) : trial_divs(i, q_k); } if constexpr (n_test_components != 1) fe_packet.test_div = test_divs(j, q_k); } - if constexpr (Form::XprBits & fe_assembler_flags::compute_physical_quad_nodes) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_physical_quad_nodes)) { fe_packet.quad_node_id = local_cell_id * n_quadrature_nodes + q_k; } value += Quadrature::weights[q_k] * form_(fe_packet); diff --git a/fdaPDE/finite_elements/fe_function.h b/fdaPDE/finite_elements/fe_function.h index 195934c..ba3e2c0 100644 --- a/fdaPDE/finite_elements/fe_function.h +++ b/fdaPDE/finite_elements/fe_function.h @@ -18,36 +18,33 @@ #define __FE_FUNCTION_H__ #include "../fields/scalar_field.h" -#include "fe_assembler.h" #include "fe_mass_assembler.h" -namespace fdapde { - -template struct TestFunction; -template struct TrialFunction; - +namespace fdapde { namespace internals { template -struct scalar_test_function_impl : public fdapde::ScalarBase> { +struct fe_scalar_test_function_impl : + public fdapde::ScalarBase> { using TestSpace = std::decay_t; - using Base = fdapde::ScalarBase>; + using Base = fdapde::ScalarBase>; using InputType = internals::fe_assembler_packet; using Scalar = double; static constexpr int StaticInputSize = TestSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_values; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_values); - constexpr scalar_test_function_impl() noexcept = default; - constexpr scalar_test_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } + constexpr fe_scalar_test_function_impl() noexcept = default; + constexpr fe_scalar_test_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } struct FirstPartialDerivative : fdapde::ScalarBase { - using Derived = TestFunction; + using TestSpace = std::decay_t; + using Derived = TestFunction; using InputType = internals::fe_assembler_packet; using Scalar = double; static constexpr int StaticInputSize = TestSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_grad; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_grad); FirstPartialDerivative() noexcept = default; FirstPartialDerivative([[maybe_unused]] const Derived& f, int i) noexcept : i_(i) { } @@ -59,31 +56,32 @@ struct scalar_test_function_impl : public fdapde::ScalarBase -struct vector_test_function_impl : public fdapde::MatrixBase> { +struct fe_vector_test_function_impl : + public fdapde::MatrixBase> { using TestSpace = std::decay_t; - using Base = fdapde::MatrixBase>; + using Base = fdapde::MatrixBase>; using InputType = internals::fe_assembler_packet; using Scalar = double; static constexpr int StaticInputSize = TestSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_values; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_values); static constexpr int ReadOnly = 1; static constexpr int Rows = FeSpace_::n_components; static constexpr int Cols = 1; - constexpr vector_test_function_impl() noexcept = default; - constexpr vector_test_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } + constexpr fe_vector_test_function_impl() noexcept = default; + constexpr fe_vector_test_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } template struct Jacobian : fdapde::MatrixBase> { - // using TestSpace = std::decay_t; + using TestSpace = std::decay_t; using Derived = Derived_; template using Meta = Jacobian; using InputType = internals::fe_assembler_packet; @@ -92,7 +90,7 @@ struct vector_test_function_impl : public fdapde::MatrixBase::type xpr_; }; template struct Divergence : fdapde::ScalarBase> { - // using TestSpace = std::decay_t; + using TestSpace = std::decay_t; using Derived = Derived_; template using Meta = Divergence; using InputType = internals::fe_assembler_packet; using Scalar = typename Derived::Scalar; static constexpr int StaticInputSize = TestSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_div; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_div); explicit constexpr Divergence(const Derived& xpr) noexcept : xpr_(xpr) { } // fe assembly evaluation @@ -132,32 +130,34 @@ struct vector_test_function_impl : public fdapde::MatrixBase -struct scalar_trial_function_impl : public fdapde::ScalarBase> { +struct fe_scalar_trial_function_impl : + public fdapde::ScalarBase> { using TrialSpace = std::decay_t; - using Base = fdapde::ScalarBase>; + using Base = fdapde::ScalarBase>; using InputType = internals::fe_assembler_packet; using Scalar = double; static constexpr int StaticInputSize = TrialSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_values; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_values); - constexpr scalar_trial_function_impl() noexcept = default; - constexpr scalar_trial_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } + constexpr fe_scalar_trial_function_impl() noexcept = default; + constexpr fe_scalar_trial_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } struct FirstPartialDerivative : fdapde::ScalarBase { - using Derived = TrialFunction; + using TrialSpace = std::decay_t; + using Derived = TrialFunction; using InputType = internals::fe_assembler_packet; using Scalar = double; static constexpr int StaticInputSize = TrialSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_grad; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_grad); FirstPartialDerivative() noexcept = default; FirstPartialDerivative([[maybe_unused]] const Derived& f, int i) noexcept : i_(i) { } @@ -169,30 +169,32 @@ struct scalar_trial_function_impl : public fdapde::ScalarBase -struct vector_trial_function_impl : public fdapde::MatrixBase> { +struct fe_vector_trial_function_impl : + public fdapde::MatrixBase> { using TrialSpace = std::decay_t; - using Base = fdapde::MatrixBase>; + using Base = fdapde::MatrixBase>; using InputType = internals::fe_assembler_packet; using Scalar = double; static constexpr int StaticInputSize = TrialSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_values; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_values); static constexpr int ReadOnly = 1; static constexpr int Rows = FeSpace_::n_components; static constexpr int Cols = 1; - constexpr vector_trial_function_impl() noexcept = default; - constexpr vector_trial_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } + constexpr fe_vector_trial_function_impl() noexcept = default; + constexpr fe_vector_trial_function_impl(FeSpace_& fe_space) noexcept : fe_space_(&fe_space) { } template struct Jacobian : fdapde::MatrixBase> { + using TrialSpace = std::decay_t; using Derived = Derived_; template using Meta = Jacobian; using InputType = internals::fe_assembler_packet; @@ -201,7 +203,7 @@ struct vector_trial_function_impl : public fdapde::MatrixBase::type xpr_; }; template struct Divergence : fdapde::ScalarBase> { + using TrialSpace = std::decay_t; using Derived = Derived_; template using Meta = Divergence; using InputType = internals::fe_assembler_packet; using Scalar = typename Derived::Scalar; static constexpr int StaticInputSize = TrialSpace::local_dim; static constexpr int NestAsRef = 0; - static constexpr int XprBits = 0 | fe_assembler_flags::compute_shape_div; + static constexpr int XprBits = 0 | int(fe_assembler_flags::compute_shape_div); explicit constexpr Divergence(const Derived_& xpr) noexcept : xpr_(xpr) { } // fe assembly evaluation @@ -241,8 +244,8 @@ struct vector_trial_function_impl : public fdapde::MatrixBase -struct TestFunction : + requires(std::is_same_v::space_category, finite_element>) +struct TestFunction : public std::conditional_t< - FeSpace_::n_components == 1, internals::scalar_test_function_impl, - internals::vector_test_function_impl> { + FeSpace_::n_components == 1, internals::fe_scalar_test_function_impl, + internals::fe_vector_test_function_impl> { using Base = std::conditional_t< - FeSpace_::n_components == 1, internals::scalar_test_function_impl, - internals::vector_test_function_impl>; + FeSpace_::n_components == 1, internals::fe_scalar_test_function_impl, + internals::fe_vector_test_function_impl>; constexpr TestFunction() = default; constexpr TestFunction(FeSpace_& fe_space) : Base(fe_space) { } }; // partial derivative of scalar test function template -struct PartialDerivative, 1> : public TestFunction::FirstPartialDerivative { +struct PartialDerivative, 1> : + public TestFunction::FirstPartialDerivative { + using Base = TestFunction::FirstPartialDerivative; PartialDerivative() = default; - PartialDerivative(const TestFunction& f, int i) : TestFunction::FirstPartialDerivative(f, i) { } + PartialDerivative(const TestFunction& f, int i) : Base(f, i) { } }; // gradient of vectorial test function template requires(FeSpace_::n_components > 1) -constexpr auto grad(const TestFunction& xpr) { - return typename TestFunction::template Jacobian>(xpr); +constexpr auto grad(const TestFunction& xpr) { + return + typename TestFunction::template Jacobian>(xpr); } // divergence of vectorial test function template requires(FeSpace_::n_components > 1) -constexpr auto div(const TestFunction& xpr) { - return typename TestFunction::template Divergence>(xpr); +constexpr auto div(const TestFunction& xpr) { + return + typename TestFunction::template Divergence>(xpr); } template -struct TrialFunction : public std::conditional_t< - FeSpace_::n_components == 1, internals::scalar_trial_function_impl, - internals::vector_trial_function_impl> { + requires(std::is_same_v::space_category, finite_element>) +struct TrialFunction : + public std::conditional_t< + FeSpace_::n_components == 1, internals::fe_scalar_trial_function_impl, + internals::fe_vector_trial_function_impl> { using Base = std::conditional_t< - FeSpace_::n_components == 1, internals::scalar_trial_function_impl, - internals::vector_trial_function_impl>; + FeSpace_::n_components == 1, internals::fe_scalar_trial_function_impl, + internals::fe_vector_trial_function_impl>; using TrialSpace = typename Base::TrialSpace; static constexpr int local_dim = FeSpace_::local_dim; static constexpr int embed_dim = FeSpace_::embed_dim; @@ -314,22 +324,26 @@ struct TrialFunction : public std::conditional_t< // partial derivative of scalar trial function template -struct PartialDerivative, 1> : public TrialFunction::FirstPartialDerivative { +struct PartialDerivative, 1> : + public TrialFunction::FirstPartialDerivative { + using Base = TrialFunction::FirstPartialDerivative; PartialDerivative() = default; - PartialDerivative(const TrialFunction& f, int i) : - TrialFunction::FirstPartialDerivative(f, i) { } + PartialDerivative(const TrialFunction& f, int i) : Base(f, i) { } }; // gradient of vectorial trial function template -constexpr auto grad(const TrialFunction& xpr) - requires(FeSpace_::n_components > 1) { - return typename TrialFunction::template Jacobian>(xpr); + requires(FeSpace_::n_components > 1) +constexpr auto grad(const TrialFunction& xpr) { + return + typename TrialFunction::template Jacobian>(xpr); } // divergence of vectorial trial function template requires(FeSpace_::n_components > 1) -constexpr auto div(const TrialFunction& xpr) { - return typename TrialFunction::template Divergence>(xpr); +constexpr auto div(const TrialFunction& xpr) { + return + typename TrialFunction::template Divergence>( + xpr); } // representation of u(x) = \sum_{i=1}^{n_dofs} u_i \psi_i(x) with \{ \psi_i \}_i a finite element basis system @@ -435,18 +449,18 @@ class FeFunction : // getters const DVector& coeff() const { return coeff_; } - constexpr FeSpace& fe_space() { return *fe_space_; } - constexpr const FeSpace& fe_space() const { return *fe_space_; } + constexpr FeSpace& function_space() { return *fe_space_; } + constexpr const FeSpace& function_space() const { return *fe_space_; } constexpr int rows() const { return Rows; } constexpr int cols() const { return Cols; } constexpr int input_size() const { return StaticInputSize; } void set_coeff(const DVector& coeff) { coeff_ = coeff; } // linear algebra between fe functions friend constexpr FeFunction operator+(FeFunction& lhs, FeFunction& rhs) { - return FeFunction(lhs.fe_space(), lhs.coeff() + rhs.coeff()); + return FeFunction(lhs.function_space(), lhs.coeff() + rhs.coeff()); } friend constexpr FeFunction operator-(FeFunction& lhs, FeFunction& rhs) { - return FeFunction(lhs.fe_space(), lhs.coeff() - rhs.coeff()); + return FeFunction(lhs.function_space(), lhs.coeff() - rhs.coeff()); } // assignment from expansion coefficient vector FeFunction& operator=(const DVector& coeff) { diff --git a/fdaPDE/finite_elements/fe_map.h b/fdaPDE/finite_elements/fe_map.h index 1400d91..4e396e8 100644 --- a/fdaPDE/finite_elements/fe_map.h +++ b/fdaPDE/finite_elements/fe_map.h @@ -119,12 +119,12 @@ struct FeMap> : int cell_id = begin.index(); for (int i = 0, n = nodes.rows(); i < n; ++i) { // map node to reference cell and evaluate - auto cell = xpr_->fe_space().dof_handler().cell(cell_id + i / n_quad_nodes); + auto cell = xpr_->function_space().dof_handler().cell(cell_id + i / n_quad_nodes); SVector ref_node = cell.invJ() * (nodes.row(i).transpose() - cell.node(0)); DVector active_dofs = cell.dofs(); Scalar value = 0; - for (int i = 0, n = xpr_->fe_space().n_shape_functions(); i < n; ++i) { - value += xpr_->coeff()[active_dofs[i]] * xpr_->fe_space().eval(i, ref_node); + for (int i = 0, n = xpr_->function_space().n_shape_functions(); i < n; ++i) { + value += xpr_->coeff()[active_dofs[i]] * xpr_->function_space().eval(i, ref_node); } mapped(i, 0) = value; } diff --git a/fdaPDE/splines/bs_assembler_base.h b/fdaPDE/splines/bs_assembler_base.h index e156d24..5891164 100644 --- a/fdaPDE/splines/bs_assembler_base.h +++ b/fdaPDE/splines/bs_assembler_base.h @@ -16,12 +16,14 @@ #ifndef __BS_ASSEMBLER_BASE_H__ #define __BS_ASSEMBLER_BASE_H__ - + +#include "bs_integration.h" + namespace fdapde{ template struct BsMap; -enum bs_assembler_flags { +enum class bs_assembler_flags { compute_shape_values = 0x0001, compute_shape_dx = 0x0002, compute_shape_ddx = 0x0004, @@ -78,84 +80,110 @@ struct bs_assembler_base { bs_assembler_base() = default; bs_assembler_base( - const Form_& form, typename fe_traits::geo_iterator begin, typename fe_traits::geo_iterator end, - const Quadrature_&... quadrature) - requires(sizeof...(quadrature) <= 1): + const Form_& form, const geo_iterator& begin, const geo_iterator& end, const Quadrature_&... quadrature) + requires(sizeof...(quadrature) <= 1) + : form_(meta::xpr_wrap() { - return !( - std::is_invocable_v>); - })>(form)), + return !(std::is_invocable_v>); + })>(form)), dof_handler_(std::addressof(internals::test_space(form_).dof_handler())), - test_space_ (std::addressof(internals::test_space(form_))), + test_space_(std::addressof(internals::test_space(form_))), begin_(begin), end_(end) { fdapde_assert(dof_handler_->n_dofs() > 0); // copy quadrature rule + Eigen::Matrix quad_nodes__; if constexpr (sizeof...(quadrature) == 1) { auto quad_rule = std::get<0>(std::make_tuple(quadrature...)); + fdapde_assert(local_dim == quad_rule.local_dim); + quad_nodes__.resize(quad_rule.order, quad_rule.local_dim); + quad_weights_.resize(quad_rule.order, 1); for (int i = 0; i < quad_rule.order; ++i) { - quad_nodes_ [i] = quad_rule.nodes [i]; - quad_weights_[i] = quad_rule.weights[i]; - } + quad_weights_(i, 0) = quad_rule.weights[i]; + for (int j = 0; j < local_dim; ++j) { quad_nodes__(i, j) = quad_rule.nodes(i, j); } + } } else { - internals::get_bs_quadrature(test_space_->order(), quad_nodes_, quad_weights_); + internals::get_bs_quadrature(test_space_->order(), quad_nodes__, quad_weights_); + } + // build grid of quadrature nodes on reference domain + n_quadrature_nodes_ = quad_nodes__.rows(); + int n_cells = end_.index() - begin_.index(), n_src_points = quad_nodes__.rows(); // -------------------------------------- + quad_nodes_.resize(n_cells * n_src_points, local_dim); + int i = 0; + for (auto it = begin_; it != end_; ++it) { + // map src nodes on cell it + double a = it->nodes()[0], b = it->nodes()[1]; + for (int j = 0; j < quad_nodes__.rows(); ++j) { + quad_nodes_(i, 0) = (b - a) / 2 * quad_nodes__(j, 0) + (b + a) / 2; + i++; + } } } const TestSpace& test_space() const { return *test_space_; } protected: // evaluation of \psi_i(q_j), i = 1, ..., n_basis, j = 1, ..., n_quadrature_nodes - template - MdArray> eval_shape_values(BasisType__&& basis) { + template + MdArray> + eval_shape_values(BasisType__&& basis, const std::vector& active_dofs, IteratorType cell) const { + std::cout << "###########################" << std::endl; using BasisType = std::decay_t; - int n_quadrature_nodes = quad_nodes_.rows(); - int n_basis = basis.size(); - - MdArray> shape_values_(n_basis, n_quadrature_nodes); + int n_basis = active_dofs.size(); + MdArray> shape_values_(n_basis, n_quadrature_nodes_); for (int i = 0; i < n_basis; ++i) { // evaluation of \psi_i at q_j, j = 1, ..., n_quadrature_nodes - for (int j = 0; j < n_quadrature_nodes; ++j) { - shape_values_(i, j) = basis[i](quad_nodes_.row(j).transpose()); + for (int j = 0; j < n_quadrature_nodes_; ++j) { + std::cout << "cell_id: " << cell->id() << std::endl; + std::cout << cell->id() * n_quadrature_nodes_ + j << std::endl; + std::cout << quad_nodes_.rows() << std::endl; + shape_values_(i, j) = + basis[active_dofs[i]](quad_nodes_.row(cell->id() * n_quadrature_nodes_ + j).transpose()); + std::cout << i << ", " << j << std::endl; } } + std::cout << "FATTO EVAL VALUES" << std::endl; return shape_values_; } // evaluation of k-th order derivative of basis function - template + template requires(requires(BasisType__ basis, int i, int k) { basis[i].gradient(k); }) - MdArray> eval_shape_derivative(BasisType__&& basis, int order) { + MdArray> eval_shape_derivative( + BasisType__&& basis, const std::vector& active_dofs, IteratorType cell, int order) const { using BasisType = std::decay_t; using DerivativeType = decltype(std::declval()[std::declval()].gradient(std::declval())); - int n_quadrature_nodes = quad_nodes_.rows(); - int n_basis = basis.size(); + int n_basis = active_dofs.size(); // instantiate derivative functors once - std::vector der(n_basis); - for (int i = 0; i < n_basis; ++i) { ders[i] = basis[i].gradient(order); } + std::vector basis_der(n_basis); + for (int i = 0; i < n_basis; ++i) { basis_der[i] = basis[i].gradient(order); } - MdArray> shape_derivatives_(n_basis, n_quadrature_nodes); + MdArray> shape_derivatives_(n_basis, n_quadrature_nodes_); for (int i = 0; i < n_basis; ++i) { // evaluation of d^k\psi_i/dx^k at q_j, j = 1, ..., n_quadrature_nodes - for (int j = 0; j < n_quadrature_nodes; ++j) { - shape_derivatives_(i, j) = der[i](quad_nodes_.row(j).transpose()); + for (int j = 0; j < n_quadrature_nodes_; ++j) { + shape_derivatives_(i, j) = + basis_der[i](quad_nodes_.row(cell->id() * n_quadrature_nodes_ + j).transpose()); } } return shape_derivatives_; } - template MdArray> eval_shape_dx (BasisType__&& basis) { - return eval_shape_derivatives(basis, /* order = */ 1); + template + MdArray> + eval_shape_dx(BasisType__&& basis, const std::vector& active_dofs, IteratorType cell) const { + return eval_shape_derivative(basis, active_dofs, cell, /* order = */ 1); } - template MdArray> eval_shape_ddx(BasisType__&& basis) { - return eval_shape_derivatives(basis, /* order = */ 2); + template + MdArray> + eval_shape_ddx(BasisType__&& basis, const std::vector& active_dofs, IteratorType cell) const { + return eval_shape_derivative(basis, active_dofs, cell, /* order = */ 2); } void distribute_quadrature_nodes( std::unordered_map>& bs_map_buff, dof_iterator begin, dof_iterator end) { Eigen::Matrix quad_nodes; - int n_quadrature_nodes = quad_nodes_.rows(); - quad_nodes.resize(n_quadrature_nodes * (end_.begin() - begin_.index()), embed_dim); + quad_nodes.resize(n_quadrature_nodes_ * (end_.index() - begin_.index()), embed_dim); int local_cell_id = 0; for (geo_iterator it = begin_; it != end_; ++it) { - for (int q_k = 0; q_k < n_quadrature_nodes; ++q_k) { - quad_nodes.row(local_cell_id * n_quadrature_nodes + q_k) = + for (int q_k = 0; q_k < n_quadrature_nodes_; ++q_k) { + quad_nodes.row(local_cell_id * n_quadrature_nodes_ + q_k) = it->J() * quad_nodes_.row(q_k).transpose() + it->node(0); } local_cell_id++; @@ -173,10 +201,11 @@ struct bs_assembler_base { } Form form_; - Eigen::Matrix quad_nodes_, quad_weights_; + Eigen::Matrix quad_nodes_, quad_weights_; const DofHandlerType* dof_handler_; const TestSpace* test_space_; - typename fe_traits::geo_iterator begin_, end_; + geo_iterator begin_, end_; + int n_quadrature_nodes_; }; } // namespace internals diff --git a/fdaPDE/splines/bs_bilinear_form_assembler.h b/fdaPDE/splines/bs_bilinear_form_assembler.h index db94051..ac46cf3 100644 --- a/fdaPDE/splines/bs_bilinear_form_assembler.h +++ b/fdaPDE/splines/bs_bilinear_form_assembler.h @@ -35,51 +35,37 @@ class bs_bilinear_form_assembly_loop : static_assert(TrialSpace::local_dim == TestSpace::local_dim && TrialSpace::embed_dim == TestSpace::embed_dim); static constexpr bool is_galerkin = std::is_same_v; static constexpr bool is_petrov_galerkin = !is_galerkin; - using Base = fe_assembler_base; + using Base = bs_assembler_base; using Form = typename Base::Form; using DofHandlerType = typename Base::DofHandlerType; static constexpr int local_dim = Base::local_dim; static constexpr int embed_dim = Base::embed_dim; using Base::form_; + using Base::test_space_; // private data members const DofHandlerType* trial_dof_handler_; - Quadrature quadrature_ {}; constexpr const DofHandlerType* test_dof_handler() const { return Base::dof_handler_; } constexpr const DofHandlerType* trial_dof_handler() const { return is_galerkin ? Base::dof_handler_ : trial_dof_handler_; } const TrialSpace* trial_space_; - MdArray> shape_values_, shape_dx_, shape_ddx_; public: bs_bilinear_form_assembly_loop() = default; bs_bilinear_form_assembly_loop( - const Form_& form, typename Base::fe_traits::geo_iterator begin, typename Base::fe_traits::geo_iterator end, + const Form_& form, typename Base::geo_iterator begin, typename Base::geo_iterator end, const Quadrature_&... quadrature) requires(sizeof...(quadrature) <= 1) - : Base(form, begin, end, quadrature...), trial_space_(std::adressof(internals::trial_space(form_))) { + : Base(form, begin, end, quadrature...), trial_space_(std::addressof(internals::trial_space(form_))) { if constexpr (is_petrov_galerkin) { trial_dof_handler_ = std::addressof(internals::trial_space(form_).dof_handler()); } fdapde_assert(test_dof_handler()->n_dofs() != 0 && trial_dof_handler()->n_dofs() != 0); - if constexpr (sizeof...(Quadrature_) == 0) { - // default to higher-order quadrature - internals::get_bs_quadrature( - test_space_->order() > trial_space_->order() ? test_space_->order() : trial_space_->order(), - Base::quad_nodes_, Base::quad_weights_); - } - // evaluate basis system on reference interval at quadrature nodes according to form requests - if constexpr (Form::XprBits & bs_assembler_flags::compute_shape_values) { - test_shape_values_ = Base::eval_shape_values(test_space_ ->basis()); - trial_shape_values_ = Base::eval_shape_values(trial_space_->basis()); - } - if constexpr (Form::XprBits & bs_assembler_flags::compute_shape_dx) { - test_shape_dx_ = Base::eval_shape_dx(test_space_ ->basis()); - trial_shape_dx_ = Base::eval_shape_dx(trial_space_->basis()); - } - if constexpr (Form::XprBits & bs_assembler_flags::compute_shape_ddx) { - test_shape_ddx_ = Base::eval_shape_ddx(test_space_ ->basis()); - trial_shape_ddx_ = Base::eval_shape_ddx(trial_space_->basis()); - } + // if constexpr (sizeof...(Quadrature_) == 0) { + // // default to higher-order quadrature + // internals::get_bs_quadrature( + // test_space_->order() > trial_space_->order() ? test_space_->order() : trial_space_->order(), + // Base::quad_nodes_, Base::quad_weights_); // ------------------------------------------------- this should be executed only if trial space order != test space order + // } } SpMatrix assemble() const { @@ -94,51 +80,70 @@ class bs_bilinear_form_assembly_loop : void assemble(std::vector>& triplet_list) const { using iterator = typename Base::dof_iterator; iterator begin(Base::begin_.index(), test_dof_handler(), Base::begin_.marker()); - iterator begin(Base::end_.index(), test_dof_handler(), Base::end_.marker()); + iterator end (Base::end_.index() , test_dof_handler(), Base::end_.marker() ); // prepare assembly loop - Eigen::Matrix test_active_dofs, trial_active_dofs; - int n_trial_basis = trial_space_->n_shape_functions(); - int n_test_basis = test_space_->n_shape_functions(); - int n_quadrature_nodes = Base::quad_nodes_.rows(); + std::vector test_active_dofs, trial_active_dofs; + int n_quadrature_nodes = Base::quad_nodes_.rows(); + MdArray> test_shape_values_, test_shape_dx_, test_shape_ddx_; + MdArray> trial_shape_values_, trial_shape_dx_, trial_shape_ddx_; std::unordered_map> bs_map_buff; - if constexpr (Form::XprBits & bs_assembler_flags::compute_physical_quad_nodes) { + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_physical_quad_nodes)) { Base::distribute_quadrature_nodes( bs_map_buff, begin, end); // distribute quadrature nodes on physical mesh (if required) } - + // start assembly loop - internals::bs_assembler_packet bs_packet(); + internals::bs_assembler_packet bs_packet {}; int local_cell_id = 0; for (iterator it = begin; it != end; ++it) { bs_packet.cell_measure = it->measure(); // perform integration of weak form for (i, j)-th basis pair test_active_dofs = it->dofs(); if constexpr (is_petrov_galerkin) { trial_active_dofs = trial_dof_handler()->active_dofs(it->id()); } - for (int i = 0; i < n_trial_basis; ++i) { // trial function loop - for (int j = 0; j < n_test_basis; ++j) { // test function loop + + // compute values of active shape functions at quadrature nodes mapped on element *it + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_shape_values)) { + test_shape_values_ = Base::eval_shape_values(test_space_->basis(), test_active_dofs, it); + trial_shape_values_ = Base::eval_shape_values( + trial_space_->basis(), is_petrov_galerkin ? trial_active_dofs : test_active_dofs, it); + } + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_shape_dx)) { + test_shape_dx_ = Base::eval_shape_dx(test_space_->basis(), test_active_dofs, it); + trial_shape_dx_ = Base::eval_shape_dx( + trial_space_->basis(), is_petrov_galerkin ? trial_active_dofs : test_active_dofs, it); + } + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_shape_ddx)) { + test_shape_ddx_ = Base::eval_shapeddx(test_space_->basis(), test_active_dofs, it); + trial_shape_ddx_ = Base::eval_shape_ddx( + trial_space_->basis(), is_petrov_galerkin ? trial_active_dofs : test_active_dofs, it); + } + + for (int i = 0, n_trial_basis = is_petrov_galerkin ? trial_active_dofs.size() : test_active_dofs.size(); + i < n_trial_basis; ++i) { // trial function loop + for (int j = 0, n_test_basis = test_active_dofs.size(); j < n_test_basis; ++j) { // test function loop double value = 0; for (int q_k = 0; q_k < n_quadrature_nodes; ++q_k) { - if constexpr (Form::XprBits & bs_assembler_flags::compute_shape_values) { - bs_packet.trial_value = Base::trial_shape_values_(i, q_k); - bs_packet.test_value = Base::test_shape_values_ (j, q_k); + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_shape_values)) { + bs_packet.trial_value = trial_shape_values_(i, q_k); + bs_packet.test_value = test_shape_values_ (j, q_k); } - if constexpr (Form::XprBits & bs_assembler_flags::compute_shape_dx) { - bs_packet.trial_dx = Base::trial_shape_dx_(i, q_k); - bs_packet.test_dx = Base::test_shape_dx_ (j, q_k); + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_shape_dx)) { + bs_packet.trial_dx = trial_shape_dx_(i, q_k); + bs_packet.test_dx = test_shape_dx_ (j, q_k); } - if constexpr (Form::XprBits & bs_assembler_flags::compute_shape_ddx) { - bs_packet.trial_ddx = Base::trial_shape_ddx_(i, q_k); - bs_packet.test_ddx = Base::test_shape_ddx_ (j, q_k); + if constexpr (Form::XprBits & int(bs_assembler_flags::compute_shape_ddx)) { + bs_packet.trial_ddx = trial_shape_ddx_(i, q_k); + bs_packet.test_ddx = test_shape_ddx_ (j, q_k); } - if constexpr (Form::XprBits & fe_assembler_flags::compute_physical_quad_nodes) { + if constexpr (Form::XprBits & int(fe_assembler_flags::compute_physical_quad_nodes)) { bs_packet.quad_node_id = local_cell_id * n_quadrature_nodes + q_k; } - value += Base::quad_weights_[q_k] * form_(bs_packet); + value += Base::quad_weights_(q_k, 0) * form_(bs_packet); } triplet_list.emplace_back( test_active_dofs[j], is_galerkin ? test_active_dofs[i] : trial_active_dofs[i], - value * fe_packet.cell_measure); + value * bs_packet.cell_measure); } } } diff --git a/fdaPDE/splines/bs_function.h b/fdaPDE/splines/bs_function.h index e69de29..e0d6da1 100644 --- a/fdaPDE/splines/bs_function.h +++ b/fdaPDE/splines/bs_function.h @@ -0,0 +1,326 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __BS_FUNCTION_H__ +#define __BS_FUNCTION_H__ + +#include "../fields/scalar_field.h" +#include "bs_assembler_base.h" + +namespace fdapde { +namespace internals { + +template +struct bs_scalar_test_function_impl : + public fdapde::ScalarBase> { + fdapde_static_assert(BsSpace_::local_dim == 1, THIS_CLASS_IS_FOR_INTERVAL_MESHES_ONLY); + using TestSpace = std::decay_t; + using Base = fdapde::ScalarBase>; + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = TestSpace::local_dim; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = 0 | int(bs_assembler_flags::compute_shape_values); + private: + template + struct FirstDerivative_ : fdapde::ScalarBase> { + using Derived = Derived_; + template using Meta = FirstDerivative_; + using TestSpace = std::decay_t; // required from xpr_query<> + using Base = fdapde::ScalarBase>; + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = TestSpace::local_dim; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = 0 | int(bs_assembler_flags::compute_shape_dx); + + FirstDerivative_() noexcept = default; + FirstDerivative_(const Derived_& xpr) noexcept : xpr_(xpr) { } + // assembly evaluation + constexpr Scalar operator()(const InputType& bs_packet) const { return bs_packet.test_dx; } + constexpr TestSpace& function_space() { return *(xpr_.bs_space_); } + constexpr const TestSpace& function_space() const { return *(xpr_.bs_space_); } + constexpr int input_size() const { return StaticInputSize; } + constexpr const Derived& derived() const { return xpr_; } + private: + typename internals::ref_select::type xpr_; + }; + template + struct SecondDerivative_ : fdapde::ScalarBase> { + using Derived = Derived_; + template using Meta = SecondDerivative_; + using TestSpace = std::decay_t; // required from xpr_query<> + using Base = fdapde::ScalarBase>; + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = TestSpace::local_dim; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = 0 | int(bs_assembler_flags::compute_shape_ddx); + + SecondDerivative_() noexcept = default; + SecondDerivative_(const Derived_& xpr) noexcept : xpr_(xpr) { } + // assembly evaluation + constexpr Scalar operator()(const InputType& bs_packet) const { return bs_packet.test_ddx; } + constexpr TestSpace& function_space() { return *(xpr_.bs_space_); } + constexpr const TestSpace& function_space() const { return *(xpr_.bs_space_); } + constexpr int input_size() const { return StaticInputSize; } + constexpr const Derived& derived() const { return xpr_; } + private: + typename internals::ref_select::type xpr_; + }; + public: + // expose derivative types + using FirstDerivative = FirstDerivative_ >; + using SecondDerivative = SecondDerivative_>; + + constexpr bs_scalar_test_function_impl() noexcept = default; + constexpr bs_scalar_test_function_impl(BsSpace_& bs_space) noexcept : bs_space_(std::addressof(bs_space)) { } + // assembly evaluation + constexpr Scalar operator()(const InputType& bs_packet) const { return bs_packet.test_value; } + constexpr TestSpace& function_space() { return *bs_space_; } + constexpr const TestSpace& function_space() const { return *bs_space_; } + constexpr int input_size() const { return StaticInputSize; } + private: + TestSpace* bs_space_; +}; + +template +struct bs_scalar_trial_function_impl : + public fdapde::ScalarBase> { + fdapde_static_assert(BsSpace_::local_dim == 1, THIS_CLASS_IS_FOR_INTERVAL_MESHES_ONLY); + using TrialSpace = std::decay_t; + using Base = fdapde::ScalarBase>; + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = TrialSpace::local_dim; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = 0 | int(bs_assembler_flags::compute_shape_values); + private: + // definition of derivative functors + template + struct FirstDerivative_ : fdapde::ScalarBase> { + using Derived = Derived_; + template using Meta = FirstDerivative_; + using TrialSpace = std::decay_t; // required from xpr_query<> + using Base = fdapde::ScalarBase>; + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = TrialSpace::local_dim; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = 0 | int(bs_assembler_flags::compute_shape_dx); + + FirstDerivative_() noexcept = default; + FirstDerivative_(const Derived_& xpr) noexcept : xpr_(xpr) { } + // assembly evaluation + constexpr Scalar operator()(const InputType& bs_packet) const { return bs_packet.trial_dx; } + constexpr TrialSpace& function_space() { return *(xpr_.bs_space_); } + constexpr const TrialSpace& function_space() const { return *(xpr_.bs_space_); } + constexpr int input_size() const { return StaticInputSize; } + constexpr const Derived& derived() const { return xpr_; } + private: + typename internals::ref_select::type xpr_; + }; + template + struct SecondDerivative_ : fdapde::ScalarBase> { + using Derived = Derived_; + template using Meta = SecondDerivative_; + using TrialSpace = std::decay_t; // required from xpr_query<> + using Base = fdapde::ScalarBase>; + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = TrialSpace::local_dim; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = 0 | int(bs_assembler_flags::compute_shape_ddx); + + SecondDerivative_() noexcept = default; + SecondDerivative_(const Derived_& xpr) noexcept : xpr_(xpr) { } + // assembly evaluation + constexpr Scalar operator()(const InputType& bs_packet) const { return bs_packet.trial_ddx; } + constexpr TrialSpace& function_space() { return *(xpr_.bs_space_); } + constexpr const TrialSpace& function_space() const { return *(xpr_.bs_space_); } + constexpr int input_size() const { return StaticInputSize; } + constexpr const Derived& derived() const { return xpr_; } + private: + typename internals::ref_select::type xpr_; + }; + public: + // expose derivative types + using FirstDerivative = FirstDerivative_ >; + using SecondDerivative = SecondDerivative_>; + + constexpr bs_scalar_trial_function_impl() noexcept = default; + constexpr bs_scalar_trial_function_impl(BsSpace_& bs_space) noexcept : bs_space_(std::addressof(bs_space)) { } + // assembly evaluation + constexpr Scalar operator()(const InputType& bs_packet) const { return bs_packet.trial_value; } + constexpr TrialSpace& function_space() { return *bs_space_; } + constexpr const TrialSpace& function_space() const { return *bs_space_; } + constexpr int input_size() const { return StaticInputSize; } + private: + TrialSpace* bs_space_; +}; + +} // namespace internals + +template + requires(std::is_same_v::space_category, bspline>) +struct TestFunction : public internals::bs_scalar_test_function_impl { + using Base = internals::bs_scalar_test_function_impl; + constexpr TestFunction() = default; + constexpr TestFunction(BsSpace_& bs_space) : Base(bs_space) { } +}; +// partial derivatives of scalar test function +template +struct PartialDerivative, 1> : public TestFunction::FirstDerivative { + PartialDerivative() = default; + PartialDerivative(const TestFunction& f, [[maybe_unused]] int i) : + TestFunction::FirstDerivative(f) { } +}; +template +struct PartialDerivative, 2> : + public TestFunction::SecondDerivative { + PartialDerivative() = default; + PartialDerivative(const TestFunction& f, [[maybe_unused]] int i, [[maybe_unused]] int j) : + TestFunction::SecondDerivative(f) { } +}; +template constexpr auto dx (const TestFunction& xpr) { + return typename TestFunction::FirstDerivative(xpr); +} +template constexpr auto ddx(const TestFunction& xpr) { + return typename TestFunction::SecondDerivative(xpr); +} + +template + requires(std::is_same_v::space_category, bspline>) +struct TrialFunction : public internals::bs_scalar_trial_function_impl { + using Base = internals::bs_scalar_trial_function_impl; + using TrialSpace = typename Base::TrialSpace; + static constexpr int local_dim = BsSpace_::local_dim; + static constexpr int embed_dim = BsSpace_::embed_dim; + constexpr TrialFunction() = default; + constexpr TrialFunction(BsSpace_& bs_space) : Base(bs_space) { } + // norm evaluation + // double l2_squared_norm() { + // internals::fe_mass_assembly_loop assembler(Base::fe_space_->dof_handler()); + // return coeff_.dot(assembler.assemble() * coeff_); + // } + // double l2_norm() { return std::sqrt(l2_squared_norm()); } + const DVector& coeff() const { return coeff_; } + void set_coeff(const DVector& coeff) { coeff_ = coeff; } + private: + Eigen::Matrix coeff_; +}; +// partial derivative of scalar trial function +template +struct PartialDerivative, 1> : + public TrialFunction::FirstDerivative { + PartialDerivative() = default; + PartialDerivative(const TrialFunction& f, [[maybe_unused]] int i) : + TrialFunction::FirstDerivative(f, i) { } +}; +template +struct PartialDerivative, 2> : + public TrialFunction::SecondDerivative { + PartialDerivative() = default; + PartialDerivative(const TrialFunction& f, [[maybe_unused]] int i, [[maybe_unused]] int j) : + TrialFunction::SecondDerivative(f, i, j) { } +}; +template constexpr auto dx (const TrialFunction& xpr) { + return typename TrialFunction::FirstDerivative(xpr); +} +template constexpr auto ddx(const TrialFunction& xpr) { + return typename TrialFunction::SecondDerivative(xpr); +} + +// representation of u(x) = \sum_{i=1}^{n_dofs} u_i \psi_i(x) with \{ \psi_i \}_i a B-Spline basis system +template class BsFunction : public fdapde::ScalarBase> { + using Triangulation = typename BsSpace_::Triangulation; + fdapde_static_assert(BsSpace_::local_dim == 1, THIS_CLASS_IS_FOR_INTERVAL_MESHES_ONLY); + public: + using BsSpace = std::decay_t; + using Base = fdapde::ScalarBase>; + using DofHandlerType = typename BsSpace::DofHandlerType; + using InputType = Eigen::Matrix; + using Scalar = double; + static constexpr int StaticInputSize = BsSpace::local_dim; + static constexpr int Rows = 1; + static constexpr int Cols = 1; + static constexpr int NestAsRef = 1; + static constexpr int local_dim = Triangulation::local_dim; + static constexpr int embed_dim = Triangulation::embed_dim; + static constexpr int XprBits = 0; + + BsFunction() = default; + explicit BsFunction(BsSpace_& bs_space) : bs_space_(std::addressof(bs_space)) { + coeff_ = Eigen::Matrix::Zero(bs_space_->n_dofs()); + } + BsFunction(BsSpace_& bs_space, const Eigen::Matrix& coeff) : + bs_space_(std::addressof(bs_space)), coeff_(coeff) { + fdapde_assert(coeff.size() > 0 && coeff.size() == bs_space_->n_dofs()); + } + Scalar operator()(const InputType& p) const { + int e_id = bs_space_->triangulation().locate(p); + if (e_id == -1) return std::numeric_limits::quiet_NaN(); // return NaN if point lies outside domain + // map p to reference interval [-1, 1] + double a = bs_space_->triangulation().range()[0], b = bs_space_->triangulation().range()[1]; + double ref_p; + if constexpr (fdapde::is_subscriptable) { ref_p = p[0]; } + else { ref_p = p; } + ref_p = ((b - a) / 2) * ref_p + (b + a) / 2; + // get active dofs + typename DofHandlerType::CellType cell = bs_space_->dof_handler().cell(e_id); + std::vector active_dofs = cell.dofs(); + Scalar value = 0; + for (int i = 0, n = active_dofs.size(); i < n; ++i) { + value += coeff_[active_dofs[i]] * bs_space_->eval_shape_value(active_dofs[i], ref_p); + } + return value; + } + // norms of fe functions + // double l2_squared_norm() { + // internals::fe_mass_assembly_loop assembler(fe_space_->dof_handler()); + // return coeff_.dot(assembler.assemble() * coeff_); + // } + // double l2_norm() { return std::sqrt(l2_squared_norm()); } + + // getters + const Eigen::Matrix& coeff() const { return coeff_; } + constexpr BsSpace& function_space() { return *bs_space_; } + constexpr const BsSpace& function_space() const { return *bs_space_; } + constexpr int rows() const { return Rows; } + constexpr int cols() const { return Cols; } + constexpr int input_size() const { return StaticInputSize; } + void set_coeff(const Eigen::Matrix& coeff) { coeff_ = coeff; } + // linear algebra between fe functions + friend constexpr BsFunction operator+(BsFunction& lhs, BsFunction& rhs) { + return BsFunction(lhs.function_space(), lhs.coeff() + rhs.coeff()); + } + friend constexpr BsFunction operator-(BsFunction& lhs, BsFunction& rhs) { + return BsFunction(lhs.function_space(), lhs.coeff() - rhs.coeff()); + } + // assignment from expansion coefficient vector + BsFunction& operator=(const Eigen::Matrix& coeff) { + fdapde_assert(coeff.size() > 0 && coeff.size() == bs_space_->n_dofs()); + coeff_ = coeff; + return *this; + } + private: + Eigen::Matrix coeff_; + BsSpace* bs_space_; +}; + +} // namespace fdapde + +#endif // __BS_FUNCTION_H__ diff --git a/fdaPDE/splines/bs_integration.h b/fdaPDE/splines/bs_integration.h index e69de29..b4188eb 100644 --- a/fdaPDE/splines/bs_integration.h +++ b/fdaPDE/splines/bs_integration.h @@ -0,0 +1,134 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __BS_INTEGRATION_H__ +#define __BS_INTEGRATION_H__ + +#include "../linear_algebra/constexpr_matrix.h" + +namespace fdapde { +namespace internals { + +struct bs_quadrature_gauss_base { }; +template +concept is_bs_quadrature_gauss = std::is_base_of_v; +template [[maybe_unused]] static constexpr bool is_bs_quadrature_gauss_v = is_bs_quadrature_gauss; + +template + requires( + LhsQuadrature::local_dim == RhsQuadrature::local_dim && is_bs_quadrature_gauss_v && + is_bs_quadrature_gauss_v) +struct higher_degree_bs_quadrature : + std::type_identity< + std::conditional_t<(LhsQuadrature::order > RhsQuadrature::order), LhsQuadrature, RhsQuadrature>> { }; +template +using higher_degree_bs_quadrature_t = higher_degree_bs_quadrature::type; + +// quadrature points and weights at: https://people.sc.fsu.edu/~jburkardt/datasets/datasets.html +template struct bs_quadrature_gauss_legendre; + +// degree: highest polynomial degree correctly integrated by the formula +// order : number of quadrature nodes + +// 1D 1 point formula +template <> struct bs_quadrature_gauss_legendre<1, 1> : public bs_quadrature_gauss_base { + static constexpr int local_dim = 1; + static constexpr int order = 1; + static constexpr int degree = 2 * order - 1; // 1 + + static constexpr cexpr::Vector nodes { + std::array {0.000000000000000} + }; + static constexpr cexpr::Vector weights { + std::array {2.000000000000000} + }; +}; + +// 1D 2 point formula +template <> struct bs_quadrature_gauss_legendre<1, 2> : public bs_quadrature_gauss_base { + static constexpr int local_dim = 1; + static constexpr int order = 2; + static constexpr int degree = 2 * order - 1; // 3 + + static constexpr cexpr::Vector nodes { + std::array {-0.5773502691896257, 0.5773502691896257} + }; + static constexpr cexpr::Vector weights { + std::array { 0.9999999999999998, 0.9999999999999998} + }; +}; + +// 1D 3 point formula +template <> struct bs_quadrature_gauss_legendre<1, 3> : public bs_quadrature_gauss_base { + static constexpr int local_dim = 1; + static constexpr int order = 3; + static constexpr int degree = 2 * order - 1; // 5 + + static constexpr cexpr::Vector nodes { + std::array {-0.7745966692414834, 0.0000000000000000, 0.7745966692414834} + }; + static constexpr cexpr::Vector weights { + std::array { 0.5555555555555555, 0.8888888888888888, 0.5555555555555555} + }; +}; + +// 1D 4 point formula +template <> struct bs_quadrature_gauss_legendre<1, 4> : public bs_quadrature_gauss_base { + static constexpr int local_dim = 1; + static constexpr int order = 4; + static constexpr int degree = 2 * order - 1; // 7 + + static constexpr cexpr::Vector nodes { + std::array {-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526} + }; + static constexpr cexpr::Vector weights { + std::array { 0.3478548451374538, 0.6521451548625461, 0.6521451548625461, 0.3478548451374538} + }; +}; + +// copy weights and nodes of quadrature rule based on run-time polynomial order +template + requires(requires(T t, int i, int j) { + { t(i, j) } -> std::same_as; + { t.resize(i, j) } -> std::same_as; + }) +void get_bs_quadrature(int order, T& quad_nodes, T& quad_weights) { + fdapde_assert(order >= 0 && order <= 7); + auto copy_ = [](QuadRule q, T& quad_nodes_, T& quad_weights_) { + quad_nodes_ .resize(q.order, q.local_dim); + quad_weights_.resize(q.order, q.local_dim); + for (int i = 0; i < q.order; ++i) { + quad_nodes_ (i, 0) = q.nodes [i]; + quad_weights_(i, 0) = q.weights[i]; + } + }; + if (order == 0 || order == 1) { copy_(bs_quadrature_gauss_legendre<1, 1> {}, quad_nodes, quad_weights); } + if (order > 1 && order <= 3) { copy_(bs_quadrature_gauss_legendre<1, 3> {}, quad_nodes, quad_weights); } + if (order > 3 && order <= 5) { copy_(bs_quadrature_gauss_legendre<1, 3> {}, quad_nodes, quad_weights); } + if (order > 5 && order <= 7) { copy_(bs_quadrature_gauss_legendre<1, 4> {}, quad_nodes, quad_weights); } +} + +} // namespace internals + +// 1D formulas +[[maybe_unused]] static struct QGL1D1P_ : internals::bs_quadrature_gauss_legendre<1, 1> { } QGL1D1P; +[[maybe_unused]] static struct QGL1D2P_ : internals::bs_quadrature_gauss_legendre<1, 2> { } QGL1D2P; +[[maybe_unused]] static struct QGL1D3P_ : internals::bs_quadrature_gauss_legendre<1, 3> { } QGL1D3P; +[[maybe_unused]] static struct QGL1D4P_ : internals::bs_quadrature_gauss_legendre<1, 4> { } QGL1D4P; + +} // namespace fdapde + +#endif // __BS_INTEGRATION_H__ diff --git a/fdaPDE/splines/bs_map.h b/fdaPDE/splines/bs_map.h new file mode 100644 index 0000000..baeea0c --- /dev/null +++ b/fdaPDE/splines/bs_map.h @@ -0,0 +1,95 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __BS_MAP_H__ +#define __BS_MAP_H__ + +#include "../fields/scalar_field.h" +#include "../fields/meta.h" +#include "bs_function.h" + +namespace fdapde { + +// anytime you compose a trial or test function with a functor which is not callable at fe_assembler_packet, we wrap it +// into a BsMap, a fe_assembler_packet callable type encoding the functor evaluated at a fixed set of (quadrature) nodes +template +struct BsMap : + public std::conditional_t< + meta::is_scalar_field_v, ScalarBase>, + MatrixBase>> { + private: + using OutputType = decltype(std::declval().operator()(std::declval())); + static constexpr bool is_scalar = meta::is_scalar_field_v; + using Derived = std::decay_t; + public: + using InputType = internals::bs_assembler_packet; + using Scalar = double; + static constexpr int StaticInputSize = Derived::StaticInputSize; + using Base = std::conditional_t< + is_scalar, ScalarBase>, MatrixBase>>; + static constexpr int NestAsRef = 0; + static constexpr int XprBits = Derived::XprBits | int(bs_assembler_flags::compute_physical_quad_nodes); + static constexpr int ReadOnly = 1; + static constexpr int Rows = []() { if constexpr(is_scalar) return 1; else return Derived::Rows; }(); + static constexpr int Cols = []() { if constexpr(is_scalar) return 1; else return Derived::Cols; }(); + + constexpr BsMap() = default; + constexpr BsMap(const Derived_& xpr) : xpr_(&xpr) { } + template + void init( + std::unordered_map>& buff, const DMatrix& nodes, + [[maybe_unused]] CellIterator begin, [[maybe_unused]] CellIterator end) const { + const void* ptr = reinterpret_cast(xpr_); + if (buff.find(ptr) == buff.end()) { + DMatrix mapped(nodes.rows(), Rows * Cols); + if constexpr (is_scalar) { + for (int i = 0, n = nodes.rows(); i < n; ++i) { mapped(i, 0) = xpr_->operator()(nodes.row(i)); } + } else { + for (int i = 0, n = nodes.rows(); i < n; ++i) { + auto tmp = xpr_->operator()(nodes.row(i)); + for (int j = 0; j < tmp.size(); ++j) { mapped(i, j) = tmp[j]; } + } + } + buff[ptr] = mapped; + map_ = &buff[ptr]; + } else { + map_ = &buff[ptr]; + } + } + // fe assembler evaluation + + // here we need to better state what is the scalar field interface and what the matrix field one + + constexpr OutputType operator()(const InputType& fe_packet) const { + if constexpr (is_scalar) { + return map_->operator()(fe_packet.quad_node_id, 0); + } else { + return map_->row(fe_packet.quad_node_id); + } + } + constexpr auto eval(int i, const InputType& fe_packet) const { return map_->operator()(fe_packet.quad_node_id, i); } + constexpr const Derived& derived() const { return xpr_; } + constexpr int input_size() const { return StaticInputSize; } + constexpr int rows() const { return Rows; } + constexpr int cols() const { return Cols; } + private: + const Derived* xpr_; + mutable const DMatrix* map_; +}; + +} // namespace fdapde + +#endif // __BS_MAP_H__ diff --git a/fdaPDE/splines/bs_space.h b/fdaPDE/splines/bs_space.h index dc0e940..786538e 100644 --- a/fdaPDE/splines/bs_space.h +++ b/fdaPDE/splines/bs_space.h @@ -18,14 +18,24 @@ #define __BS_SPACE_H__ #include "../utils/symbols.h" +#include "bs_assembler_base.h" #include "dof_handler.h" namespace fdapde { +// forward declarations template class BsFunction; +namespace internals { +template +class bs_bilinear_form_assembly_loop; +template +class bs_linear_form_assembly_loop; + +} // namespace internals + template class BsSpace { - fdapde_assert( + fdapde_static_assert( Triangulation_::local_dim == 1 && Triangulation_::embed_dim == 1, THIS_CLASS_IS_FOR_INTERVAL_MESHES_ONLY); template struct subscript_t_impl { using type = std::decay_t().operator[](std::declval()))>; @@ -35,21 +45,30 @@ template class BsSpace { using Triangulation = std::decay_t; static constexpr int local_dim = Triangulation::local_dim; static constexpr int embed_dim = Triangulation::embed_dim; - using FeType = std::decay_t; using BasisType = SplineBasis; using ShapeFunctionType = subscript_t; using DofHandlerType = DofHandler; + using space_category = fdapde::bspline; + template + using bilinear_form_assembly_loop = + internals::bs_bilinear_form_assembly_loop; + template + using linear_form_assembler_loop = + internals::bs_linear_form_assembly_loop ; BsSpace() = default; BsSpace(const Triangulation_& interval, int order) : - interval_(std::addressof(interval)), dof_handler_(interval), order_(order) { + triangulation_(std::addressof(interval)), dof_handler_(interval), order_(order) { dof_handler_.enumerate(BasisType(interval, order)); // build reference [-1, 1] interval with nodes mapped from physical interval [a, b] - Eigen::Matrix ref_nodes(triangulation.n_nodes()); - for (int i = 0; i < nodes.rows(); ++i) { ref_nodes[i] = map_to_reference(triangulation.nodes(i, 0)); } - ref_interval_ = Triangulation(ref_nodes); + Eigen::Matrix ref_nodes(triangulation_->n_nodes()); + for (int i = 0; i < triangulation_->n_nodes(); ++i) { + ref_nodes[i] = map_to_reference(triangulation_->nodes()(i, 0)); + std::cout << triangulation_->nodes()(i, 0) << " : " << ref_nodes[i] << " -- "; + } + std::cout << std::endl; // generate basis on reference [-1, 1] interval - basis_ = BasisType(ref_interval, order); + basis_ = BasisType(Triangulation(ref_nodes), order); } // observers const Triangulation& triangulation() const { return *triangulation_; } @@ -60,7 +79,6 @@ template class BsSpace { int n_dofs() const { return dof_handler_.n_dofs(); } const BasisType& basis() const { return basis_; } int order() const { return order_; } - // evaluation template requires(std::is_invocable_v) @@ -90,13 +108,13 @@ template class BsSpace { // evaluation on physical domain // evaluate value of the i-th shape function defined on physical domain [a, b] template auto eval_cell_value(int i, const InputType& p) const { - if constexpr (fdapde::is_subscriptable) { - if (p[0] < triangulation.range()[0] || p[0] > triangulation.range()[1]) return 0.0; - return eval_shape_value(i, map_to_reference(p)); - } else { - if (p < triangulation.range()[0] || p > triangulation.range()[1]) return 0.0; - return eval_shape_value(i, map_to_reference(p)); + double p_; + if constexpr (fdapde::is_subscriptable) { p_ = p[0]; } + else { p_ = p; } + if (p_ < triangulation_->range()[0] || p_ > triangulation_->range()[1]) { + return std::numeric_limits::quiet_NaN(); // return NaN if point lies outside domain } + return eval_shape_value(i, map_to_reference(p_)); } // need to return something which represent a basis function on the whole physical domain @@ -108,8 +126,8 @@ template class BsSpace { private: // linear mapping from p \in [a, b] to \hat p \in [-1, +1] double map_to_reference(double p) { - double a = triangulation.range()[0], b = triangulation.range()[1]; - return ((b - a) / 2) * p + (b + a) / 2; + double a = triangulation_->range()[0], b = triangulation_->range()[1]; + return (p - (b - a) / 2) * 2 / (b + a); }; const Triangulation* triangulation_; diff --git a/fdaPDE/splines/dof_handler.h b/fdaPDE/splines/dof_handler.h index 72db6db..5e1e385 100644 --- a/fdaPDE/splines/dof_handler.h +++ b/fdaPDE/splines/dof_handler.h @@ -64,10 +64,10 @@ class DofHandler<1, 1, fdapde::bspline> { }; // constructor DofHandler() = default; - DofHandler(const TriangulationType& triangulation) : triangulation_(std::addressof(triangulation)) { } + DofHandler(const TriangulationType& triangulation) : triangulation_(std::addressof(triangulation)), order_() { } // getters CellType cell(int id) const { return CellType(id, this); } - Eigen::Map> dofs() const { + Eigen::Map> dofs() const { // dofs active on cell return Eigen::Map>( dofs_.data(), dofs_.size() / 2, 2); } @@ -90,16 +90,14 @@ class DofHandler<1, 1, fdapde::bspline> { } return result; } + // In any given knot span [u_i, u_{i+1}) at most p+1 basis functions are non zero, namely N_{i-p,p}, ..., N_{i,p} + // (property P2.2, pag 55, Piegl, L., & Tiller, W. (2012). The NURBS book. Springer Science & Business Media.) std::vector active_dofs(int i) const { - int cell_id = 2 * i; - int n_dofs = dofs_[cell_id + 1] - dofs_[cell_id] + 1; - std::vector tmp(n_dofs); - for (int j = 0; j <= n_dofs; ++j) { tmp[j] = dofs_[cell_id] + j; } - return tmp; - } - template void active_dofs(int cell_id, ContainerT& dst) const { - for (int i = dofs_[cell_id], n = dofs_[cell_id + 1]; i <= n; ++i) { dst.push_back(i); } + std::vector dofs; + for (int j = 0; j < order_ + 1; ++j) { dofs.push_back(i + j); } + return dofs; } + template void active_dofs(int i, ContainerT& dst) const { dst = active_dofs(i); } operator bool() const { return n_dofs_ != 0; } // iterates over geometric cells coupled with dofs informations (possibly filtered by marker) @@ -193,16 +191,14 @@ class DofHandler<1, 1, fdapde::bspline> { template void enumerate(BsType&& bs) { n_dofs_ = bs.size(); + order_ = bs.order(); dofs_coords_.resize(n_dofs_); for (int i = 0; i < n_dofs_; ++i) { dofs_coords_[i] = bs[i].knot(); } - int i = 0, n_cells = triangulation()->n_cells(); + int n_cells = triangulation()->n_cells(); for (int j = 0; j < n_cells; ++j) { - dofs_.push_back(i); - while (i < n_dofs_ && dofs_coords_[i] == dofs_coords_[i + 1]) { i++; } - i++; - while (i < n_dofs_ && dofs_coords_[i] == dofs_coords_[i + 1]) { i++; } - dofs_.push_back(i); - } + dofs_.push_back(j); + dofs_.push_back(j + order_); + } // Regardless of the number of physical dofs at the interval boundary, only the basis functions associated with // the first and last dofs are non-zero at the boundary nodes. Hence, we treat only these dofs as boundary dofs boundary_dofs_.resize(n_dofs_); @@ -219,6 +215,7 @@ class DofHandler<1, 1, fdapde::bspline> { int n_dofs_; std::vector dofs_markers_; const TriangulationType* triangulation_; + int order_; }; } // namespace fdapde diff --git a/fdaPDE/splines/spline_basis.h b/fdaPDE/splines/spline_basis.h index 06a9fdb..dcec852 100644 --- a/fdaPDE/splines/spline_basis.h +++ b/fdaPDE/splines/spline_basis.h @@ -77,6 +77,7 @@ class SplineBasis { constexpr int size() const { return basis_.size(); } constexpr const std::vector& knots_vector() const { return knots_; } int n_knots() const { return knots_.size(); } + int order() const { return order_; } }; } // namespace fdapde diff --git a/fdaPDE/utils/traits.h b/fdaPDE/utils/traits.h index 67e1b41..0dac0f5 100644 --- a/fdaPDE/utils/traits.h +++ b/fdaPDE/utils/traits.h @@ -28,7 +28,6 @@ namespace fdapde { template struct subscript_result_of { using type = decltype(std::declval().operator[](std::declval()...)); }; - template concept is_subscriptable = requires(Fn fn, Arg arg) { { fn.operator[](arg) }; diff --git a/fdapde_spack.yaml b/spack.yaml similarity index 100% rename from fdapde_spack.yaml rename to spack.yaml