From d9380ac02fcf7744695e2047586dcbee1a1db30e Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Thu, 30 Nov 2023 16:23:27 +0000 Subject: [PATCH] Updated DESCRIPTION. Minor corrections. --- DESCRIPTION | 3 +- R/pde.R | 15 +-- tests/adv.diff.K.R | 6 +- tests/adv.diff.R | 6 +- tests/functional_space_interface.R | 2 +- tests/mixed_bc.R | 8 +- tests/parabolic.R | 10 +- tests/pde.2.R | 7 +- vignettes/Introduction.Rmd | 12 ++- vignettes/Lake_Como.Rmd | 162 +++++++++++++++++++++++++++++ vignettes/Parabolic.Rmd | 9 +- 11 files changed, 205 insertions(+), 35 deletions(-) create mode 100644 vignettes/Lake_Como.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index 698eed3..0450e49 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -28,5 +28,6 @@ Encoding: UTF-8 Suggests: knitr, rmarkdown, - RTriangle + RTriangle, + sf VignetteBuilder: knitr diff --git a/R/pde.R b/R/pde.R index 2a1b1a7..10ec482 100644 --- a/R/pde.R +++ b/R/pde.R @@ -26,28 +26,31 @@ if(!is_init) pde_$init() pde_$solve() }, - solution = function(){ + get_solution = function(){ pde_$solution() }, get_dofs_coordinates = function(){ pde_$get_dofs_coordinates() }, - set_dirichletBC = function(dirichletBC){ + set_boundary_condition = function(fun, type="dirichlet", on=NULL){ + if(!any(type == c("dirichlet", "Dirichlet", "d"))) + stop("Only Dirichlet boundary condtions allowed.") + if(!is_parabolic){ - dirichletBC_ <- as.matrix(dirichletBC(pde_$get_dofs_coordinates())) + dirichletBC_ <- as.matrix(fun(pde_$get_dofs_coordinates())) }else{ - dirichletBC_ <- dirichletBC(pde_$get_dofs_coordinates(),times) + dirichletBC_ <- fun(pde_$get_dofs_coordinates(),times) pde_$set_dirichlet_bc(dirichletBC_) } is_dirichletBC_set <<- TRUE pde_$set_dirichlet_bc(dirichletBC_) }, - set_initialCondition = function(initialCondtion){ + set_initialCondition = function(fun){ if(!is_parabolic) stop("Cannot set initial condition for elliptic problem.") is_initialCondition_set <<- TRUE - pde_$set_initial_condition(initialCondition(pde_$get_dofs_coordinates())) + pde_$set_initial_condition(fun(pde_$get_dofs_coordinates())) }, get_mass = function(){ pde_$get_mass() diff --git a/tests/adv.diff.K.R b/tests/adv.diff.K.R index 1b92fdd..ee72fe1 100644 --- a/tests/adv.diff.K.R +++ b/tests/adv.diff.K.R @@ -40,20 +40,20 @@ u <- function(points){ return(gamma_ * sin(pi * points[,2])) } ## Dirichlet BC -dirichletBC <- function(points){ +g <- function(points){ return(matrix(0,nrow=nrow(points), ncol=1)) } ## create pde pde <- Pde(L, u) -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(g) ## solve problem 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 - pde$solution())^2)) +error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2)) cat("L2 error = ", error.L2, "\n") ## perform evaluation at single point diff --git a/tests/adv.diff.R b/tests/adv.diff.R index 597f929..717243b 100644 --- a/tests/adv.diff.R +++ b/tests/adv.diff.R @@ -36,18 +36,18 @@ u <- function(points){ return(gamma_ * sin(pi * points[,2])) } ## Dirichlet BC -dirichletBC <- function(points){ +g <- function(points){ return(matrix(0,nrow=nrow(points), ncol=1)) } ## create pde pde <- Pde(L, u) -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(g) ## solve problem 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 - pde$solution())^2)) +error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2)) cat("L2 error = ", error.L2, "\n") ## perform evaluation at single point diff --git a/tests/functional_space_interface.R b/tests/functional_space_interface.R index adaa7d1..062adca 100644 --- a/tests/functional_space_interface.R +++ b/tests/functional_space_interface.R @@ -22,7 +22,7 @@ dirichletBC <- function(points){ } ## 5. Building the PDE object pde <- Pde(L, u) -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(dirichletBC) ## 6. computing the discrete solution pde$solve() ## 7. Plots diff --git a/tests/mixed_bc.R b/tests/mixed_bc.R index 9ab23fe..3e8e730 100644 --- a/tests/mixed_bc.R +++ b/tests/mixed_bc.R @@ -63,7 +63,7 @@ plot(st_as_sfc(mesh)) points(mesh$get_nodes()[dirichlet_id,], pch=16) # Dirichlet boundary conditions -dirichletBC <- function(points, times){ +g <- function(points, times){ res <- matrix(0, nrow=nrow(points), ncol=length(times)) for( t in 1:length(times)){ res[dirichlet_id,t] <- rep(Q*exp(-times[t]), times=length(dirichlet_id)) @@ -71,7 +71,7 @@ dirichletBC <- function(points, times){ return(res) } -initialCondition <- function(points){ +u0 <- function(points){ res <- matrix(0, nrow=nrow(points)) res[dirichlet_id,] <- rep(Q, times=length(dirichlet_id)) return(res) @@ -79,10 +79,10 @@ initialCondition <- function(points){ ## 5. Building the PDE object pde <- Pde(Lu, f) # setting initial conditions -pde$set_initialCondition(initialCondition) +pde$set_initialCondition(u0) # setting boundary conditions -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(g) ## 7. computing the discrete solution pde$solve() diff --git a/tests/parabolic.R b/tests/parabolic.R index ba60a2e..9a809a3 100644 --- a/tests/parabolic.R +++ b/tests/parabolic.R @@ -44,19 +44,19 @@ f <- function(points,times){ return(res) } ## Dirichlet BC -dirichletBC <- function(points, times){ +g <- function(points, times){ return(matrix(0, nrow=nrow(points), ncol=length(times))) } # initial condition -initialCondition <- function(points){ +u0 <- function(points){ return(sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2])) } ## create pde pde <- Pde(L, f) -pde$set_dirichletBC(dirichletBC) -pde$set_initialCondition(initialCondition) +pde$set_boundary_condition(g) +pde$set_initialCondition(u0) pde$solve() ## compute L2 norm of the error @@ -64,7 +64,7 @@ u_ex <- exact_solution(pde$get_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] - pde$solution()[,t])^2)) + error.L2[t] <- sqrt(sum(pde$get_mass() %*% (u_ex[,t] - pde$get_solution()[,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 6c617b5..d08c0d9 100644 --- a/tests/pde.2.R +++ b/tests/pde.2.R @@ -23,19 +23,20 @@ u <- function(points){ return(8.*pi^2* sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]) ) } ## Dirichlet BC -dirichletBC <- function(points){ +g <- function(points){ return(rep(0, times=nrow(points))) } ## Pde constructor pde <- Pde(L,u) -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(fun=g, + type = "dirichlet") ## solve problem 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 - pde$solution())^2)) +error.L2 <- sqrt(sum(pde$get_mass() %*% (u_ex - pde$get_solution())^2)) cat("L2 error = ", error.L2, "\n") ## perform evaluation at single point diff --git a/vignettes/Introduction.Rmd b/vignettes/Introduction.Rmd index 2975f5d..3c497dd 100644 --- a/vignettes/Introduction.Rmd +++ b/vignettes/Introduction.Rmd @@ -99,13 +99,14 @@ f <- function(points){ return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) } # Dirichlet boundary conditions -dirichletBC <- function(points){ +g <- function(points){ return(matrix(0,nrow=nrow(points), ncol=1)) } ## 5. Building the PDE object pde <- Pde(Lu, f) # setting boundary conditions -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(type="dirichlet", + fun=g) ## 7. computing the discrete solution pde$solve() @@ -242,7 +243,7 @@ f <- function(points){ } ## dirichletBC -dirichletBC <- function(points){ +g <- function(points){ return(matrix(0,nrow=nrow(points), ncol=1)) } ``` @@ -253,7 +254,8 @@ Finally, we build a `pde` object passing the differential operator `L`, as first pde <- Pde(Lu, f) ## set Dirichlet boundary conditions -pde$set_dirichletBC(dirichletBC) +pde$set_boundary_condition(type="dirichlet", + fun=g) ``` We can compute the discrete solution of the Poisson problem calling the `solve` method: ```{r Solving the problem} @@ -267,7 +269,7 @@ 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-pde$solution())^2)) +error.L2 <- sqrt(sum(pde$get_mass()%*%(u_ex-pde$get_solution())^2)) cat("L2 error = ", error.L2, "\n") ``` diff --git a/vignettes/Lake_Como.Rmd b/vignettes/Lake_Como.Rmd new file mode 100644 index 0000000..713cf46 --- /dev/null +++ b/vignettes/Lake_Como.Rmd @@ -0,0 +1,162 @@ +--- +title: 'Concentration of pollutant spreading in Lake Como' +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Concentration of pollutant spreading in Lake Como} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE,fig.width= 7., + fig.height= 5., + fig.align= "center", + collapse = TRUE, + comment = "#>") + +options(warn = -1) +``` + +```{r} +library(sf, quietly = T) +library(RTriangle, quietly = T) +library(mapview, quietly = T) +library(femR, quietly = T) + +lake_shp <- sf::read_sf("/home/aldo/Documents/como_lake_boundaries/deims_sites_boundariesPolygon.shp") +st_crs(lake_shp) <- 4326 +mapview(lake_shp, col.regions = "white", lwd = 2) +``` + +```{r} +lake_boundary <- st_geometry(lake_shp) +nodes <- unique(st_coordinates(lake_boundary)) +edges <- cbind( 1:nrow(nodes),c(2:nrow(nodes),1)) +pslg <- pslg(P=nodes[,1:2], S = edges) +triangulation <- triangulate(pslg, a = 1e-5, q=20) +mesh <- Mesh(triangulation) + +poly_list <- list() +for(e in 1:nrow(mesh$get_elements())){ + element_sf <- st_polygon(list(rbind(mesh$get_nodes()[mesh$get_elements()[e,]+1,], + mesh$get_nodes()[mesh$get_elements()[e,1]+1,]))) + poly_list[[e]] <- element_sf + +} +poly_list <- st_sfc(poly_list, crs=4326) +plot(poly_list) +mesh_sf <- data.frame(elem_id=rep(NA,times=nrow(mesh$get_elements()))) +mesh_sf <- st_sf(mesh_sf, geometry = poly_list) + +map.type <- "Esri.WorldTopoMap" +mapview(mesh_sf, col.region="white",alpha.regions=0.2,legend=F,map.type=map.type) + + mapview(lake_shp, col.regions = "white", alpha.regions=0.2, lwd = 2, legend=F) +``` + + +```{r} +set.seed(0) +num_stations <- 20 +boundary_points <- which(mesh$boundary() == 1) +boundary_id <- sample(x = 1:length(boundary_points),size = num_stations) +boundary_id + +sources <- list() +for(i in 1:length(boundary_id)) + sources[[i]] <- st_point(x=mesh$get_nodes()[boundary_id[i],]) + +sources <- st_sfc(sources, crs=4326) + +``` +```{r} +initialCondition <- function(points){ + initialCoeff +} + +dirichletBC <- function(points, times){ + res <- matrix(nrow=nrow(points), ncol=length(times)) + res[,1] <- initialCoeff + for(t in 2:length(times)){ + res[,t] <- initialCoeff*(1-times[t]) + } + return(res) +} + +f <- function(points, times){ + res <- matrix(nrow=nrow(points), ncol=length(times)) + for(t in 2:length(times)){ + res[,t] <- exp(sqrt(points[,1]-x0)^2 + (sqrt(points[,1]-x0)^2)^2)*exp(-times[t]) + } + return(res) +} + + +# PDEs Solution +Vh <- FunctionSpace(mesh) +u <- Function(Vh) +Lu <- dt(u) - laplace(u) + dot(c(1,0), grad(u)) + +# PDE object +pde <- Pde(Lu,f) +pde$set_dirichletBC(dirichletBC) +pde$set_initialCondition(initialCondition) +pde$solve() +``` + +```{r} +set.seed(0) +num_stations <- 20 +boundary_points <- which(mesh$boundary() == 1) +boundary_id <- sample(x = 1:length(boundary_points),size = num_stations) +stations <- list() +for(i in 1:length(boundary_id)) + stations[[i]] <- st_point(x=mesh$get_nodes()[boundary_id[i],]) + +stations <- st_sfc(stations, crs=4326) +station_name <- as.character(1:num_stations) +stations <- st_as_sf(data.frame(station_name=station_name), + geometry = stations) + +times <-seq(from=0, to=1, by=0.125) +mesh <- mesh %X% times +observations <- matrix(round(rnorm(num_stations*length(times),mean = 0.13,sd=0.01), digits = 3), + nrow=num_stations, ncol=length(times)) + +mapview(stations, legend=F, cex=5, col.region="black") +observations <- matrix(round(rnorm(num_stations,mean = 0.13,sd=0.01), digits = 2), + nrow=num_stations, ncol=1) +``` + +```{r} +initialCoeff <- round(rnorm(nrow(mesh$get_nodes()), + mean=0.13, sd=0.01), digits = 2) + +# dummmy function +initialCondition <- function(points){ + initialCoeff +} + +dirichletBC <- function(points, times){ + res <- matrix(nrow=nrow(points), ncol=length(times)) + res[,1] <- initialCoeff + for(t in 2:length(times)){ + res[,t] <- initialCoeff*(1-times[t]) + } + return(res) +} + +f <- function(points, times){ + matrix(0, nrow=nrow(points), ncol=length(times)) +} + +# PDEs Solution +Vh <- FunctionSpace(mesh) +u <- Function(Vh) +Lu <- dt(u) - laplace(u) + dot(c(1,0), grad(u)) + +# PDE object +pde <- Pde(Lu,f) +pde$set_dirichletBC(dirichletBC) +pde$set_initialCondition(initialCondition) +pde$solve() +``` diff --git a/vignettes/Parabolic.Rmd b/vignettes/Parabolic.Rmd index ebb742f..86600d2 100644 --- a/vignettes/Parabolic.Rmd +++ b/vignettes/Parabolic.Rmd @@ -106,12 +106,12 @@ f <- function(points,times){ return(res) } ## Dirichlet BC -dirichletBC <- function(points, times){ +g <- function(points, times){ return(matrix(0, nrow=nrow(points), ncol=length(times))) } ## initial condition -initialCondition <- function(points){ +u0 <- function(points){ return(sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2])) } @@ -119,8 +119,9 @@ initialCondition <- function(points){ pde <- Pde(Lu, f) # set initial condition and dirichlet BC -pde$set_initialCondition(initialCondition) -pde$set_dirichletBC(dirichletBC) +pde$set_initialCondition(u0) +pde$set_boundary_condition(type="dirichlet", + fun=g) ## solve problem pde$solve()