Skip to content

Commit

Permalink
Introduced private attributes within R6 classes.
Browse files Browse the repository at this point in the history
  • Loading branch information
@aldoclemente
aldoclemente committed Jan 8, 2024
1 parent 656e45c commit eecdb9a
Show file tree
Hide file tree
Showing 9 changed files with 284 additions and 215 deletions.
82 changes: 37 additions & 45 deletions R/domain.R
View file
Original file line number Diff line number Diff line change
@@ -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
}
)
)
Expand Down Expand Up @@ -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")
Expand All @@ -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
Expand Down Expand Up @@ -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}.
Expand All @@ -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)
})

Expand All @@ -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
})
106 changes: 59 additions & 47 deletions R/function_space.R
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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))
}
)
)
Expand All @@ -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)
}
)
)

Expand Down Expand Up @@ -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
}
)
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
Loading

0 comments on commit eecdb9a

Please sign in to comment.