diff --git a/R/domain.R b/R/domain.R index c69fcc4..58c0686 100644 --- a/R/domain.R +++ b/R/domain.R @@ -1,13 +1,17 @@ .DomainCtr <- R6Class("Domain", - public = list( + private = list( geometry = "ANY", - coords = "data.frame", # (x,y,boundary) + coords = "matrix", # (x,y) + boundary = "matrix", # time_interval = vector(mode="numeric", length = 0), - crs = "ANY", - initialize = function(geometry, coords, crs){ - self$geometry <- geometry - self$coords <- coords - self$crs <- crs + crs = "ANY" + ), + public = list( + initialize = function(geometry, coords, boundary, crs){ + private$geometry <- geometry + private$coords <- coords + private$boundary <- boundary + private$crs <- crs } ) ) @@ -51,9 +55,10 @@ setMethod("Domain", signature = "list", function(x){ edges=x$edges, edges_group=x$edges_group, edges_boundary=x$edges_boundary, edges_ring = x$edges_ring) - coords <- data.frame(x=x$nodes[,1], y=x$nodes[,2], boundary=x$nodes_boundary) - crs = NA - .DomainCtr$new(geometry = geometry, coords=coords, crs = crs) + coords <- cbind(x$nodes[,1], x$nodes[,2]) + boundary <- x$nodes_boundary + crs = NA + .DomainCtr$new(geometry = geometry, coords=coords, boundary = boundary, crs = crs) }) setOldClass("pslg") @@ -78,9 +83,10 @@ setMethod("Domain", signature = "pslg", function(x){ edges=edges, edges_group=edges_group, edges_boundary=edges_boundary, edges_ring=edges_ring) - coords <- data.frame(x=nodes[,1], y=nodes[,2], boundary=nodes_boundary) + coords <- cbind(nodes[,1], nodes[,2]) + boundary <- nodes_boundary - .DomainCtr$new(geometry = geometry, coords=coords, crs = NA) + .DomainCtr$new(geometry = geometry, coords=coords, boundary=boundary, crs = NA) }) #' @rdname Domain @@ -216,32 +222,18 @@ setMethod("Domain", signature="sfc", function(x){ edges_ring = edges_ring) } - # da salvare + coords <- unique(st_coords[order(st_coords[,(ncol(st_coords)-1)]),][,c(1:2,ncol(st_coords))]) - coords <- as.data.frame(coords) - storage.mode(coords$boundary) <- "integer" + + boundary <- coords[,3] + storage.mode(boundary) <- "integer" + coords <- coords[,1:2] crs <- NA if(!is.na(st_crs(x))) crs <- st_crs(x)$input - .DomainCtr$new(geometry = geometry, coords=coords, crs = crs) + .DomainCtr$new(geometry = geometry, coords=coords, boundary=boundary, crs = crs) }) -# setMethod("Domain", signature = "sfc", function(x){ -# if(!any(class(x) %in% c("sfc_POLYGON", "sfc_MULTIPOLYGON"))) -# stop("Error!") -# x_boundary <- x %>% st_union() # return a sfc_POLYGON with 1 feature -# nodes <- unique(st_coordinates(x_boundary))[,1:2] -# edges <- cbind(1:nrow(nodes),c(2:nrow(nodes),1)) -# nodes_group <- rep(1, times=nrow(nodes)) -# edges_group <- rep(1, times=nrow(edges)) -# geometry <- list(coords = nodes, edges = edges, -# nodes_group = nodes_group, edges_group = edges_group) -# crs <- NA_crs_ -# if(!is.na(st_crs(x))) crs <- st_crs(x) -# .DomainCtr$new(geometry = geometry, time_interval = vector(mode="numeric", length = 0), -# coords=geometry$nodes, crs = crs) -# }) - #' Delaunay triangulation of the spatial domain #' #' @param domain a domain object returned by \code{Domain}. @@ -262,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(domain$geometry)){ - groups_id <- unique(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(domain$geometry[[sub_id]]$edges_ring==holes_id[i]) - if( ! all(as.integer(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(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(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, domain$geometry[[sub_id]]$edges) - SB = rbind(SB, as.matrix(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 = as.matrix(domain$coords[,1:2]), PB = as.matrix(domain$coords[,3]), + 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) - res$geometry <- domain$geometry - res$time_interval <- domain$time_interval - res$crs <- domain$crs + set_geometry(res, extract_private(domain)$geometry) + set_time_interval(res, extract_private(domain)$time_interval) + set_crs(res, extract_private(domain)$crs) return(res) }) @@ -312,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.") - op1$time_interval <- op2 + set_time_interval(op1, op2) op1 }) \ No newline at end of file diff --git a/R/function_space.R b/R/function_space.R index 18eaeac..b9ea828 100644 --- a/R/function_space.R +++ b/R/function_space.R @@ -1,20 +1,22 @@ .BasisFunctionCtr <- R6Class("BasisFunction", - public = list( + private = list( cpp_handler = "ANY", ## pointer to cpp backend - mesh = "Mesh", + mesh = "Mesh" + ), + public = list( initialize = function(cpp_handler, mesh){ - self$cpp_handler = cpp_handler - self$mesh = mesh + private$cpp_handler = cpp_handler + private$mesh = mesh }, size = function(){ - self$cpp_handler$size() + private$cpp_handler$size() }, get_dofs_coordinates = function(){ - self$cpp_handler$get_dofs_coordinates() + private$cpp_handler$get_dofs_coordinates() }, ## evaluates the basis system on the given set of locations eval = function(locations) { - return(self$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 @@ -23,9 +25,9 @@ ## required dof coordinates... c <- func if (is.function(func)) { ## recover fem expansion - c <- func(self$cpp_handler$get_dofs_coordinates()) + c <- func(private$cpp_handler$get_dofs_coordinates()) } - return(self$cpp_handler$integrate(c)) + return(private$cpp_handler$integrate(c)) } ) ) @@ -43,22 +45,26 @@ setMethod("BasisFunction", signature = signature("Mesh","integer"), function(mesh, fe_order){ 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 = ""))), extract_private(mesh)$cpp_handler, 0), mesh=mesh) }) #' @export .FunctionSpaceCtr <- R6Class("FunctionSpace", - public = list( - mesh = "ANY", + private= list( + mesh = "Mesh", basis_function = "BasisFunction", ## cpp backend - fe_order = "integer", ## this is specific for fem, must be generalized + fe_order = "integer" ## this is specific for fem, must be generalized + ), + public = list( initialize = function(mesh, basis_function, fe_order){ - self$mesh <- mesh - self$basis_function <- basis_function - self$fe_order <- fe_order + private$mesh <- mesh + private$basis_function <- basis_function + private$fe_order <- fe_order }, - get_basis = function() { return(self$basis_function) } + get_basis = function(){ + return(private$basis_function) + } ) ) @@ -114,34 +120,38 @@ setMethod("FunctionSpace", signature = c(mesh="Mesh", fe_order="missing"), ## finite element function .FunctionCtr <- R6Class("Function", + private = list( + FunctionSpace = "FunctionSpace" + ), public = list( - FunctionSpace = "FunctionSpace", coeff = "matrix", initialize = function(FunctionSpace, coeff){ - self$FunctionSpace <- FunctionSpace + private$FunctionSpace <- FunctionSpace self$coeff <- coeff }, eval = function(X) { - M = dim(self$FunctionSpace$mesh$get_nodes())[2] + M = dim(extract_private(private$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, + if(length(extract_private(private$FunctionSpace)$mesh$get_times()) == 0){ + evals <- apply(private$FunctionSpace$get_basis()$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], + 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], MARGIN=1, FUN=sum) } } return(evals) }, - set_coeff = function(coefficients){ + 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 } ) @@ -185,16 +195,17 @@ Function <- function(FunctionSpace) { #' @importFrom graphics plot #' @export plot.Function <- function(x, ...){ - times <- x$FunctionSpace$mesh$times + mesh <- extract_private(extract_private(x)$FunctionSpace)$mesh + times <- mesh$get_times() is_parabolic <- FALSE if(length(times)!=0) is_parabolic <- TRUE if(!is_parabolic){ - plot_data <- data.frame(X=x$FunctionSpace$mesh$get_nodes()[,1], - Y=x$FunctionSpace$mesh$get_nodes()[,2], - coeff=x$coeff[1:nrow(x$FunctionSpace$mesh$get_nodes())]) - I=x$FunctionSpace$mesh$get_elements()[,1]-1 - J=x$FunctionSpace$mesh$get_elements()[,2]-1 - K=x$FunctionSpace$mesh$get_elements()[,3]-1 + 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 fig<- plot_ly(plot_data, x=~X, y=~Y, z=~coeff, i = I, j = J, k = K, intensity=~coeff, color=~coeff, type="mesh3d", @@ -211,24 +222,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(x$FunctionSpace$mesh$get_nodes()[,1], times=length(times)), - Y=rep(x$FunctionSpace$mesh$get_nodes()[,2], times=length(times)), + plot_data <- data.frame(X=rep(mesh$get_nodes()[,1], times=length(times)), + Y=rep(mesh$get_nodes()[,2], times=length(times)), Z=rep(0, times=length(times)), - coeff=as.vector(x$coeff[1:nrow(x$FunctionSpace$mesh$get_nodes()),]), - times = round(rep(times, each=nrow(x$FunctionSpace$mesh$get_nodes())),3)) - coeff <- matrix(nrow=nrow(x$FunctionSpace$mesh$get_elements()), + 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()), ncol=length(times)) for(t in 1:length(times)){ - coeff[,t] <- apply(x$FunctionSpace$mesh$get_elements(), MARGIN=1, FUN= + coeff[,t] <- apply(mesh$get_elements(), MARGIN=1, FUN= function(edge){ mean(x$coeff[edge,t]) }) } limits = c(min(coeff), max(coeff)) cmin = limits[1]; cmax=limits[2] - I=x$FunctionSpace$mesh$get_elements()[,1]-1 - J=x$FunctionSpace$mesh$get_elements()[,2]-1 - K=x$FunctionSpace$mesh$get_elements()[,3]-1 + I=mesh$get_elements()[,1]-1 + J=mesh$get_elements()[,2]-1 + K=mesh$get_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", @@ -259,13 +270,14 @@ plot.Function <- function(x, ...){ #' @return A plotly object #' @export setMethod("contour", signature=c(x="Function"), function(x, ...){ - times <- x$FunctionSpace$mesh$times + mesh <- extract_private(extract_private(x)$FunctionSpace)$mesh + times <- mesh$get_times() is_parabolic <- FALSE if(length(times)!=0) is_parabolic <- TRUE if(!is_parabolic){ - xrange <- range(x$FunctionSpace$mesh$get_nodes()[,1]) - yrange <- range(x$FunctionSpace$mesh$get_nodes()[,2]) + xrange <- range(mesh$get_nodes()[,1]) + yrange <- range(mesh$get_nodes()[,2]) Nx <- 40 Ny <- 40 eval_x <- seq(from=xrange[1], to=xrange[2], length.out=Nx) @@ -280,8 +292,8 @@ setMethod("contour", signature=c(x="Function"), function(x, ...){ yaxis = list(title = "", showgrid=F, zeroline=F, ticks="", showticklabels=F)) }else{ - xrange <- range(x$FunctionSpace$mesh$get_nodes()[,1]) - yrange <- range(x$FunctionSpace$mesh$get_nodes()[,2]) + xrange <- range(mesh$get_nodes()[,1]) + yrange <- range(mesh$get_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 cb94265..afee0bc 100644 --- a/R/mesh.R +++ b/R/mesh.R @@ -2,37 +2,39 @@ ## 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 + ), 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 + private$cpp_handler <- cpp_handler + private$m <- m + private$n <- n }, get_nodes = function(){ - self$cpp_handler$nodes() + private$cpp_handler$nodes() }, get_elements = function(){ - self$cpp_handler$elements()+1 + private$cpp_handler$elements()+1 }, get_boundary = function(){ - self$cpp_handler$boundary() + private$cpp_handler$boundary() }, get_neighbors = function(){ - self$cpp_handler$neighbors() + private$cpp_handler$neighbors() }, get_times = function(){ - return(self$times) + private$times }, - set_deltaT = function(deltaT){ - if(length(self$time_interval)==0) + set_time_step = function(time_step){ + if(length(private$time_interval)==0) stop("Error!") - self$deltaT <- deltaT - self$times <- seq(self$time_interval[1], self$time_interval[2], by=deltaT) + private$time_step <- time_step + private$times <- seq(private$time_interval[1], private$time_interval[2], by=time_step) } ) ) @@ -112,44 +114,50 @@ 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.") - op1$times <- times - op1$time_interval <- c(times[1], times[length(times)]) - op1$deltaT <- times[2] - times[1] + set_times(op1, op2) + set_time_interval(op1, c(op2[1], op2[length(op2)])) + set_time_step(op1, (op2[2] - op2[1])) op1 }) ## Mesh - auxiliary methods -unroll_edges_aux <- function(Mesh){ - mesh <- Mesh$cpp_handler - 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)] + 1 - edges[(3*(i-1) + 2),] = mesh$elements()[i,c(2,3)] + 1 - edges[(3*(i-1) + 3),] = mesh$elements()[i,c(3,1)] + 1 +# 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 +# } + +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 -} - -setGeneric("unroll_edges", function(Mesh) standardGeneric("unroll_edges")) -setMethod("unroll_edges", "Mesh", function(Mesh){ - unroll_edges_aux(Mesh) }) plot_mesh_aux <- function(x, ...){ edges <- unroll_edges(x) plot_ly(...) %>% - add_markers(x = x$cpp_handler$nodes()[,1], - y = x$cpp_handler$nodes()[,2], + add_markers(x = x$get_nodes()[,1], + y = x$get_nodes()[,2], color = I('black'), size = I(1), hoverinfo = 'text', - text = paste('
Coordinates:', round(x$cpp_handler$nodes()[,1],2), - round(x$cpp_handler$nodes()[,2],2)), + text = paste('
Coordinates:', round(x$get_nodes()[,1],2), + round(x$get_nodes()[,2],2)), showlegend = T, visible = T) %>% - add_segments(x = x$cpp_handler$nodes()[edges[,1],1], - y = x$cpp_handler$nodes()[edges[,1],2], - xend = x$cpp_handler$nodes()[edges[,2],1], - yend = x$cpp_handler$nodes()[edges[,2],2], + 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], color = I('black'), size = I(1), showlegend = F) %>% layout( @@ -167,11 +175,6 @@ plot_mesh_aux <- function(x, ...){ )) } - -# setMethod("plot", signature=c(x="Mesh"), function(x, ...){ -# plot_mesh_aux(x, ...) -# }) - #' Plot a Mesh object #' #' @param x A \code{Mesh} object defining the triangular mesh, as generated by \code{Mesh} diff --git a/R/pde.R b/R/pde.R index 178096c..cdc1d51 100644 --- a/R/pde.R +++ b/R/pde.R @@ -3,82 +3,98 @@ pde_type_list <- list("laplacian" = 1, "elliptic" = 2, "parabolic" = 3) # Pde Class wraps C++ R_PDE class .PdeCtr <- R6Class("Pde", - public = list( + private = list( is_dirichletBC_set = "logical", is_initialCondition_set = "logical", is_parabolic ="logical", cpp_handler = "ANY", - DifferentialOp = "DiffOpObject", + DifferentialOp = "DiffOpObject" + ), + public = list( 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 + 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 }, solve = function(){ - 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(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(!private$is_dirichletBC_set){ + if(!private$is_parabolic){ + dirichletBC_ = as.matrix(rep(0,times=nrow(private$cpp_handler$get_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()), + ncol=length(times)) + } + private$cpp_handler$set_dirichlet_bc(dirichletBC_) + private$is_dirichletBC_set <- TRUE } - if(self$is_parabolic & (!self$is_initialCondition_set)){ + if(private$is_parabolic & (!private$is_initialCondition_set)){ stop("initialCondition must be provided.") } - self$cpp_handler$solve() + private$cpp_handler$solve() + private$DifferentialOp$f$coeff <- private$cpp_handler$solution() + invisible(self) }, get_dofs_coordinates = function(){ - self$cpp_handler$get_dofs_coordinates() + private$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(any(typeof(fun) == c("function", "closure"))){ - if(!self$is_parabolic){ - dirichletBC_ <- as.matrix(fun(self$cpp_handler$get_dofs_coordinates())) + if(!private$is_parabolic){ + dirichletBC_ <- as.matrix(fun(private$cpp_handler$get_dofs_coordinates())) }else{ - dirichletBC_ <- fun(self$cpp_handler$get_dofs_coordinates(), - self$DifferentialOp$f$FunctionSpace$mesh$times) + times <- extract_private( + extract_private( + private$DifferentialOp$f)$FunctionSpace)$mesh$get_times() + dirichletBC_ <- fun(private$cpp_handler$get_dofs_coordinates(), times) } }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){ - if(nrow(as.matrix(fun)) == nrow(self$cpp_handler$get_dofs_coordinates())){ + if(nrow(as.matrix(fun)) == nrow(private$cpp_handler$get_dofs_coordinates())){ dirichletBC_ <- fun }else if (nrow(as.matrix(fun)) == 1L){ - if(!self$is_parabolic) - dirichletBC_ <- matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), ncol=1) - else - dirichletBC_ <- matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), - ncol=length(self$DifferentialOp$f$FunctionSpace$mesh$times)) - } + if(!private$is_parabolic) + dirichletBC_ <- matrix(fun, nrow=nrow(private$cpp_handler$get_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()), + ncol=length(times)) + } + } } - self$is_dirichletBC_set <- TRUE - self$cpp_handler$set_dirichlet_bc(dirichletBC_) + private$is_dirichletBC_set <- TRUE + private$cpp_handler$set_dirichlet_bc(dirichletBC_) + invisible(self) }, set_initial_condition = function(fun){ - if(!self$is_parabolic) + if(!private$is_parabolic) stop("Cannot set initial condition for elliptic problem.") - self$is_initialCondition_set <- TRUE + private$is_initialCondition_set <- TRUE if(typeof(fun) == "closure"){ - self$cpp_handler$set_initial_condition(fun(self$cpp_handler$get_dofs_coordinates())) + private$cpp_handler$set_initial_condition(fun(private$cpp_handler$get_dofs_coordinates())) }else if(any(typeof(fun) == c("matrix","vector", "numeric" ,"double"))){ - if(nrow(as.matrix(fun)) == nrow(self$cpp_handler$get_dofs_coordinates())){ - self$cpp_handler$set_initial_condition(fun) + if(nrow(as.matrix(fun)) == nrow(private$cpp_handler$get_dofs_coordinates())){ + private$cpp_handler$set_initial_condition(fun) }else if(nrow(as.matrix(fun)) == 1L){ - self$cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(self$cpp_handler$get_dofs_coordinates()), + private$cpp_handler$set_initial_condition(matrix(fun, nrow=nrow(private$cpp_handler$get_dofs_coordinates()), ncol=1)) } } + invisible(self) }, get_mass = function(){ - self$cpp_handler$mass() + private$cpp_handler$mass() }, get_stiff = function(){ - self$cpp_handler$stiff() + private$cpp_handler$stiff() } ) ) @@ -113,7 +129,8 @@ parse_pde_parameters <- function(DifferentialOp){ } if(pde_type == pde_type_list$parabolic) - pde_parameters[["time_mesh"]] <- as.vector(DifferentialOp$f$FunctionSpace$mesh$times) + pde_parameters[["time_mesh"]] <- as.vector( + extract_private(extract_private(DifferentialOp$f)$FunctionSpace)$mesh$get_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) @@ -129,16 +146,20 @@ make_pde <- function(DifferentialOp) { pde_parameters <- parse_pde_parameters(DifferentialOp) ## define Rcpp module - D <- DifferentialOp$f$FunctionSpace$mesh$cpp_handler ## domain - fe_order <- DifferentialOp$f$FunctionSpace$fe_order ## finite element order + 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 cpp_handler <- NULL if (fe_order == 1) { ## linear finite elements - cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters, DifferentialOp$f) + cpp_handler <- new(cpp_pde_2d_fe1, D, pde_type - 1, pde_parameters) } if (fe_order == 2) { ## quadratic finite elements - cpp_handler <- new(cpp_pde_2d_fe2, D, pde_type - 1, pde_parameters, DifferentialOp$f) + cpp_handler <- new(cpp_pde_2d_fe2, D, pde_type - 1, pde_parameters) } - return(cpp_handler) } @@ -161,7 +182,8 @@ setMethod("Pde", signature=c(DifferentialOp="DiffOpObject", f="function"), quad_nodes <- cpp_handler$get_quadrature_nodes() ## evaluate forcing term on quadrature nodes if(pde_type == pde_type_list$parabolic) { - times <- DifferentialOp$f$FunctionSpace$mesh$times + times <- extract_private( + extract_private(DifferentialOp$f)$FunctionSpace)$mesh$get_times() 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 544aa85..72a2d39 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(x$geometry)){ + for(sub_id in 1:length(extract_private(x)$geometry)){ - geometry <- 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(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(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 <- 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(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(x$geometry)) - for(sub_id in 1:length(x$geometry)){ - geometry <- 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(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 <- x$crs + crs <- extract_private(x)$crs if(is.na(crs)) crs <- NA_crs_ - if(length(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(x){ - st_cast(st_linestring(mesh$get_nodes()[x,]), + polygon_list <- apply(x$get_elements(), MARGIN=1, FUN=function(elem){ + st_cast(st_linestring(x$get_nodes()[elem,]), to ="POLYGON" ) }) - crs <- x$crs + crs <- extract_private(x)$crs if(is.na(crs)) crs <- NA_crs_ mesh_sf <- st_sfc(polygon_list, crs = crs) mesh_sf @@ -117,7 +117,6 @@ st_as_sfc.Mesh <- function(x, ...){ #' @rdname sf #' @export st_crs.Domain <- function(x, ...){ - #st_crs(x$crs, ...) st_crs(st_as_sfc(x), ...) } @@ -125,9 +124,9 @@ st_crs.Domain <- function(x, ...){ #' @rdname sf #' @export `st_crs<-.Domain` = function(x, value) { - #st_crs(x$crs) = value sfc <- st_as_sfc(x) st_crs(sfc) <- value - crs <- st_crs(sfc)$input + set_crs(x, value) + #crs <- st_crs(sfc)$input } diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..5196860 --- /dev/null +++ b/R/utils.R @@ -0,0 +1,41 @@ +# returns the private attributes of the R6 class passed as parameter +setGeneric("extract_private", function(x) standardGeneric("extract_private")) +setMethod("extract_private", signature = "R6", + function(x){ + 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) +}) + +# extract_private <- function(x){ +# if(!is(x,"R6")) stop("Input parameter must be an R6 class.") +# x$.__enclos_env__$private +# } diff --git a/src/include/r_pde.h b/src/include/r_pde.h index 554f7f9..466df24 100644 --- a/src/include/r_pde.h +++ b/src/include/r_pde.h @@ -63,11 +63,10 @@ template class R_PDE : public PDEWrapper { DomainType domain_ {}; // triangulation QuadratureRule integrator_ {}; // quadrature rule (exact for the provided fem order) int pde_type_; // one of the pde_type enum values - Rcpp::Environment solution_; public: // constructor - R_PDE(Rcpp::Environment mesh, int pde_type, const Rcpp::Nullable& pde_parameters, Rcpp::Environment solution) : - pde_type_(pde_type), solution_(solution){ + R_PDE(Rcpp::Environment mesh, int pde_type, const Rcpp::Nullable& pde_parameters) : + pde_type_(pde_type){ // set domain SEXP meshptr = mesh[".pointer"]; R_Mesh* ptr = reinterpret_cast*>(R_ExternalPtrAddr(meshptr)); @@ -114,10 +113,9 @@ template class R_PDE : public PDEWrapper { // initialize internal pde status void init() { pde_.init(); } //solve - void solve() { - pde_.solve(); - solution_["coeff"] = pde_.solution(); - } + void solve() { pde_.solve(); } + // returns solution + const DMatrix& solution() const { return pde_.solution(); } // destructor ~R_PDE() = default; }; diff --git a/src/r_pde.cpp b/src/r_pde.cpp index 2039ffc..7e9fa49 100644 --- a/src/r_pde.cpp +++ b/src/r_pde.cpp @@ -21,7 +21,7 @@ using cpp_pde_2d_fe1 = R_PDE<2,2,1>; RCPP_MODULE(cpp_pde_2d_fe1) { Rcpp::class_>("cpp_pde_2d_fe1") - .constructor, Rcpp::Environment>() + .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("mass" , &R_PDE<2,2,1>::mass ) @@ -31,12 +31,13 @@ RCPP_MODULE(cpp_pde_2d_fe1) { .method("set_forcing" , &R_PDE<2,2,1>::set_forcing ) .method("set_initial_condition", &R_PDE<2,2,1>::set_initial_condition) .method("init" , &R_PDE<2,2,1>::init ) - .method("solve" , &R_PDE<2,2,1>::solve ); + .method("solve" , &R_PDE<2,2,1>::solve ) + .method("solution" , &R_PDE<2,2,1>::solution ); } using cpp_pde_2d_fe2 = R_PDE<2,2,2>; RCPP_MODULE(cpp_pde_2d_fe2) { Rcpp::class_>("cpp_pde_2d_fe2") - .constructor, Rcpp::Environment>() + .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("mass" , &R_PDE<2,2,2>::mass ) @@ -46,5 +47,6 @@ RCPP_MODULE(cpp_pde_2d_fe2) { .method("set_forcing" , &R_PDE<2,2,2>::set_forcing ) .method("set_initial_condition", &R_PDE<2,2,2>::set_initial_condition) .method("init" , &R_PDE<2,2,2>::init ) - .method("solve" , &R_PDE<2,2,2>::solve ); + .method("solve" , &R_PDE<2,2,2>::solve ) + .method("solution" , &R_PDE<2,2,2>::solution ); } diff --git a/tests/parabolic.R b/tests/parabolic.R index 6ce7336..a943dd7 100644 --- a/tests/parabolic.R +++ b/tests/parabolic.R @@ -12,7 +12,7 @@ times = seq(t0,t_max,length=M) # Spatio-temporal domain ( Mesh inherits from DOMAIN :) ) mesh = Mesh(unit_square) %X% c(t0, t_max) #times -mesh$set_deltaT(dT) # NUOVO +mesh$set_time_step(dT) # NUOVO class(mesh) plot(mesh) @@ -31,9 +31,9 @@ exact_solution <- function(points,t){ u <- Function(Vh) ## define differential operator in its strong formulation -#L <- dt(f) + (-1)*laplace(f) +L <- dt(u) -laplace(u) -L <- dt(u) - laplace(u) + dot(c(0,1),grad(u)) +#L <- dt(u) - laplace(u) + dot(c(0,1),grad(u)) ## forcing term f <- function(points,times){