From 3d329b6af1e4b06a52c7313f02c090a54f492380 Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Wed, 3 Jan 2024 14:32:44 +0100 Subject: [PATCH] femR relies on R6 classes. --- DESCRIPTION | 2 +- NAMESPACE | 22 ++- R/domain.R | 44 +++-- R/function_space.R | 174 ++++++++++------- R/mesh.R | 76 ++++---- R/operators.R | 190 ++++++++++++------- R/pde.R | 151 ++++++++------- man/div.Rd | 28 --- man/dot_product.Rd | 33 ---- man/dt-Function-missing-missing-method.Rd | 31 --- man/laplace.Rd | 27 --- man/plus-DiffOpObject-DiffOpObject-method.Rd | 19 -- man/times-matrix-FunctionGrad-method.Rd | 19 -- man/times-numeric-DiffOpObject-method.Rd | 19 -- man/times-numeric-Function-method.Rd | 29 --- tests/adv.diff.K.R | 8 +- tests/adv.diff.R | 8 +- tests/parabolic.R | 11 +- 18 files changed, 408 insertions(+), 483 deletions(-) delete mode 100644 man/div.Rd delete mode 100644 man/dot_product.Rd delete mode 100644 man/dt-Function-missing-missing-method.Rd delete mode 100644 man/laplace.Rd delete mode 100644 man/plus-DiffOpObject-DiffOpObject-method.Rd delete mode 100644 man/times-matrix-FunctionGrad-method.Rd delete mode 100644 man/times-numeric-DiffOpObject-method.Rd delete mode 100644 man/times-numeric-Function-method.Rd diff --git a/DESCRIPTION b/DESCRIPTION index cd786a0..c6e8aed 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,7 +15,7 @@ Authors@R: c( person("Aldo", "Clemente", role = c("aut", "cre"), email = "laura.sangalli@polimi.it")) Maintainer: Aldo Clemente 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 index 106ec36..4299a24 100644 --- 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 index a16e768..c69fcc4 100644 --- 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 index d5b228c..9f7174f 100644 --- a/R/function_space.R +++ b/R/function_space.R @@ -1,49 +1,68 @@ - -.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}: @@ -51,7 +70,7 @@ setMethod("BasisFunction", signature = signature("Mesh","integer"), #' @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 diff --git a/R/mesh.R b/R/mesh.R index 2084944..cb94265 100644 --- a/R/mesh.R +++ b/R/mesh.R @@ -1,37 +1,40 @@ ## Mesh Class Definition -.MeshCtr <- setRefClass( - Class = "Mesh", contains = "Domain", - fields = c( - cpp_handler = "ANY", ## cpp backend - m = "integer", ## local dim - n = "integer", ## embedding dim - times = "numeric", - deltaT = "numeric" # for parabolic problems - ), - methods = c( - get_nodes = function(){ - cpp_handler$nodes() - }, - get_elements = function(){ - cpp_handler$elements()+1 - }, - get_boundary = function(){ - cpp_handler$boundary() - }, - get_neighbors = function(){ - cpp_handler$neighbors() - }, - get_times = function(){ - return(times) - }, - set_deltaT = function(deltaT){ - if(length(time_interval)==0) - stop("Error!") - deltaT <<- deltaT - times <<- seq(time_interval[1], time_interval[2], by=deltaT) - } - ) +#' @export +.MeshCtr <- R6Class("Mesh", inherit = .DomainCtr, + public = list( + cpp_handler = "ANY", ## cpp backend + m = 0L, ## local dim + n = 0L, ## embedding dim + times = vector(mode="double", length = 0L), + deltaT = NA, # for parabolic problems + initialize = function(cpp_handler, m, n){ + self$cpp_handler <- cpp_handler + self$m <- m + self$n <- n + }, + get_nodes = function(){ + self$cpp_handler$nodes() + }, + get_elements = function(){ + self$cpp_handler$elements()+1 + }, + get_boundary = function(){ + self$cpp_handler$boundary() + }, + get_neighbors = function(){ + self$cpp_handler$neighbors() + }, + get_times = function(){ + return(self$times) + }, + set_deltaT = function(deltaT){ + if(length(self$time_interval)==0) + stop("Error!") + self$deltaT <- deltaT + self$times <- seq(self$time_interval[1], self$time_interval[2], by=deltaT) + } + ) ) #' Create mesh object @@ -67,13 +70,12 @@ setMethod("Mesh", signature = c(domain="list"), m <- ncol(domain$elements) - 1 n <- ncol(domain$nodes) if(m == 2 & n == 2) - .MeshCtr(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n), - times=vector(mode="double")) + .MeshCtr$new(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n)) else stop("wrong input argument provided.") }) -# @importClassesFrom RTriangle triangulation +setOldClass("triangulation") #' @rdname Mesh setMethod("Mesh", signature=c(domain="triangulation"), @@ -91,9 +93,7 @@ setMethod("Mesh", signature=c(domain="triangulation"), domain <- list(elements = elements, nodes = nodes, boundary = boundary) if(m == 2 & n == 2) - .MeshCtr(cpp_handler=new(Mesh_2D, domain), m=as.integer(m),n=as.integer(n), - times=vector(mode="double"), - time_interval = vector(mode="numeric", length = 0), crs = NA_crs_) + .MeshCtr$new(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n)) else stop("wrong input argument provided.") }) diff --git a/R/operators.R b/R/operators.R index 85a2964..82059bf 100644 --- a/R/operators.R +++ b/R/operators.R @@ -1,12 +1,22 @@ # gradient of Function -.FunctionGradCtr <- setRefClass( - Class = "FunctionGrad", - fields = c( +.FunctionGradCtr <- R6Class("FunctionGrad", + public = list( f = "Function", ## Function of which the gradient is taken - K = "ANY" ## set by product operator + K = "ANY", ## set by product operator + initialize = function(f, K=matrix(c(1,0,0,1), nrow=2,ncol=2,byrow=T)){ + self$f <- f + self$K <- K + } ) ) +#' Function gradient +#' +#' @name grad +#' +#' @exportClass FunctionGrad +setOldClass(c("FunctionGrad", "R6")) + ## take gradient of Function #' compute gradient of Function @@ -28,31 +38,30 @@ setGeneric("grad", function(f) standardGeneric("grad")) #' @rdname grad setMethod("grad", signature(f = "Function"), function(f) { - .FunctionGradCtr(f = f) + .FunctionGradCtr$new(f = f) }) +# setMethod("*", signature=c(e1="matrix", e2="FunctionGrad") + #' product overload for FunctionGradObejct #' #' @param e1 a numeric matrix. -#' @param e2 a FunctionGrad created by \code{grad} function. -#' @return A FunctionGrad. +#' @param e2 a FunctionGrad class created by \code{grad} function. +#' @return A FunctionGrad class. #' @export -setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"), - function(e1,e2){ - .FunctionGradCtr(f = e2$f, K = e1) -}) - -# ## diffusion tensor - FunctionGrad product overload -# `*.FunctionGrad` <- function(op1, op2) { -# if (!is.matrix(op1) && !is.function(op1)) -# stop("bad diffusion tensor type") -# .FunctionGradCtr(f = op2$f, K = op1) -# } +`*.FunctionGrad` <- function(e1, e2){ + if (!is.matrix(e1)) stop("First argument must be a matrix!") + .FunctionGradCtr$new(f = e2$f, K = e1) +} ## base class for differential operators -.DiffOpCtr <- setRefClass( - Class = "DiffOpObject", - fields = c( +.DiffOpCtr <- R6Class("DiffOpObject", + public = list( + initialize = function(tokens, params, f){ + self$tokens <- tokens + self$params <- params + self$f <- f + }, ## single blocks composing the overall operator tokens = "vector", params = "list", @@ -61,6 +70,15 @@ setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"), ) ) +#' R6 class representing a differential operator +#' +#' @name Differential Operator +#' +#' @exportClass DiffOpObject +setOldClass(c("DiffOpObject", "R6")) + +##setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"), + ## sum of differential operators #' plus operator overload for DiffOpObject @@ -68,23 +86,24 @@ setMethod("*", signature=c(e1="matrix", e2="FunctionGrad"), #' @param e1 a DiffOpObject. #' @param e2 a DiffOpObject. #' @return A S4 object representing the sum of two differential operators. +#' @name DifferentialOps #' @export -setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"), - function(e1, e2) { - if (tracemem(e1$f) != tracemem(e2$f)) { +`+.DiffOpObject` <- function(e1, e2){ + if (!identical(e1$f, e2$f)) { stop("operator arguments must be the same") } ## check for duplicated operator tokens if (any(duplicated(c(e1$tokens, e2$tokens)))) { stop("duplicated operator") } - .DiffOpCtr( + .DiffOpCtr$new( tokens = c(e1$tokens, e2$tokens), params = c(e1$params, e2$params), f = e1$f ) - } -) +} + +# setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"), #' difference of differential operators #' @@ -92,10 +111,13 @@ setMethod("+", signature = c(e1="DiffOpObject", e2="DiffOpObject"), #' @param e2 a DiffOpObject. #' @return A S4 object representing the difference of two differential operators. #' @rdname minus_DiffOb_op +#' @name DifferentialOps #' @export -setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"), - function(e1, e2) { - if (tracemem(e1$f) != tracemem(e2$f)) { +`-.DiffOpObject` <- function(e1, e2){ + if(missing(e2)){ + e1$params[[1]] <- -e1$params[[1]] + return(e1) + }else if (!identical(e1$f, e2$f)) { stop("operator arguments must be the same") } ## check for duplicated operator tokens @@ -105,55 +127,62 @@ setMethod("-", signature = c(e1="DiffOpObject", e2="DiffOpObject"), for(i in 1:length(e2$params)){ e2$params[[i]] = -e2$params[[i]] } - .DiffOpCtr( + .DiffOpCtr$new( tokens = c(e1$tokens, e2$tokens), params = c(e1$params, e2$params), f = e1$f ) -}) + } + +#setMethod("-", signature =c(e1 = "DiffOpObject", e2 = "missing"), - -# minus (unary) operator for DiffOpObject +# minus (unary) operator for differential operators # # @param e1 a DiffOpObject. # @return A S4 object differential operator whose paramter has changed sign. # @export -#' @rdname minus_DiffOb_op -setMethod("-", signature =c(e1 = "DiffOpObject", e2 = "missing"), - function(e1, e2){ - e1$params[[1]] <- -e1$params[[1]] - e1 - } -) +# @rdname minus_DiffOb_op +# `-.DiffOpObject` <- function(e1){ +# e1$params[[1]] <- -e1$params[[1]] +# e1 +# } + +# setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"), #' product by scalar for differential operators #' #' @param e1 a numeric. #' @param e2 a DiffOpObject. #' @return A S4 object representing the product by scalar for a differential operator. +#' @name DifferentialOps #' @export -setMethod("*", signature=c(e1="numeric", e2="DiffOpObject"), - function(e1,e2){ - #if (!is.numeric(e1)) stop("bad product") +`*.DiffOpObject` <- function(e1,e2){ + if (!is.numeric(e1)) stop("bad product") e2$params[[1]] <- e1*e2$params[[1]] e2 - -}) - + } + ## diffusion term -.DiffusionCtr <- setRefClass( - Class = "DiffusionOperator", - contains = "DiffOpObject" +.DiffusionCtr <- R6Class("DiffusionOperator", + inherit = .DiffOpCtr ) +#' Diffusion Operator +#' +#' @name DifferentialOperator +#' +#' @exportClass DiffusionOperator +setOldClass(c("DiffusionOperator", "DiffOpObject")) + ## laplace() returns a special operator for the case of ## isotropic and stationary diffusion -#' laplace operator for Function +#' laplace operator for Function class #' #' @param f a Function. #' @return A S4 object representing the laplace operator applied to the function passed as parameter. +#' @rdname DiffusionOperator #' @export #' @examples #' \dontrun{ @@ -168,7 +197,7 @@ laplace <- function(f) { if (!is(f, "Function")) { stop("wrong argument type") } - .DiffusionCtr( + .DiffusionCtr$new( tokens = "diffusion", params = list(diffusion = 1), f = f @@ -181,6 +210,7 @@ laplace <- function(f) { #' #' @param f a Function. #' @return A S4 object representing the diffusion term of a second order linear differential operator. +#' @rdname DiffusionOperator #' @export #' @examples #' \dontrun{ @@ -195,7 +225,7 @@ laplace <- function(f) { div <- function(f) { if (is(f, "FunctionGrad")) { if (!is.null(f$K)) { - return(.DiffusionCtr( + return(.DiffusionCtr$new( tokens = "diffusion", params = list(diffusion = f$K), f = f$f) @@ -206,17 +236,23 @@ div <- function(f) { } ## transport term -.TransportCtr <- setRefClass( - Class = "TransportOperator", - contains = "DiffOpObject" +.TransportCtr <- R6Class("TransportOperator", + inherit = .DiffOpCtr ) +#' Transport Operator +#' +#' @name TransportOperator +#' +#' @exportClass TransportOperator +setOldClass(c("TransportOperator", "DiffOpObject")) + #' dot product between vector and FunctionGrad #' #' @param op1 a numeric vector. #' @param op2 a FunctionGrad. #' @return A S4 object representing the advection term of a second order linear differential operator. -#' @rdname dot_product +#' @rdname TransportOperator #' @export #' @examples #' \dontrun{ @@ -230,10 +266,10 @@ div <- function(f) { #' } setGeneric("dot", function(op1, op2) standardGeneric("dot")) -#' @rdname dot_product +#' @rdname TransportOperator setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"), function(op1, op2) { - .TransportCtr( + .TransportCtr$new( tokens = "transport", params = list(transport = as.matrix(op1)), f = op2$f @@ -241,16 +277,25 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"), }) ## reaction term -.ReactionCtr <- setRefClass( - Class = "ReactionOperator", - contains = "DiffOpObject" +.ReactionCtr <- R6Class("ReactionOperator", + inherit = .DiffOpCtr ) +#' Reaction Operator +#' +#' @name ReactionOperator +#' +#' @exportClass ReactionOperator +setOldClass(c("ReactionOperator", "DiffOpObject")) + +# setMethod("*", signature = c(e1="numeric", e2="Function"), + #' product by scalar for FunctionObejct #' #' @param e1 a numeric. #' @param e2 a FunctioObject created by \code{Function}. #' @return A S4 object representing the reaction term of a second order linear differential operator. +#' @rdname ReactionOperator #' @export #' @examples #' \dontrun{ @@ -261,20 +306,24 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"), #' f <- Function(Vh) #' reaction <- 2*f #' } -setMethod("*", signature = c(e1="numeric", e2="Function"), - function(e1,e2){ - .ReactionCtr( +`*.Function` <- function(e1,e2){ + if(!is(e1,"numeric")) stop("First argument must be a scalar!") + .ReactionCtr$new( tokens = "reaction", params = list(reaction = e1), f = e2 ) -}) +} -.TimeDerivativeCtr <- setRefClass( - Class = "TimeDerivative", - contains ="DiffOpObject" +.TimeDerivativeCtr <- R6Class("TimeDerivative", + inherit = .DiffOpCtr ) +#' @name TimeDerivative +#' +#' @exportClass TimeDerivative +setOldClass(c("TimeDerivative", "DiffOpObject")) + # overloading stats::dt function #' time derivate of FunctionObejct @@ -283,6 +332,7 @@ setMethod("*", signature = c(e1="numeric", e2="Function"), #' @param df missing. #' @param ncp missing. #' @return A S4 object representing the time derivative of a Function. +#' @rdname TimeDerivative #' @export #' @examples #' \dontrun{ @@ -295,7 +345,7 @@ setMethod("*", signature = c(e1="numeric", e2="Function"), #' } setMethod("dt", signature = c(x="Function", df="missing", ncp="missing"), function(x,df,ncp){ - .TimeDerivativeCtr( + .TimeDerivativeCtr$new( tokens="time", params = list(time=1L), f=x diff --git a/R/pde.R b/R/pde.R index 877f7bd..178096c 100644 --- a/R/pde.R +++ b/R/pde.R @@ -2,123 +2,141 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3) # Pde Class wraps C++ R_PDE class -.PdeCtr <- setRefClass( - Class = "PdeObject", - fields = c( +.PdeCtr <- R6Class("Pde", + public = list( is_dirichletBC_set = "logical", is_initialCondition_set = "logical", is_parabolic ="logical", - cpp_handler = "ANY"), - methods = c( + cpp_handler = "ANY", + DifferentialOp = "DiffOpObject", + initialize = function(cpp_handler, DifferentialOp, is_parabolic, is_dirichletBC_set, is_initialCondition_set){ + self$cpp_handler <- cpp_handler + self$DifferentialOp <- DifferentialOp + self$is_parabolic <- is_parabolic + self$is_dirichletBC_set <- is_dirichletBC_set + self$is_initialCondition_set <- is_initialCondition_set + }, solve = function(){ - if(!is_dirichletBC_set){ - if(!is_parabolic){ - dirichletBC_ = as.matrix(rep(0,times=nrow(cpp_handler$get_dofs_coordinates()))) + if(!self$is_dirichletBC_set){ + if(!self$is_parabolic){ + dirichletBC_ = as.matrix(rep(0,times=nrow(self$cpp_handler$get_dofs_coordinates()))) }else - dirichletBC_ = matrix(0, nrow=nrow(cpp_handler$get_dofs_coordinates()), ncol=length(times)) - cpp_handler$set_dirichlet_bc(dirichletBC_) - is_dirichletBC_set <<- TRUE + dirichletBC_ = matrix(0, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), + ncol=length(self$DifferentialOp$f$FunctionSpace$mesh$times)) + self$cpp_handler$set_dirichlet_bc(dirichletBC_) + self$is_dirichletBC_set <- TRUE } - if(is_parabolic & (!is_initialCondition_set)){ + if(self$is_parabolic & (!self$is_initialCondition_set)){ stop("initialCondition must be provided.") } - cpp_handler$solve() - }, - get_solution = function(){ - cpp_handler$solution() + self$cpp_handler$solve() }, get_dofs_coordinates = function(){ - cpp_handler$get_dofs_coordinates() + self$cpp_handler$get_dofs_coordinates() }, set_boundary_condition = function(fun, type="dirichlet", on=NULL){ if(!any(type == c("dirichlet", "Dirichlet", "d"))) stop("Only Dirichlet boundary condtions allowed.") dirichletBC_ <- NULL - if(typeof(fun) == "closure"){ - if(!is_parabolic){ - dirichletBC_ <- as.matrix(fun(cpp_handler$get_dofs_coordinates())) + if(any(typeof(fun) == c("function", "closure"))){ + if(!self$is_parabolic){ + dirichletBC_ <- as.matrix(fun(self$cpp_handler$get_dofs_coordinates())) }else{ - dirichletBC_ <- fun(cpp_handler$get_dofs_coordinates(),times) + dirichletBC_ <- fun(self$cpp_handler$get_dofs_coordinates(), + self$DifferentialOp$f$FunctionSpace$mesh$times) } }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){ - if(nrow(as.matrix(fun)) == nrow(cpp_handler$get_dofs_coordinates())){ + if(nrow(as.matrix(fun)) == nrow(self$cpp_handler$get_dofs_coordinates())){ dirichletBC_ <- fun }else if (nrow(as.matrix(fun)) == 1L){ - if(!is_parabolic) - dirichletBC_ <- matrix(fun, nrow=nrow(cpp_handler$get_dofs_coordinates()), ncol=1) + if(!self$is_parabolic) + dirichletBC_ <- matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), ncol=1) else - dirichletBC_ <- matrix(fun, nrow=nrow(cpp_handler$get_dofs_coordinates()), ncol=length(times)) + dirichletBC_ <- matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), + ncol=length(self$DifferentialOp$f$FunctionSpace$mesh$times)) } } - is_dirichletBC_set <<- TRUE - cpp_handler$set_dirichlet_bc(dirichletBC_) + self$is_dirichletBC_set <- TRUE + self$cpp_handler$set_dirichlet_bc(dirichletBC_) }, set_initial_condition = function(fun){ - if(!is_parabolic) + if(!self$is_parabolic) stop("Cannot set initial condition for elliptic problem.") - is_initialCondition_set <<- TRUE + self$is_initialCondition_set <- TRUE if(typeof(fun) == "closure"){ - cpp_handler$set_initial_condition(fun(cpp_handler$get_dofs_coordinates())) + self$cpp_handler$set_initial_condition(fun(self$cpp_handler$get_dofs_coordinates())) }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){ - if(nrow(as.matrix(fun)) == nrow(cpp_handler$get_dofs_coordinates())){ - cpp_handler$set_initial_condition(fun) + if(nrow(as.matrix(fun)) == nrow(self$cpp_handler$get_dofs_coordinates())){ + self$cpp_handler$set_initial_condition(fun) }else if(nrow(as.matrix(fun)) == 1L){ - cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(cpp_handler$get_dofs_coordinates()), + self$cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), ncol=1)) } } }, get_mass = function(){ - cpp_handler$mass() + self$cpp_handler$mass() }, get_stiff = function(){ - cpp_handler$stiff() + self$cpp_handler$stiff() } ) ) +#' @name pde +#' +#' @exportClass Pde +setOldClass(c("Pde", "R6")) + ## infers the type of a pde -extract_pde_type <- function(L) { - if ("time" %in% names(L$params)) { +extract_pde_type <- function(DifferentialOp) { + if ("time" %in% names(DifferentialOp$params)) { return(pde_type_list$parabolic) } - if ("diffusion" %in% names(L$params) && !is.matrix(L$params$diffusion)) { + if ("diffusion" %in% names(DifferentialOp$params) && !is.matrix(DifferentialOp$params$diffusion)) { return(pde_type_list$laplacian) } return(pde_type_list$elliptic) } ## parse pde parameters -parse_pde_parameters <- function(L){ +parse_pde_parameters <- function(DifferentialOp){ + pde_type <- extract_pde_type(DifferentialOp) + pde_parameters <- NULL pde_parameters$diffusion <- 1.0 pde_parameters$transport <- matrix(0, nrow = 2, ncol = 1) pde_parameters$reaction <- 0.0 - for (i in length(L$params)) { - pde_parameters[[paste(L$tokens[i])]] <- L$params[[paste(L$tokens[i])]] + for (i in 1:length(DifferentialOp$params)) { + pde_parameters[[DifferentialOp$tokens[i]]] <- DifferentialOp$params[[DifferentialOp$tokens[i]]] } + + if(pde_type == pde_type_list$parabolic) + pde_parameters[["time_mesh"]] <- as.vector(DifferentialOp$f$FunctionSpace$mesh$times) + + if(pde_type == pde_type_list$parabolic & is(pde_parameters[["diffusion"]], "numeric")) + pde_parameters[["diffusion"]] <- pde_parameters[["diffusion"]]*matrix(c(1,0,0,1), nrow=2, ncol=2, byrow=T) + return(pde_parameters) } ## build cpp object -make_pde <- function(L) { - pde_type = extract_pde_type(L) +make_pde <- function(DifferentialOp) { + pde_type <- extract_pde_type(DifferentialOp) ## pde parameters - pde_parameters <- parse_pde_parameters(L) - - if(pde_type == pde_type_list$parabolic) pde_parameters["time_mesh"] <- L$f$FunctionSpace$mesh$times + pde_parameters <- parse_pde_parameters(DifferentialOp) ## define Rcpp module - D <- L$f$FunctionSpace$mesh$cpp_handler ## domain - fe_order <- L$f$FunctionSpace$fe_order ## finite element order + D <- DifferentialOp$f$FunctionSpace$mesh$cpp_handler ## domain + fe_order <- DifferentialOp$f$FunctionSpace$fe_order ## finite element order cpp_handler <- NULL if (fe_order == 1) { ## linear finite elements - cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters, L$f) + cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters, DifferentialOp$f) } if (fe_order == 2) { ## quadratic finite elements - cpp_handler <- new(cpp_pde_2d_fe2, D, pde_type - 1, pde_parameters, L$f) + cpp_handler <- new(cpp_pde_2d_fe2, D, pde_type - 1, pde_parameters, DifferentialOp$f) } return(cpp_handler) @@ -126,24 +144,24 @@ make_pde <- function(L) { #' A PDEs object #' -#' @param L a differential operator. +#' @param DifferentialOp a differential operator. #' @param f a standard R function representing the forcing term of the PDE or a numeric value, in case of constant forcing term. -#' @return A S4 object representing a PDE. +#' @return A R6 class representing a PDE. #' @rdname pde #' @export -setGeneric("Pde", function(L,f) standardGeneric("Pde")) +setGeneric("Pde", function(DifferentialOp,f) standardGeneric("Pde")) #' @rdname pde -setMethod("Pde", signature=c(L="DiffOpObject", f="ANY"), - function(L,f){ +setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"), + function(DifferentialOp,f){ - pde_type = extract_pde_type(L) - cpp_handler <- make_pde(L) + pde_type = extract_pde_type(DifferentialOp) + cpp_handler <- make_pde(DifferentialOp) quad_nodes <- cpp_handler$get_quadrature_nodes() ## evaluate forcing term on quadrature nodes if(pde_type == pde_type_list$parabolic) { - times <- L$f$FunctionSpace$mesh$times + times <- DifferentialOp$f$FunctionSpace$mesh$times cpp_handler$set_forcing(as.matrix(f(quad_nodes, times))) }else{ cpp_handler$set_forcing(as.matrix(f(quad_nodes))) @@ -155,16 +173,17 @@ setMethod("Pde", signature=c(L="DiffOpObject", f="ANY"), is_parabolic = pde_type == pde_type_list$parabolic ## return - .PdeCtr(is_dirichletBC_set = FALSE, - is_initialCondition_set = FALSE, - is_parabolic = is_parabolic, - cpp_handler = cpp_handler) + .PdeCtr$new(is_dirichletBC_set = FALSE, + is_initialCondition_set = FALSE, + is_parabolic = is_parabolic, + cpp_handler = cpp_handler, + DifferentialOp = DifferentialOp) }) #' @rdname pde -setMethod("Pde", signature=c(L="DiffOpObject", f="numeric"), - function(L,f){ - pde_type = extract_pde_type(L) +setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="numeric"), + function(DifferentialOp,f){ + pde_type = extract_pde_type(DifferentialOp) fun <- NULL if(pde_type!=pde_type_list$parabolic){ @@ -177,5 +196,5 @@ setMethod("Pde", signature=c(L="DiffOpObject", f="numeric"), } } - return(Pde(L,fun)) + return(Pde(DifferentialOp,fun)) }) diff --git a/man/div.Rd b/man/div.Rd deleted file mode 100644 index d679390..0000000 --- a/man/div.Rd +++ /dev/null @@ -1,28 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{div} -\alias{div} -\title{divergence operator FunctionGrad} -\usage{ -div(f) -} -\arguments{ -\item{f}{a Function.} -} -\value{ -A S4 object representing the diffusion term of a second order linear differential operator. -} -\description{ -divergence operator FunctionGrad -} -\examples{ -\dontrun{ -library(femR) -data("unit_square") -mesh <- Mesh(unit_square) -Vh <- FunctionSpace(mesh) -f <- Function(Vh) -K <- matrix(c(1,2,1,0),nrow=2,ncol=2) -diffusion <- div(K*grad(f)) -} -} diff --git a/man/dot_product.Rd b/man/dot_product.Rd deleted file mode 100644 index 30d7bb7..0000000 --- a/man/dot_product.Rd +++ /dev/null @@ -1,33 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{dot} -\alias{dot} -\alias{dot,vector,FunctionGrad-method} -\title{dot product between vector and FunctionGrad} -\usage{ -dot(op1, op2) - -\S4method{dot}{vector,FunctionGrad}(op1, op2) -} -\arguments{ -\item{op1}{a numeric vector.} - -\item{op2}{a FunctionGrad.} -} -\value{ -A S4 object representing the advection term of a second order linear differential operator. -} -\description{ -dot product between vector and FunctionGrad -} -\examples{ -\dontrun{ -library(femR) -data("unit_square") -mesh <- Mesh(unit_square) -Vh <- FunctionSpace(mesh) -f <- Function(Vh) -b <- c(1,1) -advection <- dot(b,grad(f)) -} -} diff --git a/man/dt-Function-missing-missing-method.Rd b/man/dt-Function-missing-missing-method.Rd deleted file mode 100644 index f9be2af..0000000 --- a/man/dt-Function-missing-missing-method.Rd +++ /dev/null @@ -1,31 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{dt,Function,missing,missing-method} -\alias{dt,Function,missing,missing-method} -\title{time derivate of FunctionObejct} -\usage{ -\S4method{dt}{Function,missing,missing}(x, df, ncp) -} -\arguments{ -\item{x}{a Function.} - -\item{df}{missing.} - -\item{ncp}{missing.} -} -\value{ -A S4 object representing the time derivative of a Function. -} -\description{ -time derivate of FunctionObejct -} -\examples{ -\dontrun{ -library(femR) -data("unit_square") -mesh <- Mesh(unit_square) -Vh <- FunctionSpace(mesh) -f <- Function(Vh) -dt(f) -} -} diff --git a/man/laplace.Rd b/man/laplace.Rd deleted file mode 100644 index 5c6e12f..0000000 --- a/man/laplace.Rd +++ /dev/null @@ -1,27 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{laplace} -\alias{laplace} -\title{laplace operator for Function} -\usage{ -laplace(f) -} -\arguments{ -\item{f}{a Function.} -} -\value{ -A S4 object representing the laplace operator applied to the function passed as parameter. -} -\description{ -laplace operator for Function -} -\examples{ -\dontrun{ -library(femR) -data("unit_square") -mesh <- Mesh(unit_square) -Vh <- FunctionSpace(mesh) -f <- Function(Vh) -laplace_f <- laplace(f) -} -} diff --git a/man/plus-DiffOpObject-DiffOpObject-method.Rd b/man/plus-DiffOpObject-DiffOpObject-method.Rd deleted file mode 100644 index 659e296..0000000 --- a/man/plus-DiffOpObject-DiffOpObject-method.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{+,DiffOpObject,DiffOpObject-method} -\alias{+,DiffOpObject,DiffOpObject-method} -\title{plus operator overload for DiffOpObject} -\usage{ -\S4method{+}{DiffOpObject,DiffOpObject}(e1, e2) -} -\arguments{ -\item{e1}{a DiffOpObject.} - -\item{e2}{a DiffOpObject.} -} -\value{ -A S4 object representing the sum of two differential operators. -} -\description{ -plus operator overload for DiffOpObject -} diff --git a/man/times-matrix-FunctionGrad-method.Rd b/man/times-matrix-FunctionGrad-method.Rd deleted file mode 100644 index 1033a27..0000000 --- a/man/times-matrix-FunctionGrad-method.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{*,matrix,FunctionGrad-method} -\alias{*,matrix,FunctionGrad-method} -\title{product overload for FunctionGradObejct} -\usage{ -\S4method{*}{matrix,FunctionGrad}(e1, e2) -} -\arguments{ -\item{e1}{a numeric matrix.} - -\item{e2}{a FunctionGrad created by \code{grad} function.} -} -\value{ -A FunctionGrad. -} -\description{ -product overload for FunctionGradObejct -} diff --git a/man/times-numeric-DiffOpObject-method.Rd b/man/times-numeric-DiffOpObject-method.Rd deleted file mode 100644 index 9101d62..0000000 --- a/man/times-numeric-DiffOpObject-method.Rd +++ /dev/null @@ -1,19 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{*,numeric,DiffOpObject-method} -\alias{*,numeric,DiffOpObject-method} -\title{product by scalar for differential operators} -\usage{ -\S4method{*}{numeric,DiffOpObject}(e1, e2) -} -\arguments{ -\item{e1}{a numeric.} - -\item{e2}{a DiffOpObject.} -} -\value{ -A S4 object representing the product by scalar for a differential operator. -} -\description{ -product by scalar for differential operators -} diff --git a/man/times-numeric-Function-method.Rd b/man/times-numeric-Function-method.Rd deleted file mode 100644 index c56f218..0000000 --- a/man/times-numeric-Function-method.Rd +++ /dev/null @@ -1,29 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/operators.R -\name{*,numeric,Function-method} -\alias{*,numeric,Function-method} -\title{product by scalar for FunctionObejct} -\usage{ -\S4method{*}{numeric,Function}(e1, e2) -} -\arguments{ -\item{e1}{a numeric.} - -\item{e2}{a FunctioObject created by \code{Function}.} -} -\value{ -A S4 object representing the reaction term of a second order linear differential operator. -} -\description{ -product by scalar for FunctionObejct -} -\examples{ -\dontrun{ -library(femR) -data("unit_square") -mesh <- Mesh(unit_square) -Vh <- FunctionSpace(mesh) -f <- Function(Vh) -reaction <- 2*f -} -} diff --git a/tests/adv.diff.K.R b/tests/adv.diff.K.R index ee72fe1..8fabddd 100644 --- a/tests/adv.diff.K.R +++ b/tests/adv.diff.K.R @@ -34,7 +34,7 @@ f <- Function(Vh) #L <- -1*laplace(f) + dot(c(alpha_,0), grad(f)) # OK I <- matrix(c(1,0,0,1), byrow=T, nrow=2,ncol=2) L <- -div(I*grad(f)) + dot(c(alpha_,0), grad(f)) - +# L <- -div(grad(f)) + dot(c(alpha_,0), grad(f)) ## forcing term u <- function(points){ return(gamma_ * sin(pi * points[,2])) @@ -53,18 +53,18 @@ pde$solve() ## compute L2 norm of the error u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates())) -error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2)) +error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - f$coeff)^2)) cat("L2 error = ", error.L2, "\n") ## perform evaluation at single point point = c(0.2, 0.5) -f$eval_at(point) +f$eval(point) ## evaluate over a 10x10 grid x <- seq(0, 1, length.out = 50) y <- x points <- expand.grid(x, y) -max(abs( f$eval_at(points) - as.matrix(exact_solution(points)))) +max(abs( f$eval(points) - as.matrix(exact_solution(points)))) # plot solution options(warn=-1) diff --git a/tests/adv.diff.R b/tests/adv.diff.R index 717243b..e23d5e6 100644 --- a/tests/adv.diff.R +++ b/tests/adv.diff.R @@ -47,19 +47,19 @@ pde$solve() ## compute L2 norm of the error u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates())) -error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2)) +error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - f$coeff)^2)) cat("L2 error = ", error.L2, "\n") ## perform evaluation at single point point = c(0.2, 0.5) -f$eval_at(point) +f$eval(point) ## evaluate over a 10x10 grid x <- seq(0, 1, length.out = 50) y <- x points <- expand.grid(x, y) -f$eval_at(points) -max(abs( f$eval_at(points) - as.matrix(exact_solution(points)))) +f$eval(points) +max(abs( f$eval(points) - as.matrix(exact_solution(points)))) # plot solution options(warn=-1) diff --git a/tests/parabolic.R b/tests/parabolic.R index d0fc942..6ce7336 100644 --- a/tests/parabolic.R +++ b/tests/parabolic.R @@ -33,7 +33,7 @@ u <- Function(Vh) ## define differential operator in its strong formulation #L <- dt(f) + (-1)*laplace(f) -L <- dt(u) - laplace(u) + dot(c(0,10),grad(u)) +L <- dt(u) - laplace(u) + dot(c(0,1),grad(u)) ## forcing term f <- function(points,times){ @@ -60,6 +60,7 @@ pde$set_initial_condition(u0) pde$set_boundary_condition(g(pde$get_dofs_coordinates(),times)) pde$set_initial_condition(u0(pde$get_dofs_coordinates())) + pde$solve() ## compute L2 norm of the error @@ -67,20 +68,20 @@ u_ex <- exact_solution(pde$get_dofs_coordinates(), times) error.L2 <- matrix(0, nrow=length(times), ncol=1) for( t in 1:length(times)){ - error.L2[t] <- sqrt(sum(pde$get_mass() %*% (u_ex[,t] - pde$get_solution()[,t])^2)) + error.L2[t] <- sqrt(sum(pde$get_mass() %*% (u_ex[,t] - u$coeff[,t])^2)) cat(paste0("L2 error at time ",times[t]," ", error.L2[t], "\n")) } # ------------------------------------------------------------------------------ ## perform evaluation at single point point = c(0.2, 0.5) -u$eval_at(point) +u$eval(point) # # ## evaluate over a 10x10 grid x <- seq(0, 1, length.out = 50) y <- x points <- expand.grid(x, y) -dim( u$eval_at(points) ) -max(abs( u$eval_at(points) - as.matrix(exact_solution(points)))) +dim( u$eval(points) ) +max(abs( u$eval(points) - as.matrix(exact_solution(points)))) # # ## plot solution options(warn=-1)