From 72c540f9bf00b4228c95f46dd7a7e1058c926586 Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Fri, 22 Dec 2023 13:20:58 +0100 Subject: [PATCH] Updated submodules. --- R/function_space.R | 15 +------ src/core | 2 +- src/include/r_functional_space.h | 75 +++++++++++++++++++++++++++++++- src/include/r_pde.h | 2 +- tests/pde.2.R | 13 +++++- 5 files changed, 88 insertions(+), 19 deletions(-) diff --git a/R/function_space.R b/R/function_space.R index bb89e50..d5b228c 100644 --- a/R/function_space.R +++ b/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) ) diff --git a/src/core b/src/core index 2c61d4f..8e7e8a6 160000 --- a/src/core +++ b/src/core @@ -1 +1 @@ -Subproject commit 2c61d4f87409506119e33a981f929fbb16c9719b +Subproject commit 8e7e8a67aee741d21cf290c1b22ced11f2dc7411 diff --git a/src/include/r_functional_space.h b/src/include/r_functional_space.h index 2eb8e24..090b2c5 100644 --- a/src/include/r_functional_space.h +++ b/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 class R_FunctionalSpace { private: using DomainType = Mesh; using FunctionalSpaceType = FunctionalSpace; - using QuadratureRule = Integrator; // exact for FEM spaces (TODO: generalize) + using QuadratureRule = Integrator; // 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 dofs_; // for each element, the degrees of freedom associated to it + DMatrix 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* ptr = reinterpret_cast*>(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(domain_, size); @@ -99,8 +120,58 @@ template class R_FunctionalSpace { double integrate(const DVector& c) const { return integrator_.integrate(domain_, fun_space_(c, integrator_.quadrature_nodes(domain_))); } + // destructor ~R_FunctionalSpace() = default; }; +template +void R_FunctionalSpace::enumerate_dofs(const Mesh& 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(); + std::set 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 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::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__ diff --git a/src/include/r_pde.h b/src/include/r_pde.h index 9629c83..554f7f9 100644 --- a/src/include/r_pde.h +++ b/src/include/r_pde.h @@ -56,7 +56,7 @@ struct PDEWrapper { template class R_PDE : public PDEWrapper { private: using DomainType = Mesh; - using QuadratureRule = Integrator; + using QuadratureRule = Integrator; template using PDEType = PDE, FEM, fem_order>; // internal data diff --git a/tests/pde.2.R b/tests/pde.2.R index db7ad4d..b452b21 100644 --- a/tests/pde.2.R +++ b/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) \ No newline at end of file