Skip to content

Commit

Permalink
Updated R interface.
Browse files Browse the repository at this point in the history
  • Loading branch information
@aldoclemente
aldoclemente committed Jan 15, 2024
1 parent eecdb9a commit 252bcd0
Show file tree
Hide file tree
Showing 19 changed files with 297 additions and 311 deletions.
1 change: 0 additions & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,5 @@ Encoding: UTF-8
Suggests:
knitr,
rmarkdown,
rmapshaper,
mapview
VignetteBuilder: knitr
1 change: 0 additions & 1 deletion NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
9 changes: 8 additions & 1 deletion R/datasets.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,11 @@
#' This dataset can be used to create a mesh object with the function \code{Mesh}.
#'
#' @name unit_square
NULL
NULL

#' Lake Como boundary
#'
#' sf object storing the boundary of Lake Como.
#'
#' @name lake_como_boundary
NULL
44 changes: 22 additions & 22 deletions R/domain.R
View file
Original file line number Diff line number Diff line change
@@ -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
}
)
)
Expand Down Expand Up @@ -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)
})

Expand All @@ -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
})
134 changes: 73 additions & 61 deletions R/function_space.R
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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))
}
)
)
Expand All @@ -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_
}
)
)
Expand Down Expand Up @@ -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)
)
}
Expand All @@ -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_
}
)
)
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
Loading

0 comments on commit 252bcd0

Please sign in to comment.