Updated submodules.

R/function_space.R

@@ -1,16 +1,3 @@
-# .CppFunctionSpace <- setRefClass(
-#   Class = "CppFunctionSpace",
-#   fields = c(
-#     cpp_handler = "ANY"  ## pointer to cpp backend
-#   ),
-#   methods = c(
-#     ## evaluates the basis system on the given set of locations
-#     eval = function(type = c("pointwise", "areal"), locations) {
-#       type <- match.arg(type)
-#       return(cpp_handler$eval(match(type, c("pointwise", "areal")) - 1, locations))
-#     }
-#   )
-# )
 
 .BasisFunctionCtr <- setRefClass(
   Class = "BasisFunction",
@@ -78,7 +65,7 @@ setMethod("FunctionSpace",
           signature = c(mesh = "ANY", fe_order = "numeric"),
           function(mesh, fe_order) {
             .FunctionSpaceCtr(
-              basis_function = BasisFunction(mesh,fe_order),
+              basis_function = BasisFunction(mesh,as.integer(fe_order)),
               mesh = mesh,
               fe_order = as.integer(fe_order)
             )

src/include/r_functional_space.h

@@ -31,6 +31,10 @@ using fdapde::core::LagrangianBasis;
 using fdapde::core::pointwise_evaluation;
 using fdapde::core::areal_evaluation;
 using fdapde::core::eval;
+using fdapde::core::ct_binomial_coefficient;
+using fdapde::core::FEM;
+using fdapde::core::fem_order;
+using fdapde::core::combinations;
 
 // type-erasure wrapper for a function space object
 struct I_FunctionalSpace {
@@ -62,22 +66,39 @@ template <int M, int N, int R> class R_FunctionalSpace {
    private:
     using DomainType = Mesh<M, N>;
     using FunctionalSpaceType = FunctionalSpace;
-    using QuadratureRule = Integrator<DomainType::local_dimension, R>; // exact for FEM spaces (TODO: generalize)
+    using QuadratureRule = Integrator<FEM, DomainType::local_dimension, R>; // exact for FEM spaces (TODO: generalize)
 
     // internal data
     DomainType domain_ {};
     FunctionalSpaceType fun_space_ {}; // wrapped functional space
     QuadratureRule integrator_ {};     // quadrature rule (exact for the provided fem order)
     int space_type_;
+    std::size_t n_dofs_ = 0;        // degrees of freedom, i.e. the maximum ID in the dof_table_
+    DMatrix<int> dofs_;             // for each element, the degrees of freedom associated to it
+    DMatrix<int> boundary_dofs_;    // unknowns on the boundary of the domain, for boundary conditions prescription
+
+    void enumerate_dofs(const DomainType& mesh);
    public:
+
+    enum {
+        fem_order = R,
+        n_dof_per_element = ct_nnodes(DomainType::local_dimension, fem_order),
+        n_dof_per_edge = fem_order - 1,
+        n_dof_internal =
+          n_dof_per_element - (DomainType::local_dimension + 1) - DomainType::n_facets_per_element * (fem_order - 1)   // > 0 \iff R > 2
+    };
+
     // constructor
     R_FunctionalSpace(Rcpp::Environment mesh, int space_type) : space_type_(space_type) {
         // set domain
         SEXP meshptr = mesh[".pointer"];
         R_Mesh<M, N>* ptr = reinterpret_cast<R_Mesh<M, N>*>(R_ExternalPtrAddr(meshptr));
 	domain_ = ptr->domain();
+
+	enumerate_dofs(domain_);
+
 	// define functional space
-	std::size_t size = (R == 1 ? domain_.n_nodes() : domain_.n_nodes() + domain_.n_edges());
+	std::size_t size = n_dofs_;
 	switch(space_type_) {
 	case space_type::fem_lagrange: {
 	  fun_space_ = LagrangianBasis<DomainType, R>(domain_, size);
@@ -99,8 +120,58 @@ template <int M, int N, int R> class R_FunctionalSpace {
     double integrate(const DVector<double>& c) const {
       return integrator_.integrate(domain_, fun_space_(c, integrator_.quadrature_nodes(domain_)));
     }
+
     // destructor
     ~R_FunctionalSpace() = default;
 };
 
+template<int M, int N, int R>
+void R_FunctionalSpace<M,N,R>::enumerate_dofs(const Mesh<M, N>& mesh) {
+    if (n_dofs_ != 0) return;   // return early if dofs already computed
+    if constexpr (fem_order == 1) {
+      n_dofs_ = mesh.n_nodes();
+      dofs_ = mesh.elements();
+      boundary_dofs_ = mesh.boundary();
+    } else {
+      dofs_.resize(mesh.n_elements(), n_dof_per_element);
+      dofs_.leftCols(DomainType::n_vertices) = mesh.elements();   // copy dofs associated to geometric vertices
+
+      int next = mesh.n_nodes();   // next valid ID to assign
+      auto edge_pattern = combinations<DomainType::n_vertices_per_edge, DomainType::n_vertices>();
+      std::set<int> boundary_set;
+
+      // cycle over mesh edges
+      for (auto edge = mesh.facet_begin(); edge != mesh.facet_end(); ++edge) {
+            for (std::size_t i = 0; i < DomainType::n_elements_per_facet; ++i) {
+                int element_id = (*edge).adjacent_elements()[i];
+                if (element_id >= 0) {
+                    // search for dof insertion point
+                    std::size_t j = 0;
+                    for (; j < edge_pattern.rows(); ++j) {
+                        std::array<int, DomainType::n_vertices_per_edge> e {};
+                        for (std::size_t k = 0; k < DomainType::n_vertices_per_edge; ++k) {
+                            e[k] = mesh.elements()(element_id, edge_pattern(j, k));
+                        }
+                        std::sort(e.begin(), e.end());   // normalize edge ordering
+                        if ((*edge).node_ids() == e) break;
+                    }
+                    dofs_(element_id, DomainType::n_vertices + j) = next;
+                    if ((*edge).on_boundary()) boundary_set.insert(next);
+
+                    // insert any internal dofs, if any (for cubic or higher order) + insert n_dof_per_edge dofs (for
+                    // cubic or higher)
+                }
+            }
+            next++;
+      }
+
+      n_dofs_ = next;   // store number of unknowns
+      // update boundary
+      boundary_dofs_ = DMatrix<int>::Zero(n_dofs_, 1);
+      boundary_dofs_.topRows(mesh.boundary().rows()) = mesh.boundary();
+      for (auto it = boundary_set.begin(); it != boundary_set.end(); ++it) { boundary_dofs_(*it, 0) = 1; }
+    }
+    return;
+}
+
 #endif   // __R_FUNCTIONAL_SPACE_H__

src/include/r_pde.h

@@ -56,7 +56,7 @@ struct PDEWrapper {
 template <int M, int N, int R> class R_PDE : public PDEWrapper {
    private:
     using DomainType = Mesh<M, N>;
-    using QuadratureRule = Integrator<DomainType::local_dimension, R>;
+    using QuadratureRule = Integrator<FEM, DomainType::local_dimension, R>;
     template <typename L> using PDEType = PDE<DomainType, L, DMatrix<double>, FEM, fem_order<R>>;
 
     // internal data

tests/pde.2.R

@@ -8,7 +8,7 @@ class(mesh)
 plot(mesh)
 
 # create Functional Space
-fe_order = 1L
+fe_order = 2
 Vh <- FunctionSpace(mesh, fe_order)
 
 exact_solution <- function(points){
@@ -64,3 +64,14 @@ Psi <- basis_function$eval(pde$get_dofs_coordinates())
 dim(Psi)
 
 A  <- basis_function$eval(as.matrix(points))
+dim(A)
+
+R0 <- pde$get_mass()
+R1 <- pde$get_stiff()
+
+lambda <- 1.
+row_1 <- cbind( 1/n * t(Psi)%*%Psi, lambda*t(R1))
+row_2 <- cbind(lambda*R1, R0)
+
+M <- rbind(row_1,
+           row_2)