femR relies on R6 classes.

fdaPDE · Jan 3, 2024 · 3d329b6 · 3d329b6
1 parent 72c540f
commit 3d329b6
Show file tree

Hide file tree

Showing 18 changed files with 408 additions and 483 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -15,7 +15,7 @@ Authors@R: c( person("Aldo", "Clemente", role = c("aut", "cre"),
 	                                        email = "[email protected]"))
 Maintainer: Aldo Clemente <[email protected]>
 Depends: 
-    R (>= 3.5.0),
+    R (>= 3.5.0), R6,
     plotly, Matrix, methods,
     RTriangle, sf
 Description: Solving partial differential equations relying on the Finite Element Method.

diff --git a/NAMESPACE b/NAMESPACE
@@ -1,5 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method("*",DiffOpObject)
+S3method("*",Function)
+S3method("*",FunctionGrad)
+S3method("+",DiffOpObject)
+S3method("-",DiffOpObject)
 S3method("st_crs<-",Domain)
 S3method(plot,Function)
 S3method(plot,Mesh)
@@ -8,6 +13,8 @@ S3method(st_as_sfc,Domain)
 S3method(st_as_sfc,Mesh)
 S3method(st_crs,Domain)
 export("%X%")
+export(.FunctionSpaceCtr)
+export(.MeshCtr)
 export(Domain)
 export(Function)
 export(FunctionSpace)
@@ -19,9 +26,18 @@ export(dot)
 export(grad)
 export(laplace)
 export(read_mesh)
-exportMethods("*")
-exportMethods("+")
-exportMethods("-")
+exportClasses(BasisFunction)
+exportClasses(DiffOpObject)
+exportClasses(DiffusionOperator)
+exportClasses(Domain)
+exportClasses(Function)
+exportClasses(FunctionGrad)
+exportClasses(FunctionSpace)
+exportClasses(Mesh)
+exportClasses(Pde)
+exportClasses(ReactionOperator)
+exportClasses(TimeDerivative)
+exportClasses(TransportOperator)
 exportMethods(contour)
 exportMethods(dt)
 import(Matrix)

diff --git a/R/domain.R b/R/domain.R
@@ -1,12 +1,29 @@
-.DomainCtr <- setRefClass("Domain", 
-                          fields=c(
-                            geometry = "ANY",
-                            coords = "data.frame",         # (x,y,boundary)
-                            time_interval = "numeric",
-                            crs ="ANY"
-                           )
+.DomainCtr <- R6Class("Domain", 
+                      public = list(
+                        geometry = "ANY",
+                        coords = "data.frame",         # (x,y,boundary)
+                        time_interval = vector(mode="numeric", length = 0),
+                        crs = "ANY",
+                        initialize = function(geometry, coords, crs){
+                                      self$geometry <- geometry 
+                                      self$coords <- coords 
+                                      self$crs <- crs
+                                    }
+                      )
 )
 
+#' @name Domain
+#'
+#' @exportClass Domain
+setOldClass(c("Domain", "R6"))
+
+#' Mesh 
+#'
+#' @name Mesh
+#'
+#' @exportClass Mesh
+setOldClass(c("Mesh","Domain"))
+
 #' S4 class representing a spatial (spatio-temporal) domain
 #' 
 #' @param x could be a \code{pslg} returned by \code{\link[RTriangle]{triangulate}}, 
@@ -36,10 +53,11 @@ setMethod("Domain", signature = "list", function(x){
 
     coords <- data.frame(x=x$nodes[,1], y=x$nodes[,2], boundary=x$nodes_boundary)
   crs = NA
-  .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
-             coords=coords, crs = crs)
+  .DomainCtr$new(geometry = geometry, coords=coords, crs = crs)
 })
 
+setOldClass("pslg")
+
 #' @rdname Domain
 setMethod("Domain", signature = "pslg", function(x){
   nodes <- x$P
@@ -62,8 +80,7 @@ setMethod("Domain", signature = "pslg", function(x){
 
   coords <- data.frame(x=nodes[,1], y=nodes[,2], boundary=nodes_boundary)
 
-  .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0),
-             coords=coords, crs = NA)
+  .DomainCtr$new(geometry = geometry, coords=coords, crs = NA)
 })
 
 #' @rdname Domain
@@ -206,8 +223,7 @@ setMethod("Domain", signature="sfc", function(x){
 
   crs <- NA
   if(!is.na(st_crs(x))) crs <- st_crs(x)$input
-      .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
-                coords=coords, crs = crs)
+      .DomainCtr$new(geometry = geometry, coords=coords, crs = crs)
 })
 
 # setMethod("Domain", signature = "sfc", function(x){
@@ -222,7 +238,7 @@ setMethod("Domain", signature="sfc", function(x){
 #                    nodes_group = nodes_group, edges_group = edges_group)
 #   crs <- NA_crs_
 #   if(!is.na(st_crs(x))) crs <- st_crs(x)
-#   .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
+#   .DomainCtr$new(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
 #              coords=geometry$nodes, crs = crs)
 # })
 

diff --git a/R/function_space.R b/R/function_space.R
@@ -1,57 +1,76 @@
-
-.BasisFunctionCtr <- setRefClass(
-  Class = "BasisFunction",
-  fields = c(
-    cpp_handler = "ANY"  ## pointer to cpp backend
-  ),
-  methods = c(
-    ## evaluates the basis system on the given set of locations
-    eval = function(locations) {
-      return(cpp_handler$eval(0L, locations))
-    },
-    ## integrates a function expressed as basis expansion with respect to this basis system over the
-    ## whole domain of definition
-    integrate = function(func) {
-      if (!(is.function(func) || is.vector(func))) stop("invalid argument, should be either a function or a vector")
-      ## required dof coordinates...
-      c <- func
-      if (is.function(func)) { ## recover fem expansion
-        c <- func(mesh$cpp_handler$nodes())
-      }
-      return(cpp_handler$integrate(c))
-    }
+.BasisFunctionCtr <- R6Class("BasisFunction",
+                              public = list(
+                                cpp_handler = "ANY",  ## pointer to cpp backend
+                                mesh = "Mesh",
+                                initialize = function(cpp_handler, mesh){
+                                  self$cpp_handler = cpp_handler
+                                  self$mesh = mesh
+                                },
+                                ## evaluates the basis system on the given set of locations
+                                eval = function(locations) {
+                                  return(self$cpp_handler$eval(0L, locations))
+                                },
+                                ## integrates a function expressed as basis expansion with respect to this basis system over the
+                                ## whole domain of definition
+                                integrate = function(func) {
+                                  if (!(is.function(func) || is.vector(func))) stop("invalid argument, should be either a function or a vector")
+                                  ## required dof coordinates...
+                                  c <- func
+                                  if (is.function(func)) { ## recover fem expansion
+                                    c <- func(self$mesh$cpp_handler$nodes())
+                                  }
+                                  return(cpp_handler$integrate(c))
+                                }
   )
 )
 
+#' Finite Element Basis 
+#'
+#' @name Basis
+#'
+#' @exportClass BasisFunction
+setOldClass(c("BasisFunction", "R6"))
+
 setGeneric("BasisFunction", function(mesh, fe_order) standardGeneric("BasisFunction"))
 
 setMethod("BasisFunction", signature = signature("Mesh","integer"),
           function(mesh, fe_order){
-            basis_function <- .BasisFunctionCtr(cpp_handler = 
+            basis_function <- .BasisFunctionCtr$new(cpp_handler = 
                                                   new(eval(parse(text = paste("cpp_lagrange_basis_2d_fe", 
-                                                                              as.character(fe_order), sep = ""))), mesh$cpp_handler, 0))
+                                                                              as.character(fe_order), sep = ""))), mesh$cpp_handler, 0),
+                                                  mesh=mesh)
 })
 
-.FunctionSpaceCtr <- setRefClass(
-  Class = "FunctionSpaceObject",
-  fields = c(
-    mesh = "ANY",
-    basis_function = "BasisFunction",     ## cpp backend
-    fe_order = "integer"                  ## this is specific for fem, must be generalized
-  ),
-  methods = c(
-    get_basis = function() { return(basis_function) }
-  )
+#' @export 
+.FunctionSpaceCtr <- R6Class("FunctionSpace",
+                              public = list(
+                                mesh = "ANY",
+                                basis_function = "BasisFunction",     ## cpp backend
+                                fe_order = "integer",                  ## this is specific for fem, must be generalized
+                                initialize = function(mesh, basis_function, fe_order){
+                                  self$mesh <- mesh 
+                                  self$basis_function <- basis_function
+                                  self$fe_order <- fe_order
+                                },
+                                get_basis = function() { return(self$basis_function) }
+                              )
 )
 
+#' Finite Element Functional Space
+#'
+#' @name FunctionSpace
+#'
+#' @exportClass FunctionSpace
+setOldClass(c("FunctionSpace", "R6"))
+
 #' Create FunctionSpace object
 #'
 #' @param mesh A mesh object created by \code{Mesh}:
 #' @param fe_order Either '1' or '2'. It specifies the finite element order.
 #' @return An S4 object representing a Function Space.
 #' @export 
 #' @rdname FunctionSpace
-setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpace"))
+setGeneric("FunctionSpace", function(mesh, fe_order) standardGeneric("FunctionSpace"))
 
 #' @rdname FunctionSpace
 #' @examples
@@ -62,11 +81,11 @@ setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpa
 #' Vh <- FunctionSpace(mesh = mesh, fe_order = 1)
 #' }
 setMethod("FunctionSpace",
-          signature = c(mesh = "ANY", fe_order = "numeric"),
+          signature = c(mesh = "Mesh", fe_order = "numeric"),
           function(mesh, fe_order) {
-            .FunctionSpaceCtr(
-              basis_function = BasisFunction(mesh,as.integer(fe_order)),
+            .FunctionSpaceCtr$new(
               mesh = mesh,
+              basis_function = BasisFunction(mesh, as.integer(fe_order)),
               fe_order = as.integer(fe_order)
             )
           }
@@ -80,46 +99,55 @@ setMethod("FunctionSpace",
 #' mesh <- Mesh(unit_square)
 #' Vh <- FunctionSpace(mesh)
 #' }
-setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
+setMethod("FunctionSpace", signature = c(mesh="Mesh", fe_order="missing"),
           function(mesh){
-              return(.FunctionSpaceCtr(mesh=mesh, fe_order=1L))
+              return(.FunctionSpaceCtr$new(mesh=mesh, 
+                                           basis_function = BasisFunction(mesh, 1L), 
+                                           fe_order=1L))
 })
 
 ## finite element function
-.FunctionCtr <- setRefClass(
-  Class = "Function",
-  fields = c(
-    FunctionSpace  = "ANY", 
-    coeff = "matrix",
-    pde = "ANY"
-  ),
-  methods = list(
-    eval = function(X) {
-      M = dim(FunctionSpace$mesh$get_nodes())[2]
-      if(is.vector(X)) X <- t(as.matrix(X)) 
-
-      if(dim(X)[2] != M) {
-        stop(paste("matrix of evaluation points should be a 2 columns matrix"))
-      }
-      evals <- NULL  
-      if(length(FunctionSpace$mesh$times) == 0){
-        evals <- apply(FunctionSpace$basis_function$eval(as.matrix(X)) %*% coeff, 
-                       MARGIN=1, FUN=sum)
-      }else{
-        evals <- matrix(nrow=nrow(X),ncol=length(FunctionSpace$mesh$times))
-        for(t in 1:length(FunctionSpace$mesh$times)){
-          evals[,t] <- apply(FunctionSpace$basis_function$eval(as.matrix(X)) %*% coeff[,t], 
-                             MARGIN=1, FUN=sum) 
-        }  
-      }
-      return(evals)
-    },
-    set_coeff = function(coefficients){
-      coeff <<- coefficients
-    }
-  )
+.FunctionCtr <- R6Class("Function",
+                        public = list(
+                          FunctionSpace  = "FunctionSpace", 
+                          coeff = "matrix",
+                          initialize = function(FunctionSpace, coeff){
+                            self$FunctionSpace <- FunctionSpace
+                            self$coeff <- coeff
+                          },
+                          eval = function(X) {
+                            M = dim(self$FunctionSpace$mesh$get_nodes())[2]
+                            if(is.vector(X)) X <- t(as.matrix(X)) 
+                            
+                            if(dim(X)[2] != M) {
+                              stop(paste("matrix of evaluation points should be a 2 columns matrix"))
+                            }
+                            evals <- NULL  
+                            if(length(self$FunctionSpace$mesh$times) == 0){
+                              evals <- apply(self$FunctionSpace$basis_function$eval(as.matrix(X)) %*% self$coeff, 
+                                             MARGIN=1, FUN=sum)
+                            }else{
+                              evals <- matrix(nrow=nrow(X),ncol=length(self$FunctionSpace$mesh$times))
+                              for(t in 1:length(self$FunctionSpace$mesh$times)){
+                                evals[,t] <- apply(self$FunctionSpace$basis_function$eval(as.matrix(X)) %*% self$coeff[,t], 
+                                                   MARGIN=1, FUN=sum) 
+                              }  
+                            }
+                            return(evals)
+                          },
+                          set_coeff = function(coefficients){
+                            self$coeff <- coefficients
+                          }
+                        )
 )
 
+#' Finite Element Function
+#'
+#' @name Function
+#'
+#' @exportClass Function
+setOldClass(c("Function", "R6"))
+
 ## constructor
 
 #' Create Function object
@@ -137,7 +165,7 @@ setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
 #' }
 Function <- function(FunctionSpace) {
   coeff = matrix(ncol = 1, nrow = 0)
-  .FunctionCtr(coeff = coeff, FunctionSpace = FunctionSpace)
+  .FunctionCtr$new(coeff = coeff, FunctionSpace = FunctionSpace)
 }
 
 ## Function plot overload