From 252bcd000c12c74bc030fa334411785737ebdb69 Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Mon, 15 Jan 2024 10:34:36 +0100 Subject: [PATCH] Updated R interface. --- DESCRIPTION | 1 - NAMESPACE | 1 - R/datasets.R | 9 ++- R/domain.R | 44 +++++----- R/function_space.R | 134 +++++++++++++++++-------------- R/mesh.R | 91 ++++++++++----------- R/operators.R | 72 ++++++++++------- R/pde.R | 81 ++++++++----------- R/sf.R | 33 ++++---- R/utils.R | 32 +------- data/lake_como_boundary.RData | Bin 0 -> 11456 bytes src/include/r_functional_space.h | 2 +- src/include/r_pde.h | 8 +- src/r_functional_space.cpp | 16 ++-- src/r_pde.cpp | 8 +- tests/parabolic.R | 8 +- tests/pde.2.R | 16 ++-- vignettes/Introduction.Rmd | 16 ++-- vignettes/LakeComo.Rmd | 36 ++++----- 19 files changed, 297 insertions(+), 311 deletions(-) create mode 100644 data/lake_como_boundary.RData diff --git a/DESCRIPTION b/DESCRIPTION index c6e8aed..8635e69 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,6 +29,5 @@ Encoding: UTF-8 Suggests: knitr, rmarkdown, - rmapshaper, mapview VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 4299a24..84f46d6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,7 +14,6 @@ S3method(st_as_sfc,Mesh) S3method(st_crs,Domain) export("%X%") export(.FunctionSpaceCtr) -export(.MeshCtr) export(Domain) export(Function) export(FunctionSpace) diff --git a/R/datasets.R b/R/datasets.R index eb9382b..054c26d 100644 --- a/R/datasets.R +++ b/R/datasets.R @@ -5,4 +5,11 @@ #' This dataset can be used to create a mesh object with the function \code{Mesh}. #' #' @name unit_square -NULL \ No newline at end of file +NULL + +#' Lake Como boundary +#' +#' sf object storing the boundary of Lake Como. +#' +#' @name lake_como_boundary +NULL diff --git a/R/domain.R b/R/domain.R index 58c0686..02ff0a8 100644 --- a/R/domain.R +++ b/R/domain.R @@ -1,17 +1,17 @@ .DomainCtr <- R6Class("Domain", private = list( - geometry = "ANY", - coords = "matrix", # (x,y) - boundary = "matrix", # - time_interval = vector(mode="numeric", length = 0), - crs = "ANY" + geometry_ = "ANY", + coords_ = "matrix", # (x,y) + boundary_ = "matrix", # + time_interval_ = vector(mode="numeric", length = 0), + crs_ = "ANY" ), public = list( initialize = function(geometry, coords, boundary, crs){ - private$geometry <- geometry - private$coords <- coords - private$boundary <- boundary - private$crs <- crs + private$geometry_ <- geometry + private$coords_ <- coords + private$boundary_ <- boundary + private$crs_ <- crs } ) ) @@ -254,39 +254,39 @@ setMethod("build_mesh", signature=c("Domain", "numeric", "numeric"), SB = matrix(nrow=0, ncol=1) H = matrix(nrow=0, ncol=2) - for( sub_id in 1:length(extract_private(domain)$geometry)){ - groups_id <- unique(extract_private(domain)$geometry[[sub_id]]$edges_ring) + for( sub_id in 1:length(extract_private(domain)$geometry_)){ + groups_id <- unique(extract_private(domain)$geometry_[[sub_id]]$edges_ring) holes_id <- groups_id[ groups_id < 0 ] num_holes <- length(holes_id) if(num_holes > 0){ for(i in 1:num_holes){ - edges_id <- which(extract_private(domain)$geometry[[sub_id]]$edges_ring==holes_id[i]) - if( ! all(as.integer(extract_private(domain)$geometry[[sub_id]]$edges_boundary[edges_id])) ){ # NON VERO buco + edges_id <- which(extract_private(domain)$geometry_[[sub_id]]$edges_ring==holes_id[i]) + if( ! all(as.integer(extract_private(domain)$geometry_[[sub_id]]$edges_boundary[edges_id])) ){ # NON VERO buco next } - path <- t(extract_private(domain)$geometry[[sub_id]]$edges[edges_id,])[1,] + path <- t(extract_private(domain)$geometry_[[sub_id]]$edges[edges_id,])[1,] path <- c(path,path[1]) - nodes <- as.matrix(extract_private(domain)$coords[path,1:2]) + nodes <- as.matrix(extract_private(domain)$coords_[path,1:2]) holes <- st_point_on_surface(st_polygon(list(nodes))) H = rbind(H, holes) } } - S = rbind(S, extract_private(domain)$geometry[[sub_id]]$edges) - SB = rbind(SB, as.matrix(extract_private(domain)$geometry[[sub_id]]$edges_boundary)) + S = rbind(S, extract_private(domain)$geometry_[[sub_id]]$edges) + SB = rbind(SB, as.matrix(extract_private(domain)$geometry_[[sub_id]]$edges_boundary)) } if(nrow(H)==0){ H = NA } - pslg <- pslg(P = extract_private(domain)$coords[,1:2], PB = as.matrix(extract_private(domain)$boundary), + pslg <- pslg(P = extract_private(domain)$coords_[,1:2], PB = as.matrix(extract_private(domain)$boundary_), S = S, SB = SB, H = H) triangulation <- triangulate(p = pslg, a = maximum_area, q = minimum_angle) res <- Mesh(triangulation) - set_geometry(res, extract_private(domain)$geometry) - set_time_interval(res, extract_private(domain)$time_interval) - set_crs(res, extract_private(domain)$crs) + set_private(res, "geometry_" , extract_private(domain)$geometry_ ) + set_private(res, "time_interval_", extract_private(domain)$time_interval_) + set_private(res, "crs_" , extract_private(domain)$crs_ ) return(res) }) @@ -304,6 +304,6 @@ setMethod("%X%", signature=c(op1="Domain", op2="numeric"), function(op1, op2){ if(op2[1] > op2[length(op2)]) stop("Error! First time instant is greater than last time instant.") - set_time_interval(op1, op2) + set_private(op1, "time_interval_", op2) op1 }) \ No newline at end of file diff --git a/R/function_space.R b/R/function_space.R index b9ea828..b6401e5 100644 --- a/R/function_space.R +++ b/R/function_space.R @@ -1,22 +1,22 @@ .BasisFunctionCtr <- R6Class("BasisFunction", private = list( - cpp_handler = "ANY", ## pointer to cpp backend - mesh = "Mesh" + cpp_handler_ = "ANY", ## pointer to cpp backend + mesh_ = "Mesh" ), public = list( initialize = function(cpp_handler, mesh){ - private$cpp_handler = cpp_handler - private$mesh = mesh + private$cpp_handler_ = cpp_handler + private$mesh_ = mesh }, size = function(){ - private$cpp_handler$size() + private$cpp_handler_$size() }, - get_dofs_coordinates = function(){ - private$cpp_handler$get_dofs_coordinates() + dofs_coordinates = function(){ + private$cpp_handler_$dofs_coordinates() }, ## evaluates the basis system on the given set of locations eval = function(locations) { - return(private$cpp_handler$eval(0L, locations)) + return(private$cpp_handler_$eval(0L, locations)) }, ## integrates a function expressed as basis expansion with respect to this basis system over the ## whole domain of definition @@ -25,9 +25,9 @@ ## required dof coordinates... c <- func if (is.function(func)) { ## recover fem expansion - c <- func(private$cpp_handler$get_dofs_coordinates()) + c <- func(private$cpp_handler_$dofs_coordinates()) } - return(private$cpp_handler$integrate(c)) + return(private$cpp_handler_$integrate(c)) } ) ) @@ -43,27 +43,33 @@ setGeneric("BasisFunction", function(mesh, fe_order) standardGeneric("BasisFunct setMethod("BasisFunction", signature = signature("Mesh","integer"), function(mesh, fe_order){ - basis_function <- .BasisFunctionCtr$new(cpp_handler = + basis <- .BasisFunctionCtr$new(cpp_handler = new(eval(parse(text = paste("cpp_lagrange_basis_2d_fe", - as.character(fe_order), sep = ""))), extract_private(mesh)$cpp_handler, 0), + as.character(fe_order), sep = ""))), extract_private(mesh)$cpp_handler_, 0), mesh=mesh) }) #' @export .FunctionSpaceCtr <- R6Class("FunctionSpace", private= list( - mesh = "Mesh", - basis_function = "BasisFunction", ## cpp backend - fe_order = "integer" ## this is specific for fem, must be generalized + mesh_ = "Mesh", + basis_ = "BasisFunction", ## cpp backend + fe_order_ = "integer" ## this is specific for fem, must be generalized ), public = list( - initialize = function(mesh, basis_function, fe_order){ - private$mesh <- mesh - private$basis_function <- basis_function - private$fe_order <- fe_order + initialize = function(mesh, basis, fe_order){ + private$mesh_ <- mesh + private$basis_ <- basis + private$fe_order_ <- fe_order }, - get_basis = function(){ - return(private$basis_function) + basis = function(){ + return(private$basis_) + }, + mesh = function(){ + private$mesh_ + }, + fe_order = function(){ + private$fe_order_ } ) ) @@ -97,7 +103,7 @@ setMethod("FunctionSpace", function(mesh, fe_order) { .FunctionSpaceCtr$new( mesh = mesh, - basis_function = BasisFunction(mesh, as.integer(fe_order)), + basis = BasisFunction(mesh, as.integer(fe_order)), fe_order = as.integer(fe_order) ) } @@ -114,45 +120,51 @@ setMethod("FunctionSpace", setMethod("FunctionSpace", signature = c(mesh="Mesh", fe_order="missing"), function(mesh){ return(.FunctionSpaceCtr$new(mesh=mesh, - basis_function = BasisFunction(mesh, 1L), + basis = BasisFunction(mesh, 1L), fe_order=1L)) }) ## finite element function .FunctionCtr <- R6Class("Function", private = list( - FunctionSpace = "FunctionSpace" + FunctionSpace_ = "FunctionSpace", + coeff_ = "matrix" ), public = list( - coeff = "matrix", - initialize = function(FunctionSpace, coeff){ - private$FunctionSpace <- FunctionSpace - self$coeff <- coeff + initialize = function(FunctionSpace, coefficients){ + private$FunctionSpace_ <- FunctionSpace + private$coeff_ <- coefficients }, eval = function(X) { - M = dim(extract_private(private$FunctionSpace)$mesh$get_nodes())[2] + M = dim(extract_private(private$FunctionSpace_)$mesh_$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(extract_private(private$FunctionSpace)$mesh$get_times()) == 0){ - evals <- apply(private$FunctionSpace$get_basis()$eval(as.matrix(X)) %*% self$coeff, + if(length(extract_private(private$FunctionSpace_)$mesh_$time_nodes()) == 0){ + evals <- apply(private$FunctionSpace_$basis()$eval(as.matrix(X)) %*% private$coeff_, MARGIN=1, FUN=sum) }else{ - evals <- matrix(nrow=nrow(X),ncol=length(extract_private(private$FunctionSpace)$mesh$get_times())) - for(t in 1:length(extract_private(private$FunctionSpace)$mesh$get_times())){ - evals[,t] <- apply(private$FunctionSpace$get_basis()$eval(as.matrix(X)) %*% self$coeff[,t], + evals <- matrix(nrow=nrow(X),ncol=length(extract_private(private$FunctionSpace_)$mesh_$time_nodes())) + for(t in 1:length(extract_private(private$FunctionSpace_)$mesh_$time_nodes())){ + evals[,t] <- apply(private$FunctionSpace_$basis()$eval(as.matrix(X)) %*% private$coeff_[,t], MARGIN=1, FUN=sum) } } return(evals) }, set_coefficients = function(coefficients){ - if( nrow(coefficients) != private$FunctionSpace$get_basis$size()) - stop("Input parameter dimensions are different from the size of the function space.") - self$coeff <- coefficients + if( nrow(coefficients) != private$FunctionSpace_$basis()$size()) + stop("Input parameter dimension is different from function space size.") + private$coeff_ <- coefficients + }, + coefficients = function(){ + private$coeff_ + }, + FunctionSpace = function(){ + private$FunctionSpace_ } ) ) @@ -195,17 +207,17 @@ Function <- function(FunctionSpace) { #' @importFrom graphics plot #' @export plot.Function <- function(x, ...){ - mesh <- extract_private(extract_private(x)$FunctionSpace)$mesh - times <- mesh$get_times() + mesh <- x$FunctionSpace()$mesh() + times <- mesh$time_nodes() is_parabolic <- FALSE if(length(times)!=0) is_parabolic <- TRUE if(!is_parabolic){ - plot_data <- data.frame(X=mesh$get_nodes()[,1], - Y=mesh$get_nodes()[,2], - coeff=x$coeff[1:nrow(mesh$get_nodes())]) - I=mesh$get_elements()[,1]-1 - J=mesh$get_elements()[,2]-1 - K=mesh$get_elements()[,3]-1 + plot_data <- data.frame(X=mesh$nodes()[,1], + Y=mesh$nodes()[,2], + coeff=x$coefficients()[1:nrow(mesh$nodes())]) + I=mesh$elements()[,1]-1 + J=mesh$elements()[,2]-1 + K=mesh$elements()[,3]-1 fig<- plot_ly(plot_data, x=~X, y=~Y, z=~coeff, i = I, j = J, k = K, intensity=~coeff, color=~coeff, type="mesh3d", @@ -222,24 +234,24 @@ plot.Function <- function(x, ...){ eye = list(x = 1.25, y = -1.25, z = 1.25))) %>% colorbar(len = 1, title="") }else{ - plot_data <- data.frame(X=rep(mesh$get_nodes()[,1], times=length(times)), - Y=rep(mesh$get_nodes()[,2], times=length(times)), + plot_data <- data.frame(X=rep(mesh$nodes()[,1], times=length(times)), + Y=rep(mesh$nodes()[,2], times=length(times)), Z=rep(0, times=length(times)), - coeff=as.vector(x$coeff[1:nrow(mesh$get_nodes()),]), - times = round(rep(times, each=nrow(mesh$get_nodes())),3)) - coeff <- matrix(nrow=nrow(mesh$get_elements()), + coeff=as.vector(x$coefficients()[1:nrow(mesh$nodes()),]), + times = round(rep(times, each=nrow(mesh$nodes())),3)) + coeff <- matrix(nrow=nrow(mesh$elements()), ncol=length(times)) for(t in 1:length(times)){ - coeff[,t] <- apply(mesh$get_elements(), MARGIN=1, FUN= + coeff[,t] <- apply(mesh$elements(), MARGIN=1, FUN= function(edge){ - mean(x$coeff[edge,t]) + mean(x$coefficients()[edge,t]) }) } limits = c(min(coeff), max(coeff)) cmin = limits[1]; cmax=limits[2] - I=mesh$get_elements()[,1]-1 - J=mesh$get_elements()[,2]-1 - K=mesh$get_elements()[,3]-1 + I=mesh$elements()[,1]-1 + J=mesh$elements()[,2]-1 + K=mesh$elements()[,3]-1 fig<- plot_ly(plot_data, x=~X, y=~Y, z=~Z, frame=~times, i = I, j = J, k = K, cmin = limits[1], cmax=limits[2], intensity=~coeff, color=~coeff, type="mesh3d", @@ -270,14 +282,14 @@ plot.Function <- function(x, ...){ #' @return A plotly object #' @export setMethod("contour", signature=c(x="Function"), function(x, ...){ - mesh <- extract_private(extract_private(x)$FunctionSpace)$mesh - times <- mesh$get_times() + mesh <- x$FunctionSpace()$mesh() + times <- mesh$time_nodes() is_parabolic <- FALSE if(length(times)!=0) is_parabolic <- TRUE if(!is_parabolic){ - xrange <- range(mesh$get_nodes()[,1]) - yrange <- range(mesh$get_nodes()[,2]) + xrange <- range(mesh$nodes()[,1]) + yrange <- range(mesh$nodes()[,2]) Nx <- 40 Ny <- 40 eval_x <- seq(from=xrange[1], to=xrange[2], length.out=Nx) @@ -292,8 +304,8 @@ setMethod("contour", signature=c(x="Function"), function(x, ...){ yaxis = list(title = "", showgrid=F, zeroline=F, ticks="", showticklabels=F)) }else{ - xrange <- range(mesh$get_nodes()[,1]) - yrange <- range(mesh$get_nodes()[,2]) + xrange <- range(mesh$nodes()[,1]) + yrange <- range(mesh$nodes()[,2]) Nx <- 40 Ny <- 40 eval_x <- seq(from=xrange[1], to=xrange[2], length.out=Nx) diff --git a/R/mesh.R b/R/mesh.R index afee0bc..5f2c4d8 100644 --- a/R/mesh.R +++ b/R/mesh.R @@ -1,40 +1,37 @@ - -## Mesh Class Definition -#' @export .MeshCtr <- R6Class("Mesh", inherit = .DomainCtr, private = list( - cpp_handler = "ANY", ## cpp backend - m = 0L, ## local dim - n = 0L, ## embedding dim - times = vector(mode="double", length = 0L), ## time mesh for parabolic problems - time_step = NA ## time step for parabolic problems + cpp_handler_ = "ANY", ## cpp backend + m_ = 0L, ## local dim + n_ = 0L, ## embedding dim + time_nodes_ = vector(mode="double", length = 0L), ## time mesh for parabolic problems + time_step_ = NA ## time step for parabolic problems ), public = list( initialize = function(cpp_handler, m, n){ - private$cpp_handler <- cpp_handler - private$m <- m - private$n <- n + private$cpp_handler_ <- cpp_handler + private$m_ <- m + private$n_ <- n }, - get_nodes = function(){ - private$cpp_handler$nodes() + nodes = function(){ + private$cpp_handler_$nodes() }, - get_elements = function(){ - private$cpp_handler$elements()+1 + elements = function(){ + private$cpp_handler_$elements()+1 }, - get_boundary = function(){ - private$cpp_handler$boundary() + boundary = function(){ + private$cpp_handler_$boundary() }, - get_neighbors = function(){ - private$cpp_handler$neighbors() + neighbors = function(){ + private$cpp_handler_$neighbors() }, - get_times = function(){ - private$times + time_nodes = function(){ + private$time_nodes_ }, set_time_step = function(time_step){ - if(length(private$time_interval)==0) + if(length(private$time_interval_)==0) stop("Error!") - private$time_step <- time_step - private$times <- seq(private$time_interval[1], private$time_interval[2], by=time_step) + private$time_step_ <- time_step + private$time_nodes_ <- seq(private$time_interval_[1], private$time_interval_[2], by=time_step) } ) ) @@ -49,7 +46,7 @@ #' an entry '0' indicates that the corresponding node is not a boundary node.} #' } #' -#' @return An S4 object representing a Mesh. +#' @return A R6 class representing a Mesh. #' @rdname Mesh #' @export #' @examples @@ -104,7 +101,7 @@ setMethod("Mesh", signature=c(domain="triangulation"), # # @param op1 A mesh object created by \code{Mesh}. # @param op2 A numeric vector. -# @return An S4 object representing a spatio-temporal domain. +# @return A R6 object representing a spatio-temporal domain. # @rdname Domain_times_vector # @export #setGeneric("%X%", function(op1, op2) standardGeneric("%X%")) @@ -114,19 +111,19 @@ setMethod("%X%", signature=c(op1="Mesh", op2="numeric"), function(op1, op2){ if(op2[1] > op2[length(op2)]) stop("Error! First time instant is greater than last time instant.") - set_times(op1, op2) - set_time_interval(op1, c(op2[1], op2[length(op2)])) - set_time_step(op1, (op2[2] - op2[1])) + set_private(op1, "time_nodes_", op2) + set_private(op1, "time_interval_", c(op2[1], op2[length(op2)])) + set_private(op1, "time_step_", (op2[2] - op2[1])) op1 }) ## Mesh - auxiliary methods # unroll_edges_aux <- function(mesh){ -# edges <- matrix(nrow=3*nrow(mesh$get_elements()), ncol=2) -# for(i in 1:nrow(mesh$get_elements())){ -# edges[(3*(i-1) + 1),] = mesh$get_elements()[i,c(1,2)] -# edges[(3*(i-1) + 2),] = mesh$get_elements()[i,c(2,3)] -# edges[(3*(i-1) + 3),] = mesh$get_elements()[i,c(3,1)] +# edges <- matrix(nrow=3*nrow(mesh$elements()), ncol=2) +# for(i in 1:nrow(mesh$elements())){ +# edges[(3*(i-1) + 1),] = mesh$elements()[i,c(1,2)] +# edges[(3*(i-1) + 2),] = mesh$elements()[i,c(2,3)] +# edges[(3*(i-1) + 3),] = mesh$elements()[i,c(3,1)] # } # edges # } @@ -134,11 +131,11 @@ setMethod("%X%", signature=c(op1="Mesh", op2="numeric"), setGeneric("unroll_edges", function(mesh) standardGeneric("unroll_edges")) setMethod("unroll_edges", "Mesh", function(mesh){ #unroll_edges_aux(mesh) - edges <- matrix(nrow=3*nrow(mesh$get_elements()), ncol=2) - for(i in 1:nrow(mesh$get_elements())){ - edges[(3*(i-1) + 1),] = mesh$get_elements()[i,c(1,2)] - edges[(3*(i-1) + 2),] = mesh$get_elements()[i,c(2,3)] - edges[(3*(i-1) + 3),] = mesh$get_elements()[i,c(3,1)] + edges <- matrix(nrow=3*nrow(mesh$elements()), ncol=2) + for(i in 1:nrow(mesh$elements())){ + edges[(3*(i-1) + 1),] = mesh$elements()[i,c(1,2)] + edges[(3*(i-1) + 2),] = mesh$elements()[i,c(2,3)] + edges[(3*(i-1) + 3),] = mesh$elements()[i,c(3,1)] } edges }) @@ -146,18 +143,18 @@ setMethod("unroll_edges", "Mesh", function(mesh){ plot_mesh_aux <- function(x, ...){ edges <- unroll_edges(x) plot_ly(...) %>% - add_markers(x = x$get_nodes()[,1], - y = x$get_nodes()[,2], + add_markers(x = x$nodes()[,1], + y = x$nodes()[,2], color = I('black'), size = I(1), hoverinfo = 'text', - text = paste('
Coordinates:', round(x$get_nodes()[,1],2), - round(x$get_nodes()[,2],2)), + text = paste('
Coordinates:', round(x$nodes()[,1],2), + round(x$nodes()[,2],2)), showlegend = T, visible = T) %>% - add_segments(x = x$get_nodes()[edges[,1],1], - y = x$get_nodes()[edges[,1],2], - xend = x$get_nodes()[edges[,2],1], - yend = x$get_nodes()[edges[,2],2], + add_segments(x = x$nodes()[edges[,1],1], + y = x$nodes()[edges[,1],2], + xend = x$nodes()[edges[,2],1], + yend = x$nodes()[edges[,2],2], color = I('black'), size = I(1), showlegend = F) %>% layout( diff --git a/R/operators.R b/R/operators.R index 82059bf..522a043 100644 --- a/R/operators.R +++ b/R/operators.R @@ -1,11 +1,16 @@ # gradient of Function .FunctionGradCtr <- R6Class("FunctionGrad", + private = list( + Function_ = "Function" + ), public = list( - f = "Function", ## Function of which the gradient is taken 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 + initialize = function(Function, K=matrix(c(1,0,0,1), nrow=2,ncol=2,byrow=T)){ + private$Function_ <- Function self$K <- K + }, + Function = function(){ + private$Function_ } ) ) @@ -21,7 +26,7 @@ setOldClass(c("FunctionGrad", "R6")) #' compute gradient of Function #' -#' @param f a Function created by \code{Function}: +#' @param Function a Function created by \code{Function}: #' @return An S4 object representing the gradient of the Function passed as parameter. #' @export #' @rdname grad @@ -34,11 +39,11 @@ setOldClass(c("FunctionGrad", "R6")) #' f <- Function(Vh) #' grad_f <- grad(f) #' } -setGeneric("grad", function(f) standardGeneric("grad")) +setGeneric("grad", function(Function) standardGeneric("grad")) #' @rdname grad -setMethod("grad", signature(f = "Function"), function(f) { - .FunctionGradCtr$new(f = f) +setMethod("grad", signature(Function = "Function"), function(Function) { + .FunctionGradCtr$new(Function = Function) }) # setMethod("*", signature=c(e1="matrix", e2="FunctionGrad") @@ -51,22 +56,27 @@ setMethod("grad", signature(f = "Function"), function(f) { #' @export `*.FunctionGrad` <- function(e1, e2){ if (!is.matrix(e1)) stop("First argument must be a matrix!") - .FunctionGradCtr$new(f = e2$f, K = e1) + .FunctionGradCtr$new(Function = e2$Function(), K = e1) } ## base class for differential operators .DiffOpCtr <- R6Class("DiffOpObject", + private = list( + ## Function to which the operator is applied + Function_ = "Function" + ), public = list( - initialize = function(tokens, params, f){ + initialize = function(tokens, params, Function){ self$tokens <- tokens self$params <- params - self$f <- f + private$Function_ <- Function }, ## single blocks composing the overall operator tokens = "vector", params = "list", - ## Function to which the operator is applied - f = "Function" + Function = function(){ + private$Function_ + } ) ) @@ -89,7 +99,7 @@ setOldClass(c("DiffOpObject", "R6")) #' @name DifferentialOps #' @export `+.DiffOpObject` <- function(e1, e2){ - if (!identical(e1$f, e2$f)) { + if (!identical(e1$Function(), e2$Function())) { stop("operator arguments must be the same") } ## check for duplicated operator tokens @@ -99,7 +109,7 @@ setOldClass(c("DiffOpObject", "R6")) .DiffOpCtr$new( tokens = c(e1$tokens, e2$tokens), params = c(e1$params, e2$params), - f = e1$f + Function = e1$Function() ) } @@ -117,7 +127,7 @@ setOldClass(c("DiffOpObject", "R6")) if(missing(e2)){ e1$params[[1]] <- -e1$params[[1]] return(e1) - }else if (!identical(e1$f, e2$f)) { + }else if (!identical(e1$Function(), e2$Function())) { stop("operator arguments must be the same") } ## check for duplicated operator tokens @@ -130,7 +140,7 @@ setOldClass(c("DiffOpObject", "R6")) .DiffOpCtr$new( tokens = c(e1$tokens, e2$tokens), params = c(e1$params, e2$params), - f = e1$f + Function = e1$Function() ) } @@ -161,7 +171,7 @@ setOldClass(c("DiffOpObject", "R6")) if (!is.numeric(e1)) stop("bad product") e2$params[[1]] <- e1*e2$params[[1]] e2 - } +} ## diffusion term .DiffusionCtr <- R6Class("DiffusionOperator", @@ -180,7 +190,7 @@ setOldClass(c("DiffusionOperator", "DiffOpObject")) #' laplace operator for Function class #' -#' @param f a Function. +#' @param Function a Function. #' @return A S4 object representing the laplace operator applied to the function passed as parameter. #' @rdname DiffusionOperator #' @export @@ -193,14 +203,14 @@ setOldClass(c("DiffusionOperator", "DiffOpObject")) #' f <- Function(Vh) #' laplace_f <- laplace(f) #' } -laplace <- function(f) { - if (!is(f, "Function")) { +laplace <- function(Function) { + if (!is(Function, "Function")) { stop("wrong argument type") } .DiffusionCtr$new( tokens = "diffusion", params = list(diffusion = 1), - f = f + Function = Function ) } @@ -208,7 +218,7 @@ laplace <- function(f) { #' divergence operator FunctionGrad #' -#' @param f a Function. +#' @param Function a Function. #' @return A S4 object representing the diffusion term of a second order linear differential operator. #' @rdname DiffusionOperator #' @export @@ -222,13 +232,13 @@ laplace <- function(f) { #' K <- matrix(c(1,2,1,0),nrow=2,ncol=2) #' diffusion <- div(K*grad(f)) #' } -div <- function(f) { - if (is(f, "FunctionGrad")) { - if (!is.null(f$K)) { +div <- function(Function) { + if (is(Function, "FunctionGrad")) { + if (!is.null(Function$K)) { return(.DiffusionCtr$new( tokens = "diffusion", - params = list(diffusion = f$K), - f = f$f) + params = list(diffusion = Function$K), + Function = Function$Function()) ) } } @@ -272,7 +282,7 @@ setMethod("dot", signature(op1 = "vector", op2 = "FunctionGrad"), .TransportCtr$new( tokens = "transport", params = list(transport = as.matrix(op1)), - f = op2$f + Function = op2$Function() ) }) @@ -311,7 +321,7 @@ setOldClass(c("ReactionOperator", "DiffOpObject")) .ReactionCtr$new( tokens = "reaction", params = list(reaction = e1), - f = e2 + Function = e2 ) } @@ -346,8 +356,8 @@ setOldClass(c("TimeDerivative", "DiffOpObject")) setMethod("dt", signature = c(x="Function", df="missing", ncp="missing"), function(x,df,ncp){ .TimeDerivativeCtr$new( - tokens="time", + tokens = "time", params = list(time=1L), - f=x + Function = x ) }) diff --git a/R/pde.R b/R/pde.R index cdc1d51..c6e1144 100644 --- a/R/pde.R +++ b/R/pde.R @@ -7,13 +7,13 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3) is_dirichletBC_set = "logical", is_initialCondition_set = "logical", is_parabolic ="logical", - cpp_handler = "ANY", - DifferentialOp = "DiffOpObject" + cpp_handler_ = "ANY", + DifferentialOp_ = "DiffOpObject" ), public = list( initialize = function(cpp_handler, DifferentialOp, is_parabolic, is_dirichletBC_set, is_initialCondition_set){ - private$cpp_handler <- cpp_handler - private$DifferentialOp <- DifferentialOp + private$cpp_handler_ <- cpp_handler + private$DifferentialOp_ <- DifferentialOp private$is_parabolic <- is_parabolic private$is_dirichletBC_set <- is_dirichletBC_set private$is_initialCondition_set <- is_initialCondition_set @@ -21,26 +21,24 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3) solve = function(){ if(!private$is_dirichletBC_set){ if(!private$is_parabolic){ - dirichletBC_ = as.matrix(rep(0,times=nrow(private$cpp_handler$get_dofs_coordinates()))) + dirichletBC_ = as.matrix(rep(0,times=nrow(private$cpp_handler_$dofs_coordinates()))) }else{ - times <- extract_private( - extract_private( - private$DifferentialOp$f)$FunctionSpace)$mesh$get_times() - dirichletBC_ = matrix(0, nrow=nrow(private$cpp_handler$get_dofs_coordinates()), + times <- private$DifferentialOp_$Function()$FunctionSpace()$mesh()$time_nodes() + dirichletBC_ = matrix(0, nrow=nrow(private$cpp_handler_$dofs_coordinates()), ncol=length(times)) } - private$cpp_handler$set_dirichlet_bc(dirichletBC_) + private$cpp_handler_$set_dirichlet_bc(dirichletBC_) private$is_dirichletBC_set <- TRUE } if(private$is_parabolic & (!private$is_initialCondition_set)){ stop("initialCondition must be provided.") } - private$cpp_handler$solve() - private$DifferentialOp$f$coeff <- private$cpp_handler$solution() + private$cpp_handler_$solve() + private$DifferentialOp_$Function()$set_coefficients(private$cpp_handler_$solution()) invisible(self) }, - get_dofs_coordinates = function(){ - private$cpp_handler$get_dofs_coordinates() + dofs_coordinates = function(){ + private$cpp_handler_$dofs_coordinates() }, set_boundary_condition = function(fun, type="dirichlet", on=NULL){ if(!any(type == c("dirichlet", "Dirichlet", "d"))) @@ -48,30 +46,26 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3) dirichletBC_ <- NULL if(any(typeof(fun) == c("function", "closure"))){ if(!private$is_parabolic){ - dirichletBC_ <- as.matrix(fun(private$cpp_handler$get_dofs_coordinates())) + dirichletBC_ <- as.matrix(fun(private$cpp_handler_$dofs_coordinates())) }else{ - times <- extract_private( - extract_private( - private$DifferentialOp$f)$FunctionSpace)$mesh$get_times() - dirichletBC_ <- fun(private$cpp_handler$get_dofs_coordinates(), times) + times <- private$DifferentialOp_$Function()$FunctionSpace()$mesh()$time_nodes() + dirichletBC_ <- fun(private$cpp_handler_$dofs_coordinates(), times) } }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){ - if(nrow(as.matrix(fun)) == nrow(private$cpp_handler$get_dofs_coordinates())){ + if(nrow(as.matrix(fun)) == nrow(private$cpp_handler_$dofs_coordinates())){ dirichletBC_ <- fun }else if (nrow(as.matrix(fun)) == 1L){ if(!private$is_parabolic) - dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler$get_dofs_coordinates()), ncol=1) + dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler_$dofs_coordinates()), ncol=1) else{ - times <- extract_private( - extract_private( - private$DifferentialOp$f)$FunctionSpace)$mesh$get_times() - dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler$get_dofs_coordinates()), + times <- private$DifferentialOp_$Function()$FunctionSpace()$mesh()$time_nodes() + dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler_$dofs_coordinates()), ncol=length(times)) } } } private$is_dirichletBC_set <- TRUE - private$cpp_handler$set_dirichlet_bc(dirichletBC_) + private$cpp_handler_$set_dirichlet_bc(dirichletBC_) invisible(self) }, set_initial_condition = function(fun){ @@ -79,22 +73,22 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3) stop("Cannot set initial condition for elliptic problem.") private$is_initialCondition_set <- TRUE if(typeof(fun) == "closure"){ - private$cpp_handler$set_initial_condition(fun(private$cpp_handler$get_dofs_coordinates())) + private$cpp_handler_$set_initial_condition(fun(private$cpp_handler_$dofs_coordinates())) }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){ - if(nrow(as.matrix(fun)) == nrow(private$cpp_handler$get_dofs_coordinates())){ - private$cpp_handler$set_initial_condition(fun) + if(nrow(as.matrix(fun)) == nrow(private$cpp_handler_$dofs_coordinates())){ + private$cpp_handler_$set_initial_condition(fun) }else if(nrow(as.matrix(fun)) == 1L){ - private$cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(private$cpp_handler$get_dofs_coordinates()), + private$cpp_handler_$set_initial_condition(matrix(fun, nrow=nrow(private$cpp_handler_$dofs_coordinates()), ncol=1)) } } invisible(self) }, - get_mass = function(){ - private$cpp_handler$mass() + mass = function(){ + private$cpp_handler_$mass() }, - get_stiff = function(){ - private$cpp_handler$stiff() + stiff = function(){ + private$cpp_handler_$stiff() } ) ) @@ -129,8 +123,7 @@ parse_pde_parameters <- function(DifferentialOp){ } if(pde_type == pde_type_list$parabolic) - pde_parameters[["time_mesh"]] <- as.vector( - extract_private(extract_private(DifferentialOp$f)$FunctionSpace)$mesh$get_times()) + pde_parameters[["time_nodes"]] <- as.vector(DifferentialOp$Function()$FunctionSpace()$mesh()$time_nodes()) 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) @@ -146,13 +139,8 @@ make_pde <- function(DifferentialOp) { pde_parameters <- parse_pde_parameters(DifferentialOp) ## define Rcpp module - D <- extract_private( - extract_private( - extract_private( - DifferentialOp$f)$FunctionSpace)$mesh)$cpp_handler ## domain - fe_order <- extract_private( - extract_private( - DifferentialOp$f)$FunctionSpace)$fe_order ## finite element order + D <- extract_private(DifferentialOp$Function()$FunctionSpace()$mesh())$cpp_handler_ + fe_order <- DifferentialOp$Function()$FunctionSpace()$fe_order() cpp_handler <- NULL if (fe_order == 1) { ## linear finite elements cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters) @@ -170,7 +158,7 @@ make_pde <- function(DifferentialOp) { #' @return A R6 class representing a PDE. #' @rdname pde #' @export -setGeneric("Pde", function(DifferentialOp,f) standardGeneric("Pde")) +setGeneric("Pde", function(DifferentialOp, f) standardGeneric("Pde")) #' @rdname pde setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"), @@ -179,11 +167,10 @@ setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"), pde_type = extract_pde_type(DifferentialOp) cpp_handler <- make_pde(DifferentialOp) - quad_nodes <- cpp_handler$get_quadrature_nodes() + quad_nodes <- cpp_handler$quadrature_nodes() ## evaluate forcing term on quadrature nodes if(pde_type == pde_type_list$parabolic) { - times <- extract_private( - extract_private(DifferentialOp$f)$FunctionSpace)$mesh$get_times() + times <- DifferentialOp$Function()$FunctionSpace()$mesh()$time_nodes() cpp_handler$set_forcing(as.matrix(f(quad_nodes, times))) }else{ cpp_handler$set_forcing(as.matrix(f(quad_nodes))) diff --git a/R/sf.R b/R/sf.R index 72a2d39..24db96a 100644 --- a/R/sf.R +++ b/R/sf.R @@ -15,23 +15,23 @@ #' @export st_as_sf.Domain <- function(x, ...){ geom_sf <- list() - for(sub_id in 1:length(extract_private(x)$geometry)){ + for(sub_id in 1:length(extract_private(x)$geometry_)){ - geometry <- extract_private(x)$geometry[[sub_id]] + geometry <- extract_private(x)$geometry_[[sub_id]] path_id <- unique(geometry$edges_ring) path_list = vector(mode="list", length=length(path_id)) for(i in 1:length(path_id)){ path <- t(geometry$edges[geometry$edges_ring==path_id[i],][,1]) path <- c(path,path[1]) - nodes <- as.matrix(extract_private(x)$coords[path,1:2]) + nodes <- as.matrix(extract_private(x)$coords_[path,1:2]) path_list[[i]] <- st_linestring(nodes) } # edges to st_linestring ! (FAST) edge_list <- lapply(as.list(as.data.frame(t(geometry$edges))), FUN= function(edge){ - st_linestring(as.matrix(extract_private(x)$coords[edge,1:2])) + st_linestring(as.matrix(extract_private(x)$coords_[edge,1:2])) }) # nodes to st_point ! (FAST) @@ -46,7 +46,7 @@ st_as_sf.Domain <- function(x, ...){ 1:nrow(geometry$nodes)) group <- as.factor(group) - crs <- extract_private(x)$crs + crs <- extract_private(x)$crs_ if(is.na(crs)) crs <- NA_crs_ geom_sfc <- st_sfc( append(append(path_list, edge_list), pts_list), crs=crs) geom_sf[[sub_id]] <- st_as_sf(data.frame(label= label, @@ -54,7 +54,7 @@ st_as_sf.Domain <- function(x, ...){ geometry = geom_sfc) } - if(length(extract_private(x)$geometry)==1) geom_sf <- geom_sf[[1]] + if(length(extract_private(x)$geometry_)==1) geom_sf <- geom_sf[[1]] return(geom_sf) } @@ -67,25 +67,25 @@ st_as_sf.Domain <- function(x, ...){ #' @rdname sf #' @export st_as_sfc.Domain <- function(x, ...){ - polygon_list <- list(mode="list", length=length(extract_private(x)$geometry)) - for(sub_id in 1:length(extract_private(x)$geometry)){ - geometry <- extract_private(x)$geometry[[sub_id]] + polygon_list <- list(mode="list", length=length(extract_private(x)$geometry_)) + for(sub_id in 1:length(extract_private(x)$geometry_)){ + geometry <- extract_private(x)$geometry_[[sub_id]] path_id <- unique(geometry$edges_ring) path_list = vector(mode="list", length=length(path_id)) for(i in 1:length(path_id)){ path <- t(geometry$edges[geometry$edges_ring==path_id[i],])[1,] path <- c(path,path[1]) - nodes <- as.matrix(extract_private(x)$coords[path,1:2]) + nodes <- as.matrix(extract_private(x)$coords_[path,1:2]) path_list[[i]] <- nodes } polygon_list[[sub_id]] <- st_cast(st_polygon(path_list), to="MULTIPOLYGON") } - crs <- extract_private(x)$crs + crs <- extract_private(x)$crs_ if(is.na(crs)) crs <- NA_crs_ - if(length(extract_private(x)$geometry)==1){ + if(length(extract_private(x)$geometry_)==1){ result <- st_sfc(polygon_list[[1]], crs=crs) }else{ result <- st_sfc(polygon_list, crs=crs) @@ -101,13 +101,13 @@ st_as_sfc.Domain <- function(x, ...){ #' @export st_as_sfc.Mesh <- function(x, ...){ - polygon_list <- apply(x$get_elements(), MARGIN=1, FUN=function(elem){ - st_cast(st_linestring(x$get_nodes()[elem,]), + polygon_list <- apply(x$elements(), MARGIN=1, FUN=function(elem){ + st_cast(st_linestring(x$nodes()[elem,]), to ="POLYGON" ) }) - crs <- extract_private(x)$crs + crs <- extract_private(x)$crs_ if(is.na(crs)) crs <- NA_crs_ mesh_sf <- st_sfc(polygon_list, crs = crs) mesh_sf @@ -126,7 +126,6 @@ st_crs.Domain <- function(x, ...){ `st_crs<-.Domain` = function(x, value) { sfc <- st_as_sfc(x) st_crs(sfc) <- value - set_crs(x, value) - #crs <- st_crs(sfc)$input + set_private(x, "crs_", value) } diff --git a/R/utils.R b/R/utils.R index 5196860..9a2968f 100644 --- a/R/utils.R +++ b/R/utils.R @@ -5,34 +5,10 @@ setMethod("extract_private", signature = "R6", x$.__enclos_env__$private }) -setGeneric("set_geometry", function(x, value) standardGeneric("set_geometry")) -setMethod("set_geometry", signature = "Domain", - function(x, value){ - invisible(x$.__enclos_env__$private$geometry <- value) -}) - -setGeneric("set_time_interval", function(x, value) standardGeneric("set_time_interval")) -setMethod("set_time_interval", signature = "Domain", - function(x, value){ - invisible(x$.__enclos_env__$private$time_interval <- value) -}) - -setGeneric("set_crs", function(x, value) standardGeneric("set_crs")) -setMethod("set_crs", signature = "Domain", - function(x, value){ - invisible(x$.__enclos_env__$private$crs <- value) -}) - -setGeneric("set_times", function(x, value) standardGeneric("set_times")) -setMethod("set_times", signature = "Mesh", - function(x, value){ - invisible(x$.__enclos_env__$private$times <- value) -}) - -setGeneric("set_time_step", function(x, value) standardGeneric("set_time_step")) -setMethod("set_time_step", signature = "Mesh", - function(x, value){ - invisible(x$.__enclos_env__$private$time_step <- value) +setGeneric("set_private", function(x, attribute, value) standardGeneric("set_private")) +setMethod("set_private", signature = c("R6", "character"), + function(x, attribute, value){ + invisible(x$.__enclos_env__$private[[attribute]] <- value) }) # extract_private <- function(x){ diff --git a/data/lake_como_boundary.RData b/data/lake_como_boundary.RData new file mode 100644 index 0000000000000000000000000000000000000000..beafe0f9307ea6c297023b92da69cad836b129ff GIT binary patch literal 11456 zcmV;xEI-p9iwFP!000001D%@(IM#pr_wAWgh$vi&lx$gzuT6F-iqd5!Gb5`+QOKsq zs;ojnc99Vp8rmYGMRrEAd*0W5|GMx0{d@lZ=XjpSan#FqT%U1{_j!Jly;_#?ES4-( zR8+K7bksCdwDkCvihjSj&NgK#DmofmP*c%SF@vjJ$u15yb{=jXHntwV?)GG_0DNi< zp3vbB=YM%~l_A;5L)hL~SQ9V*?b-HYK0cn_s?yT-4$f}gQXXDT(u&G5^4sl|WVb2H zD#&hAkg=2BMwW4u+omXEt88m8Yinn3Z--B^Q$`_mM^y|C=;;&LB*X+73T> zw`M0uwZNE31#YBnI4@6-9&U_upxKY>s;JBpcL`Ffqq4!~3WC%!XP0yjb1<@S5JAb*G;>o`Gb8aaMIz!}UIdRa}m zj^Zx9u?`gbS!q2&kQ#?)R1ZfHqz3O?@=kh!bYGYClI}BtRNvh8Xx9Ufp;Fh;euO#-KAXU7ne&~qtD$H)}Ka6$WNjc5J zT}P10NJ*FHUJ|5IJMX8$76hrJLZ$d{0zta{t|oQDfglw>+xmEYGRSSeYsY7Tbc^pp z)muJ-RAgwb>@rD^3a7JVd$Fz}y7|y9Ib4(81@aFNqym*IQ3}&w`)Kbu*t4Kw@O0G# zsQE(Yc_l%*xxGGU8T+_#LxyxSn;_+XiaGtI1I)^JT02aT@_vgqrnG_yWjF8j5~LjR z85xtq1Sva2DOzBWAYG5h_B6^RNY?@xBR{NwY*HL+@ZL4KPbW8G{HsfCg8iRB9qS@S zLD0UoM<|9MWo6TH9)drz9$4&qTaEi!#Tj2UKt7x4f{Wmz4~{8M!4<|mCmF%E-j95+ z$JGz-9@%#jq|7|?xkKRuDdR00o838rl>U>G_ZjS-uGFC;8A*`RP6V1tydX$d+AD^h z^n-8GLXPgj{XK7UHJ9Q;a1ro&HZ8#Ng+<|%C8P%+22 zTrf5TXFqc`%5$8 z?{e?lyMS+C$Vrvk2MN-p62FNOcF@Iu@0Tp-_{{RKG(k#!n4y)5eI$3Yu0PL(`-%&{ z55RwycEvFt(+53ycaUm8z1Zb&f4s*pN_U@3kP-*-T>S*VWOc2jt01G`2GcTv6tDDM z`4a?{ z2J`(&o=m{{Gv&z7-}B;H`<=W)BPEW^E?>KY^;@q~Xf(vSwI4e~!%wM6L4uE7gP*T< zzNp807q3=HBK|X21nWW&A6b>RBg2n?Dev0ysz9#TdRFL1))DHbI`CVT@|;X^JJz#p z^2~|n$Rl3M=z&|vGk*Jf4TziT7FlyMAHb3)e1=)zsrM54uva#F-ANj(Gkc&@!WX)b zV>LJYfWlzK19@j~|IY1IA;{xw@r8?z!RyqjaT}07*ScwEo8mw*ccKJ&b#;51Y{e6- zqrI1&Ov%rls-mSr+^;^;dQu7RW%{S|Hc;k!?l0O0{`OC}iu}plq8XXITn)Y5k!KIrz;A^MwhvA~ z*9)D)kpJ|rG^f@nqoN(yX6J-$bAsSeY)U|tU&PU`^Qrl|E5#_ zX#v=|z*1UO_%X)W5qPQ)_9+^CXmq>@x|wvH`UB#yc;~Id40ylfjBf;&2XuE|RACeJ ztn~NoRqtWf((*LTQ24*Jx|rnzbfnB8d%y$rrF4;vsx{v8F_t=|5~^z z>ejunXG!DG+d8)t$5-%YefbXe^N7>?9!Is(QBUq4X?n96>$xBKPu5Phc>k@C{mVd z&jMOMKQj6o^_Vra<{I*@wc_0nGuG3Vx7~#c_3V+%+Kgz_`zO3Zt3D&{9{257P}+<8 zZg0Y`%-}sX6Jc}grzKu5VG?oI>_93OK|N_Yck|mb)Qv`u^Z9qbQr5lveFpxik1qbC zjB)B6q%PE4f!!XvRrI5--2J{okjfr>WT4YxQQld)WJLHT2e#MGb#p`$>&cHP|*FMYZ5OxMkEk47#3HWO=mBj#5|IX=nFA@1hlV>`RC4Puc18HBjo4;|=)$)CF6) z*XoGp1emWqJ<4UJody(&P7F!vr5mzxo6<520z%RL5 zwNZ~^-O7@xVc%F*o%hwK$FYYMy1L<~*llY?ve-ba8jq8hCzj4l{Q>GsEPsS|#APt2 zUoH1GxR8?-hCVVbd}9I`aTFKo_^Lbw_gAL%lCa*`7wJ!JkK#V(9pS?UnD=|F=s5Bt zW<0O^Tm+u)S~UBG{w;QgQW~Y+#8kyM2SMj!KI{%E z0KJZ5%aj^MJj8lb228yH@BCVzhTmfSQu;HHS1~Ihs_}@k=rF72XVCv$Bs)gVP^epU zLh%$qisIbtBnN$q^sDv>N1qWPO*gHFan93-F22P6!ewhC&3eG3f!}th`{CyFuXeyr zVI9kLYv8wY&BkK`Vg%{z_vu+7Gq5wEbp!Sl`b9UT0r?mDQjGlw>QQLKj*QMJI=bPk-x*YJ1 z-?;~UL!e$`*Ei@#PoTdUw}V zXL|iEn&yj`5F3{yGPp zJh|&vt{C#fC-dUR2I#eST~2~2*6(E&pwWZ+e4>-_?LO2E&w870+r0>qN6z=9{^JD6 zeY4oh67(HzYvU$^@Hv;I<5!uGzs}F9*X5w^a563Yrh@SuwYA&^q4$oWN&>De1j&JE zH+>i4!#-43TLAvFFEo$fNBwbReSY!{{O_2r#F~TgorFU(IWVu2;K0HJ;=w7PamOXp zO~;p+tnToGqd4`^ejR+?D|qCs5J9rLH}A)X{Iz4KV%>l|C+{5@-+}YZ(ctFF^QccY zww~U;sMAOE4t;Dw{X2YW(-kS?iS^S~m2&ixR<8D6o1u$`(&jujAwLerMev6pjx4=G zJeRS4OAn2~)wKl4BBt;yWxWR`PpZYE4?Ez&+i(y2-Tzv`zYTrEem)yB1@vL&ypz(D zK5Spabd((Gr&;Ov1~o^5WIQLj=@;_fNHi~w9pe~^p8M>9br^^aziowH=qIPmCnMi= zN@-Q7X$g|HPy%zh6+zMzn*F%2i6Cj{%}gypuQkk{3^>Am>T|1{T@erJ&vh#8p?~Uc zYqIT-$LgtC#`_U>>J8Uh3Yqb~6@vgc1+vPF%NmGBy)){%PSy@~A7W${< zJY`0Mc+n;&oz_IZuHCy&gBRx!of?-gW7u8Ka;kT31NOK6c}6VaaF0-D@L8M_^mDy7 zSV@7WUYJS1&iZ`A@^R>k^>S~WI}Ux-9oaq=_8sHAR9r`ob?7V~XSw1EvTFHW%0Jm(e)%4WivW4{1@4NrM$sgNPr}s5~zh@{~lb zju$~46&W%*!h-W3Y3%~5BJ`CcG_w?geu$KQzi5^TpEpWh2}IrBSW1p35RaQS)o30< zT@lmf4>=Wv^@&ek^Xzs;P=G}$=MZr*(S1}F4vbN#RTrRbx#9A)t4 zg*~?p?{R&P&uvpWy4x9XDs8jjs44PHMrZAm3hJLsuaM*zbXKO@$fXXtDtF1?!uuZt zNuEt>$q`hzQcACh{8H!|J2ZrO6mFKXQ{%j>==>ms3i+$pH|07C-BkP#!PgBtDXqy) zWk!8a7EtA>hYqPImhYR#=T)~RFV+9V^+0L^JKj^tX_ym2->1q)#m0+qRV!vzvQi0> z+LKFRr{J$04qc1I(D$A4ZJca4m+X8~7~G8Y?<%UxJ%e~B_Ez+Ds~|sJ%DKKFk2LZk z(%P_o%^c5fSKtTDk--m|h;wbL0n1aU4?6m%&Kkq-x(~0?IGZ9r8n1j|M4sx?iq&ge z1Y_pA^w5{>>0CoD+ygz);*~Cd?ip15u;#+JhMF;I8aO8y291~c!Y)RK7l%gBPZ{@a zjqpaFYtr=k%zC_UGPC|!kv3RAHD`^wXR_Vr>!&Z^TSemx^gSkZ7u-A1hntvI&|atm z&pF1wL0+1uid^(TJ{r$12#e#qYm(Zby9)Vjd@wV@SP?oay^@<1W-EetV4o8u#wnJ_YcJ`l^O?{cuR5>k3hd8YcsXGR{gclINddMtFv;X74dUN7x6slJ zdhN$XAIX7pvR`M6`Vuvz?@?@Jm!$M7!#w&p*9R5UvG-&Dr`;2+H(}m0##z(9ao!8r zb@pC6<_~F_3KT@$Ip;dOXbpRZ3%-n`Mg0oj&z;_id^+#3@P-|I$^|WL=)*ioMl8fw+$-*S8Ep9~H6FezhC=;Rt2hP^MF?VVW z>@P8jKlm1OAu({^vpwo&lKgKo4_qe=8eO~sezz(Ua*L11&QGRZ=U2zX|@Cxaq z@Ojj!D{q$@*1@i6CFj?9!Vl?h`eP0vE;Gzb-J+1UnYTWkDTeM};+0H4uLs z{o(a!ky=mqJDXgW_{bLLLZ8c#>FASkhNni-kWYEhMphiCK3CY3 z>^WTzF4jk{AReo?INODz4p;jv-swR;)i573Q$$}{6L%n25$CJg?z>&h@Y~%Vspq0# zmwRQC{G+IU_o}vrC8GY_bDD2_j{LrN>1+IBoPX-f4dOKLxw=BhLVNh_{*^UyPOxhu zwexKW*t=1@Nl*fHyK&h*rw4uD18zTuVvXcfFc=!~gce^S{F+I_a|KMMwM)?37a7iG^JGy$(Fyt3Z{KKk{T`6f7Q zq4Dt;p))ROxvO3d+&7(`Hvk4FUBNzd##`QB&SwX^?s}1#37rWhM-B!&pV%Nu^*sZ0 zZ91Ob3-*Lan{Wic#OlPGaBqPi9KOM=uFYRTB-5`cPiQ6`2ZeL_F4W4o<6H05(gGCv9`s4clwR; zPJ+3|Dm^2>b5Zk+E};9<#b;YVzwu6IMbP1`a@_=IOJz+&^WX29%UCCFEuSt~uci$P#;e)2k6aKN<3sgKb4wS7NHb2E6>sPz6b z@Co-+R0yFnw@Rw_kr}A5hvAtz__6ia;dkKgkLzhL-aKy&9b-Re(Z7A}Dd;TN)BF>R zy>Q3$95}hn&|{R)SSrF!dpLsmN)G3g+c$o?bNRrjKPK59c+ArE;TRhiyvVWI(%!d z1(mA31|EQJCZi3p;1?fL)q{jC4R`RjoGaiq4*jQNprTEGxHlN3{Q6ZK_(Hvao&cvL zG!1wNUD^$IH!HJ(nllM1Dxlqc-^EeTZH)n!9Qd&zcuzSvJnVfP5RRfInW<@daH zNdyB5rD-$4lKVo&PT*krxtuIQm!ALXk$xf2^m_h|7vM>;JJeHPwju|gFW50i&2)y) zW#Cgfuweic_`UY`F_1iHt5*bC?@STh2A)*ivEvYUi(hB+2C(UEll=&x%eZlaXgua) zRBtkm#XcAVlV@mGfi;XRvCM=n6XU*SCQ(qT*x6_p^zqXUodYk=ZELLuGakuLV4X~r z%5na%CsW@CChH+WmzhCnnkb@hzE% zst8?nb^+?ubf9?A$~D9jyE%{X2-eAdir&ik6d1OtqV*xz{yEr}0i50y%xFRAa>(qt zd>;1au;3Iqw;goaChkTLIwpt=ATBrp$X0rTV5K3OR61DpyX`wW*to%?6LH7Ub*s>G z1zdWvcM%mdrl!1UuC?N`AvzD171V3A+nDfpZH!GYY4vtZxD&Qat6 z`%m*cxr^Y;_Z4^8hrR3kNJAE7T=z^r1AE^G4tsM@H5Ah z=0E^1I8Jqz2Y%-GoxRccJ@z9JHfk#csvT^ZKpb$I%gJ2a1zMP@GsDiDZY^PbwqQ(5 zaTxM}v;KV{4dRcpV`&-Zb6qasyrUwMpsssL&{5D*L;NxJ$7NUe%Do*7A-NREf=BF@ z7096Vg4x3!FlbMSY9x62+Hf%Piz`7?mrf5%YZ7mPE^sxyK5Xg&j@t(R(;F@99#QVQcBcD)_bydyEr5@vZxZ|1}IUT;q4UtrZ_uhkmx%88Fr)?m0N}ne&voj-c(os8Zb2UAg;Li@>V^8 zuCEzyaSVd4tf`=z3CRIn#!4jMpEX;K8MM3ug~Z?W2Z9W9S!&3?HB2T^IZfb(*L=3L zl;?_5_hIKX+xerkuKlt)gc+Fxvef4W zLsyvVUk&-fU(BcTdMlvA%`Tc9t@x?9zmFb=cSguIm}$b9e7xpk;NgU*KvmccJm z#31yUsXA9C4DrvDSxvim67+Z7eX|2R91t-HooAd$ELcSSWNg2BGqMsqeTF3&dBiCC zqH7WMXXs|{=|~0dIH;u~e;6X_`AFTMXE&o6be}{)IWrT#wO+KD`(@BUiG0)3>v z!(h8{03LZ|3gMZC{>in>1&;A9MY z&{HR?R-o?Eef+Vr3jU?*7?u#j_;j5wdTJ48boZIYr{EvD9G`FN{VDy!R9PlJNZVI7 zl#f26I`BZb6r~Tb?g_R6wSO9`Ax>#GjWw);?$h*6T)dC`p(!;NeA5ZWpQn<7PSYf+ z)EHa`kMoe-%D`j%`F;06y_O3ysQ)zUEIgk>=cwn)xAdUSP`}!raT4)Dt;6DAf;^;V z<<8uXxT5OV%qH{%EZv;P0R5mMBo8X1Zc}YadS-nFeOv9gzgr!sGmm-2!+?8<{UKh0Z5SsXz>J};U1_#4h0^J}(T|1^#B z#<$e2=h(;GnQNudb-{k5u3A?*>)HIojZ#h70)dx_Xo4Si*=ZT0iW7JDd(~y z-;3DagPU`5*5TYSJ^X5IaWdFEJrmswdOM`fDuMeG-8?05?)#oN$gkz@$);S)TC$V1vJaBd<73ZOGb(SzS$~BL{E^EZw*l2|3t$&>7 z=NNK2z}Wj4{>z}g`G^+$HMS{IzzA_Qnx~!>0(*~oi0;_v3u-xE-h%aws=U4}h;#Yq z0am49_<#fcwP>yO0 z$r+R%x{#0a=MbxEe%N_J=liI-@FwKP_vo^*N}O}Q+c>2P_k!Aq22U=7rhU>1IPZTa zui?Kf4Vq6}*U|QdMdcV9v5Udl6?Z-Lwd;fR! zwa}ZPZPyoV--0Sh&phD2Avc~^ipblcFu9Ao(1oF^`|f5^p0^B^e`E#+8bWS+g2U6d zD{yWdX76oL=K^iR6J8tNv0Ta}VtNA*v61jg9@5xgY|5E9Vct0{)09H4i}iP&ZjWfNqS< zi>(ueu8vcU@+3pI#)h;tPT7JZtMk^uE@RVo=nb*1v5!8ktvp~)+@opA^B26&hy}2% zx4RAd7<-q@!_ELcKl6GH_8Xg=-IDqM4>l$rDK5o)`xeNB4w7<7*9C|SsO|B4xK2Kh)_2|QUlOv^bBZ$K(`V{`P zmEg@YH=|#GKiT?(!=MX70v)7K(DPuD0`z}6CR9!mc`*ITykIZPTnRm%QBRCv zfgaCX=@f~F?#(`39R55jYeV-O%mnF5oVs-l*hqYOdm9|u5MT}+n`3G- z8N~kPIQ$(i!mo2`y8P*gw>b;PnOchel4Y4hXTVT#yN)^Vj<)Uy^m?xUz19O|=cfP7ADT+q52nC*I#HfvzohzxKP1{9njs zUXDloUuge6m5=yaWR_+y#5jxVxc82GgWTn%+8AeXvu5{$6i|4Z>3hW4;=W5J-;l?P z@ds6Gu)oFOxa8{=sK?R_qh}B&ONs%CfzaC}i>vCU(6OZ}dR+AA3zjB~4}3#@E(^H7 ze}ww7y#B@N66|kTDq1w;7&vfH&>K3qBCpEiV+f`mI`0VGTzS2BhKurZr-50yd-&W+ zpYPVfJ}~Pc<0Hh;ij>e4sTH4Jd2?kF^Dn<-Jgq$n791$*2?o!%GRDKs%X$yq#X5qV z2F+TiCri^DW_8fLr6%?p8qlw$vdO(Qs3%Jf@`_9RplVZSD(cMQ!cX}+PB3pLJEh(& zo-}VPL3}M*#RXqJ1@PHpwx2_!&T6gU#ru+P6#8O`z30P5Qnn`yTYoW)3X*` zNwnDKEXOF}jlBE$MY!|$53qjH>i!BC!u4iJhw^jKG%Ex2F+UIWsfj^Peo7me@MHfo z4@sBnVYivI3VsdTpE)bm)sFFIoXMgwx!~6NGt5RH_0B{5w~(hFHrMqa&Zo~zo$rS3 zOq0VOnV@f*)}i^tM|sYr`$dGGQvbud8qmK?Wo=G5iF}^o>oMB~-I-$P`aJQ3@^k5< z@3Lblb)VXH6Z|*nO=>?u`FXYBxvYmApo-0+AM{~TSyX#nILP<8MV#{f>fo~LjnJjJ zvAAQf_&;HBi22!=}_@9>^A->Rkd*wn6O$?2IGvMP!f-U-^X{Vu%&#a=(z6A zijS1Oz`IKuejEExqxmV0LI(asdN8BAHVryDX3(o^k9o${KI(o=sr&CtXttnUk3PsD zxZ%Ii{CjoBBf#8&n)&KK|L*X&-`SwGcXk8V@CUbPIsZwjmA}14Z|6$(_WtWPHE933 zNn?5NukRR5O$-m}ni&5X)Z5YNzZrzl-NwPq(ZM_Z#7!hnK6pu&#rLy@QXlo$$T@ zZyyIYVOf=by+`@dHZ;^X-Dje&_3tsZ%PT3%$}36AsVGayDJsdy$t!M`|Mw>-UxwPo zy5@TOjrGl~NNx^3UJfKl*~5qb{SIX~ZPR_al)>b-|9N7%SKmk*3)A&-aB%l?wmU|W zl==5()Qxowwf{M!y@L~8!BZ(2StSKUc{vq1MOhUU1$p`Zb`E6~nSaf$xzEbM)z#V4 z+r!zO>?$euUp}C2sSisz{jrCuhr5%rkFP!EaUgs9{O1Nu_G)SGwURyjpLhT7miT{| zlI-)>jP4#@KF9vMDdqnErab)jkaZ`!Is6G(nm;$4-93GM{sfz!i_bqlMEm=1!S>HM z{1d_ddi-Y?R&Pf;oBu85?f$j`W=}5%J7;fa5BER6oJ94num5yq5Ic= zfWKSC%fnCVk8l3y0QKK?^=C#J zZ^u9B?_=v~WB=!h7FT~}wkP|Lr5wE|OZ(eQ-i|hQ9yKj&b;EeF@XODU6`hr72A+5JywGLxMU=}u&yKjSkw+dH@;bprke a_cJ(>?R-4E{_Q}jmH!1d!@=shHvj;04^wsk literal 0 HcmV?d00001 diff --git a/src/include/r_functional_space.h b/src/include/r_functional_space.h index ae71fde..6866bd1 100644 --- a/src/include/r_functional_space.h +++ b/src/include/r_functional_space.h @@ -109,7 +109,7 @@ template class R_FunctionalSpace { } // returns the coordinates of the dofs - DMatrix get_dofs_coordinates() { return fun_space_.dofs_coords(); } + DMatrix dofs_coordinates() { return fun_space_.dofs_coords(); } // destructor ~R_FunctionalSpace() = default; diff --git a/src/include/r_pde.h b/src/include/r_pde.h index 466df24..4ad5f91 100644 --- a/src/include/r_pde.h +++ b/src/include/r_pde.h @@ -90,9 +90,9 @@ template class R_PDE : public PDEWrapper { } break; case pde_type::second_order_parabolic: { SMatrix K = Rcpp::as>(pde_parameters_["diffusion"]); - DVector time_mesh = Rcpp::as>(pde_parameters_["time_mesh"]); + DVector time_nodes = Rcpp::as>(pde_parameters_["time_nodes"]); auto L = dt() + diffusion(K) + advection(b) + reaction(c); - pde_ = PDEType(domain_, time_mesh, L); + pde_ = PDEType(domain_, time_nodes, L); } break; } } @@ -104,8 +104,8 @@ template class R_PDE : public PDEWrapper { pde_.set_initial_condition(data); } // getters - DMatrix get_quadrature_nodes() const { return integrator_.quadrature_nodes(domain_); }; - DMatrix get_dofs_coordinates() const { return pde_.dof_coords(); }; + DMatrix quadrature_nodes() const { return integrator_.quadrature_nodes(domain_); }; + DMatrix dofs_coordinates() const { return pde_.dof_coords(); }; // avaiable only after initialization const SpMatrix& mass() const { return pde_.mass(); } const SpMatrix& stiff() const { return pde_.stiff(); } diff --git a/src/r_functional_space.cpp b/src/r_functional_space.cpp index d638646..82912e2 100644 --- a/src/r_functional_space.cpp +++ b/src/r_functional_space.cpp @@ -22,18 +22,18 @@ using cpp_lagrange_basis_2d_fe1 = R_FunctionalSpace<2,2,1>; RCPP_MODULE(cpp_lagrange_basis_2d_fe1) { Rcpp::class_>("cpp_lagrange_basis_2d_fe1") .constructor() - .method("size" , &R_FunctionalSpace<2,2,1>::size ) - .method("eval" , &R_FunctionalSpace<2,2,1>::eval ) - .method("integrate" , &R_FunctionalSpace<2,2,1>::integrate) - .method("get_dofs_coordinates", &R_FunctionalSpace<2,2,1>::get_dofs_coordinates); + .method("size" , &R_FunctionalSpace<2,2,1>::size ) + .method("eval" , &R_FunctionalSpace<2,2,1>::eval ) + .method("integrate" , &R_FunctionalSpace<2,2,1>::integrate ) + .method("dofs_coordinates", &R_FunctionalSpace<2,2,1>::dofs_coordinates); } using cpp_lagrange_basis_2d_fe2 = R_FunctionalSpace<2,2,2>; RCPP_MODULE(cpp_lagrange_basis_2d_fe2) { Rcpp::class_>("cpp_lagrange_basis_2d_fe2") .constructor() - .method("size" , &R_FunctionalSpace<2,2,2>::size ) - .method("eval" , &R_FunctionalSpace<2,2,2>::eval ) - .method("integrate" , &R_FunctionalSpace<2,2,2>::integrate) - .method("get_dofs_coordinates", &R_FunctionalSpace<2,2,2>::get_dofs_coordinates); + .method("size" , &R_FunctionalSpace<2,2,2>::size ) + .method("eval" , &R_FunctionalSpace<2,2,2>::eval ) + .method("integrate" , &R_FunctionalSpace<2,2,2>::integrate ) + .method("dofs_coordinates", &R_FunctionalSpace<2,2,2>::dofs_coordinates); } diff --git a/src/r_pde.cpp b/src/r_pde.cpp index 7e9fa49..3743086 100644 --- a/src/r_pde.cpp +++ b/src/r_pde.cpp @@ -22,8 +22,8 @@ using cpp_pde_2d_fe1 = R_PDE<2,2,1>; RCPP_MODULE(cpp_pde_2d_fe1) { Rcpp::class_>("cpp_pde_2d_fe1") .constructor>() - .method("get_quadrature_nodes" , &R_PDE<2,2,1>::get_quadrature_nodes ) - .method("get_dofs_coordinates" , &R_PDE<2,2,1>::get_dofs_coordinates ) + .method("quadrature_nodes" , &R_PDE<2,2,1>::quadrature_nodes ) + .method("dofs_coordinates" , &R_PDE<2,2,1>::dofs_coordinates ) .method("mass" , &R_PDE<2,2,1>::mass ) .method("stiff" , &R_PDE<2,2,1>::stiff ) .method("force" , &R_PDE<2,2,1>::force ) @@ -38,8 +38,8 @@ using cpp_pde_2d_fe2 = R_PDE<2,2,2>; RCPP_MODULE(cpp_pde_2d_fe2) { Rcpp::class_>("cpp_pde_2d_fe2") .constructor>() - .method("get_quadrature_nodes" , &R_PDE<2,2,2>::get_quadrature_nodes ) - .method("get_dofs_coordinates" , &R_PDE<2,2,2>::get_dofs_coordinates ) + .method("quadrature_nodes" , &R_PDE<2,2,2>::quadrature_nodes ) + .method("dofs_coordinates" , &R_PDE<2,2,2>::dofs_coordinates ) .method("mass" , &R_PDE<2,2,2>::mass ) .method("stiff" , &R_PDE<2,2,2>::stiff ) .method("force" , &R_PDE<2,2,2>::force ) diff --git a/tests/parabolic.R b/tests/parabolic.R index a943dd7..c9065cb 100644 --- a/tests/parabolic.R +++ b/tests/parabolic.R @@ -58,17 +58,17 @@ pde <- Pde(L, f) pde$set_boundary_condition(g) 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$set_boundary_condition(g(pde$dofs_coordinates(),times)) +pde$set_initial_condition(u0(pde$dofs_coordinates())) pde$solve() ## compute L2 norm of the error -u_ex <- exact_solution(pde$get_dofs_coordinates(), times) +u_ex <- exact_solution(pde$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] - u$coeff[,t])^2)) + error.L2[t] <- sqrt(sum(pde$mass() %*% (u_ex[,t] - u$coefficients()[,t])^2)) cat(paste0("L2 error at time ",times[t]," ", error.L2[t], "\n")) } # ------------------------------------------------------------------------------ diff --git a/tests/pde.2.R b/tests/pde.2.R index 26da401..5af2f04 100644 --- a/tests/pde.2.R +++ b/tests/pde.2.R @@ -35,8 +35,8 @@ pde$set_boundary_condition(fun=0., 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 - f$coeff)^2)) +u_ex <- as.matrix(exact_solution(pde$dofs_coordinates())) +error.L2 <- sqrt(sum(pde$mass() %*% (u_ex - f$coefficients())^2)) cat("L2 error = ", error.L2, "\n") ## perform evaluation at single point @@ -49,7 +49,7 @@ y <- x points <- expand.grid(x, y) max(abs( f$eval(points) - as.matrix(exact_solution(points)))) -evaluations <- max(abs(f$eval(mesh$get_nodes()) - exact_solution(mesh$get_nodes()))) +evaluations <- max(abs(f$eval(mesh$nodes()) - exact_solution(mesh$nodes()))) ## plot solution options(warn=-1) plot(f) %>% layout(scene=list(aspectmode="cube")) %>% hide_colorbar() @@ -59,16 +59,16 @@ contour(f) ### -basis_function <- Vh$get_basis() -Psi <- basis_function$eval(basis_function$get_dofs_coordinates()) +basis <- Vh$get_basis() +Psi <- basis$eval(basis$get_dofs_coordinates()) dim(Psi) -A <- basis_function$eval(as.matrix(points)) +A <- basis$eval(as.matrix(points)) dim(A) R0 <- pde$get_mass() R1 <- pde$get_stiff() -n <- basis_function$size() +n <- basis$size() lambda <- 1. row_1 <- cbind( 1/n * t(Psi)%*%Psi, lambda*t(R1)) row_2 <- cbind(lambda*R1, R0) @@ -76,3 +76,5 @@ row_2 <- cbind(lambda*R1, R0) M <- rbind(row_1, row_2) dim(M) + +y <- u_ex() \ No newline at end of file diff --git a/vignettes/Introduction.Rmd b/vignettes/Introduction.Rmd index fcc9b2a..102c9b4 100644 --- a/vignettes/Introduction.Rmd +++ b/vignettes/Introduction.Rmd @@ -56,8 +56,8 @@ p <- pslg(P=rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1)), S=rbind(c(1, 2), c(2, 3), c(3, 4), c(4,1))) unit_square <- triangulate(p, a = 0.00125, q=30) domain <- Mesh(unit_square) -xrange <- range(domain$get_nodes()[,1]) -yrange <- range(domain$get_nodes()[,2]) +xrange <- range(domain$nodes()[,1]) +yrange <- range(domain$nodes()[,2]) Nx <- 40 Ny <- 40 eval_x <- seq(from=xrange[1], to=xrange[2], length.out=Nx) @@ -194,10 +194,10 @@ mesh_aux <- Mesh(tpp) Vh_aux <- FunctionSpace(mesh_aux) hat_function <- Function(Vh_aux) -coeff_aux <- rep(0, times=nrow(mesh$get_nodes())) +coeff_aux <- rep(0, times=nrow(mesh_aux$nodes())) coeff_aux[5] = 1 -hat_function$set_coeff(as.matrix(coeff_aux)) +hat_function$set_coefficients(as.matrix(coeff_aux)) # scene = list(camera = list(eye = list(x = 0, y = 0, z = 1.25)), # yaxis = list(autorange="reversed"), @@ -268,8 +268,8 @@ Since, we are dealing with a simple differential problem and we know the analyti exact_solution <- function(points){ return( sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) } -u_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates())) -error.L2 <- sqrt(sum(pde$get_mass()%*%(u_ex-u$coeff)^2)) +u_ex <- as.matrix(exact_solution(pde$dofs_coordinates())) +error.L2 <- sqrt(sum(pde$mass()%*%(u_ex-u$coefficients())^2)) cat("L2 error = ", error.L2, "\n") ``` @@ -319,13 +319,13 @@ domain <- domain %X% c(0,1) mesh <- build_mesh(domain, maximum_area = 0.00125, minimum_angle = 20) # time step dT = 1e-2 -mesh$set_deltaT(dT) +mesh$set_time_step(dT) ## create Functional Space Vh <- FunctionSpace(mesh) ## define differential operator in its strong formulation u <- Function(Vh) ## define differential operator in its strong formulation -Lu <- dt(u) -laplace(u) +Lu <- dt(u) - laplace(u) ## forcing term f <- function(points,times){ res <- matrix(0, nrow=nrow(points), ncol=length(times)) diff --git a/vignettes/LakeComo.Rmd b/vignettes/LakeComo.Rmd index 59bd3f6..2bac927 100644 --- a/vignettes/LakeComo.Rmd +++ b/vignettes/LakeComo.Rmd @@ -20,7 +20,7 @@ knitr::opts_chunk$set(echo = TRUE, fig.align= "center", library(sf, quietly = T) library(mapview, quietly = T) library(femR, quietly = T) -library(rmapshaper, quietly = T) +#library(rmapshaper, quietly = T) ``` Lake Como, nestled in the Lombardy region, is recognized as a crucial freshwater ecosystem esteemed for its biodiversity and ecological significance. Environmental studies have suggested the presence and spreading of a substance of interest within the lake's waters. This hypothetical presence prompts a thorough examination to understand its nature, source, and potential implications on Lake Como's delicate ecosystem.
@@ -28,18 +28,16 @@ Assume that several monitoring stations have been strategically established arou The following interactive figure shows the comprehensive map of Lake Como and includes georeferenced black points denoting the strategic placement of monitoring stations across its expanse. The red points denote potential 'hotspots' where, in the hypothetical scenario under investigation, abnormal concentrations of the substance of interest have been recorded. ```{r, echo=FALSE} -lake_bd <- read_sf("data/deims_sites_boundariesPolygon.shp") -lake_bd_simp <- rmapshaper::ms_simplify(lake_bd, keep= 0.5, - keep_shapes=T) -st_crs(lake_bd_simp) <- 4326 -domain <- Domain(st_geometry(lake_bd_simp)) +data("lake_como_boundary", package="femR") +st_crs(lake_como_boundary) <- 4326 +domain <- Domain(st_geometry(lake_como_boundary)) mesh <- build_mesh(domain, maximum_area = 0.25e-5, minimum_angle = 20) mesh_sf <- st_as_sfc(mesh) set.seed(0) num_stations <- 20 -boundary_points <- which(mesh$get_boundary() == 1) +boundary_points <- which(mesh$boundary() == 1) boundary_id <- sample(x = 1:length(boundary_points),size = num_stations) -stations <- lapply(as.list(as.data.frame(t(mesh$get_nodes()[boundary_id,]))), +stations <- lapply(as.list(as.data.frame(t(mesh$nodes()[boundary_id,]))), st_point) stations <- st_sfc(stations, crs=4326) @@ -53,14 +51,14 @@ map.type <- "Esri.WorldTopoMap" mappa_completa <- mapview(stations, legend=F, cex=5, col.region="black", map.type=map.type) + mapview(sources, legend=F, cex=5, col.region="red") + - mapview(lake_bd_simp, col.regions = "white", lwd = 2, + mapview(lake_como_boundary, col.regions = "white", lwd = 2, alpha.regions=0.2,legend=F) mappa_completa mappa_mesh <- mapview(mesh_sf, col.region="white",alpha.regions=0.2,legend=F, alpha=0.2, lwd=0.4, map.type=map.type) + - mapview(lake_bd_simp, col.regions = "white", lwd = 2, + mapview(lake_como_boundary, col.regions = "white", lwd = 2, alpha.regions=0.2,legend=F) beta <- lapply(st_geometry(stations)[c(9,8)], as.numeric) @@ -92,16 +90,16 @@ library(femR, quietly = T) library(rmapshaper, quietly = T) lake_bd <- read_sf("data/deims_sites_boundariesPolygon.shp") -lake_bd_simp <- rmapshaper::ms_simplify(lake_bd, keep= 0.5, +lake_como_boundary <- rmapshaper::ms_simplify(lake_bd, keep= 0.5, keep_shapes=T) -st_crs(lake_bd_simp) <- 4326 +st_crs(lake_como_boundary) <- 4326 ``` The `femR` package includes utilities to read a `sfc` (Simple Features Collection) object, such as the boundary of the Lake, and generating meshes by specifying the `maximum_area` of each triangle and the `minimum_angle` between adjacent sides of the triangles. ```{r} # create the physical domain -domain <- Domain(st_geometry(lake_bd_simp)) +domain <- Domain(st_geometry(lake_como_boundary)) # create the mesh mesh <- build_mesh(domain, maximum_area = 0.25e-5, minimum_angle = 20) ``` @@ -111,7 +109,7 @@ Furthermore, the `Mesh` class overloads the `st_as_sfc` method, enabling visuali mesh_sf <- st_as_sfc(mesh) map.type <- "Esri.WorldTopoMap" -mapview(lake_bd_simp, col.regions = "white", lwd = 2, +mapview(lake_como_boundary, col.regions = "white", lwd = 2, alpha.regions=0.2, legend=F, map.type=map.type) + mapview(mesh_sf, col.region="white", alpha.regions=0.2, legend=F, alpha=0.4, lwd=0.4, map.type=map.type) @@ -159,9 +157,9 @@ pde$solve() ``` It is possible to visualize the computed solution on an interactive map using the `mapview` package as the following code snippet shows. ```{r} -coeff <- apply(mesh$get_elements(), MARGIN=1, +coeff <- apply(mesh$elements(), MARGIN=1, FUN= function(edge){ - mean(u$coeff[edge]) + mean(u$coefficients()[edge]) }) U <- st_sf(data.frame(coeff = coeff), geometry= mesh_sf) mapview(U, alpha.regions=0.8, alpha = 0, map.type=map.type) @@ -170,8 +168,8 @@ mapview(U, alpha.regions=0.8, alpha = 0, map.type=map.type) ```{r,eval=FALSE, include=FALSE} library(webshot2) library(leaflet) -cntr_crds <- c(mean(mesh$get_nodes()[, 1]), - mean(mesh$get_nodes()[, 2])) +cntr_crds <- c(mean(mesh$nodes()[, 1]), + mean(mesh$nodes()[, 2])) solution <- mapview(U, alpha.regions=0.8, legend=F, alpha = 0, @@ -260,7 +258,7 @@ g <- function(points, times){ return(res) } # initial condition -u0 <- g(mesh$get_nodes(), times[1]) +u0 <- g(mesh$nodes(), times[1]) # build PDE object pde <- Pde(Lu, 0.) # set initial conditions