Skip to content

Commit

Permalink
Updated submodules. Added BasisFunction classs.
Browse files Browse the repository at this point in the history
aldoclemente committed Dec 21, 2023
1 parent a07a557 commit 3e05f37
Showing 16 changed files with 598 additions and 427 deletions.
88 changes: 71 additions & 17 deletions R/function_space.R
Original file line number Diff line number Diff line change
@@ -1,8 +1,59 @@
# .CppFunctionSpace <- setRefClass(
# Class = "CppFunctionSpace",
# fields = c(
# cpp_handler = "ANY" ## pointer to cpp backend
# ),
# methods = c(
# ## evaluates the basis system on the given set of locations
# eval = function(type = c("pointwise", "areal"), locations) {
# type <- match.arg(type)
# return(cpp_handler$eval(match(type, c("pointwise", "areal")) - 1, locations))
# }
# )
# )

.BasisFunctionCtr <- setRefClass(
Class = "BasisFunction",
fields = c(
cpp_handler = "ANY" ## pointer to cpp backend
),
methods = c(
## evaluates the basis system on the given set of locations
eval = function(locations) {
return(cpp_handler$eval(0L, locations))
},
## integrates a function expressed as basis expansion with respect to this basis system over the
## whole domain of definition
integrate = function(func) {
if (!(is.function(func) || is.vector(func))) stop("invalid argument, should be either a function or a vector")
## required dof coordinates...
c <- func
if (is.function(func)) { ## recover fem expansion
c <- func(mesh$cpp_handler$nodes())
}
return(cpp_handler$integrate(c))
}
)
)

setGeneric("BasisFunction", function(mesh, fe_order) standardGeneric("BasisFunction"))

setMethod("BasisFunction", signature = signature("Mesh","integer"),
function(mesh, fe_order){
basis_function <- .BasisFunctionCtr(cpp_handler =
new(eval(parse(text = paste("cpp_lagrange_basis_2d_fe",
as.character(fe_order), sep = ""))), mesh$cpp_handler, 0))
})

.FunctionSpaceCtr <- setRefClass(
Class = "FunctionSpace",
Class = "FunctionSpaceObject",
fields = c(
mesh = "ANY",
fe_order="integer"
mesh = "ANY",
basis_function = "BasisFunction", ## cpp backend
fe_order = "integer" ## this is specific for fem, must be generalized
),
methods = c(
get_basis = function() { return(basis_function) }
)
)

@@ -23,15 +74,16 @@ setGeneric("FunctionSpace", function(mesh,fe_order) standardGeneric("FunctionSpa
#' mesh <- Mesh(unit_square)
#' Vh <- FunctionSpace(mesh = mesh, fe_order = 1)
#' }
setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="numeric"),
function(mesh,fe_order){
if(fe_order == 1){
return(.FunctionSpaceCtr(mesh=mesh, fe_order=1L))
}else if(fe_order == 2){
return(.FunctionSpaceCtr(mesh=mesh, fe_order=2L))
}

})
setMethod("FunctionSpace",
signature = c(mesh = "ANY", fe_order = "numeric"),
function(mesh, fe_order) {
.FunctionSpaceCtr(
basis_function = BasisFunction(mesh,fe_order),
mesh = mesh,
fe_order = as.integer(fe_order)
)
}
)

#' @rdname FunctionSpace
#' @examples
@@ -55,7 +107,7 @@ setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
pde = "ANY"
),
methods = list(
eval_at = function(X) {
eval = function(X) {
M = dim(FunctionSpace$mesh$get_nodes())[2]
if(is.vector(X)) X <- t(as.matrix(X))

@@ -64,11 +116,13 @@ setMethod("FunctionSpace", signature = c(mesh="ANY", fe_order="missing"),
}
evals <- NULL
if(length(FunctionSpace$mesh$times) == 0){
evals <- pde$eval(FunctionSpace$mesh$data, coeff, as.matrix(X))
evals <- apply(FunctionSpace$basis_function$eval(as.matrix(X)) %*% coeff,
MARGIN=1, FUN=sum)
}else{
evals <- matrix(nrow=nrow(X),ncol=length(FunctionSpace$mesh$times))
for(t in 1:length(FunctionSpace$mesh$times)){
evals[,t] <- pde$eval(FunctionSpace$mesh$data, as.matrix(coeff[,t]), as.matrix(X))
evals[,t] <- apply(FunctionSpace$basis_function$eval(as.matrix(X)) %*% coeff[,t],
MARGIN=1, FUN=sum)
}
}
return(evals)
@@ -196,7 +250,7 @@ setMethod("contour", signature=c(x="Function"), function(x, ...){
eval_x <- seq(from=xrange[1], to=xrange[2], length.out=Nx)
eval_y <- seq(from=yrange[1], to=yrange[2], length.out=Ny)
eval_points <- expand.grid(eval_x, eval_y)
Z <- matrix(x$eval_at(eval_points), nrow=Nx,ncol=Ny)
Z <- matrix(x$eval(eval_points), nrow=Nx,ncol=Ny)
fig <- plot_ly(type="contour", x=eval_x, y=eval_y, z=Z,
intensity=Z, color = Z,
contours=list(showlabels = TRUE),
@@ -218,7 +272,7 @@ setMethod("contour", signature=c(x="Function"), function(x, ...){
eval_points <- expand.grid(eval_x, eval_y)
X <- rep(eval_points[,1], times=length(times))
Y <- rep(eval_points[,2], times=length(times))
Z <- x$eval_at(eval_points)
Z <- x$eval(eval_points)

plot_data <- data.frame(X=X, Y=Y,Z=as.vector(Z),
times = rep(times, each=nrow(eval_points)))
44 changes: 18 additions & 26 deletions R/mesh.R
Original file line number Diff line number Diff line change
@@ -3,24 +3,24 @@
.MeshCtr <- setRefClass(
Class = "Mesh", contains = "Domain",
fields = c(
data = "ANY",
m = "integer", # local dim
n = "integer", # embedding dim
cpp_handler = "ANY", ## cpp backend
m = "integer", ## local dim
n = "integer", ## embedding dim
times = "numeric",
deltaT = "numeric" # for parabolic problems
),
methods = c(
get_nodes = function(){
data$nodes()
cpp_handler$nodes()
},
get_elements = function(){
data$elements()+1
cpp_handler$elements()+1
},
get_boundary = function(){
data$boundary()
cpp_handler$boundary()
},
get_neighbors = function(){
data$neighbors()
cpp_handler$neighbors()
},
get_times = function(){
return(times)
@@ -67,7 +67,7 @@ setMethod("Mesh", signature = c(domain="list"),
m <- ncol(domain$elements) - 1
n <- ncol(domain$nodes)
if(m == 2 & n == 2)
.MeshCtr(data=new(Mesh_2D, domain), m=as.integer(m),n=as.integer(n),
.MeshCtr(cpp_handler=new(cpp_2d_domain, domain), m=as.integer(m),n=as.integer(n),
times=vector(mode="double"))
else
stop("wrong input argument provided.")
@@ -79,16 +79,8 @@ setMethod("Mesh", signature = c(domain="list"),
setMethod("Mesh", signature=c(domain="triangulation"),
function(domain){
nodes <- domain$P
#nodes_group <- as.integer(domain$PB)
#edges <- domain$S
#edges_group <- as.integer(domain$SB)
#geometry <- list(nodes = nodes, edges = edges,
# nodes_group = nodes_group, edges_group = edges_group)
elements <- domain$T - 1
boundary <- domain$PB
## 1 exterior boundary edges
## < 0 interior boundary edges (holes)
#boundary[as.vector(domain$E[(domain$EB == 1 | domain$EB < 0),]),] = 1

storage.mode(elements) <- "integer"
storage.mode(nodes) <- "numeric"
@@ -99,7 +91,7 @@ setMethod("Mesh", signature=c(domain="triangulation"),

domain <- list(elements = elements, nodes = nodes, boundary = boundary)
if(m == 2 & n == 2)
.MeshCtr(data=new(Mesh_2D, domain), m=as.integer(m),n=as.integer(n),
.MeshCtr(cpp_handler=new(Mesh_2D, domain), m=as.integer(m),n=as.integer(n),
times=vector(mode="double"),
time_interval = vector(mode="numeric", length = 0), crs = NA_crs_)
else
@@ -128,7 +120,7 @@ setMethod("%X%", signature=c(op1="Mesh", op2="numeric"),

## Mesh - auxiliary methods
unroll_edges_aux <- function(Mesh){
mesh <- Mesh$data
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
@@ -146,18 +138,18 @@ setMethod("unroll_edges", "Mesh", function(Mesh){
plot_mesh_aux <- function(x, ...){
edges <- unroll_edges(x)
plot_ly(...) %>%
add_markers(x = x$data$nodes()[,1],
y = x$data$nodes()[,2],
add_markers(x = x$cpp_handler$nodes()[,1],
y = x$cpp_handler$nodes()[,2],
color = I('black'), size = I(1),
hoverinfo = 'text',
text = paste('</br><b> Coordinates:', round(x$data$nodes()[,1],2),
round(x$data$nodes()[,2],2)),
text = paste('</br><b> Coordinates:', round(x$cpp_handler$nodes()[,1],2),
round(x$cpp_handler$nodes()[,2],2)),
showlegend = T,
visible = T) %>%
add_segments(x = x$data$nodes()[edges[,1],1],
y = x$data$nodes()[edges[,1],2],
xend = x$data$nodes()[edges[,2],1],
yend = x$data$nodes()[edges[,2],2],
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],
color = I('black'), size = I(1),
showlegend = F) %>%
layout(
12 changes: 5 additions & 7 deletions R/modules.R
Original file line number Diff line number Diff line change
@@ -6,11 +6,9 @@
#' @import Matrix
NULL

# #' @exportPattern "^[[:alpha:]]+"


## load required modules
Rcpp::loadModule("PDE_2D_ORDER_1", TRUE)
Rcpp::loadModule("PDE_2D_ORDER_2", TRUE)
Rcpp::loadModule("Mesh_2D", TRUE)

Rcpp::loadModule("cpp_2d_domain", TRUE)
Rcpp::loadModule("cpp_lagrange_basis_2d_fe1", TRUE)
Rcpp::loadModule("cpp_lagrange_basis_2d_fe2", TRUE)
Rcpp::loadModule("cpp_pde_2d_fe1", TRUE)
Rcpp::loadModule("cpp_pde_2d_fe2", TRUE)
Loading

0 comments on commit 3e05f37

Please sign in to comment.