Updated DESCRIPTION. Minor corrections.

DESCRIPTION

@@ -28,5 +28,6 @@ Encoding: UTF-8
 Suggests: 
     knitr,
     rmarkdown,
-    RTriangle
+    RTriangle,
+    sf
 VignetteBuilder: knitr

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()

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

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

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

tests/mixed_bc.R

@@ -63,26 +63,26 @@ 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))
   }
   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)
 }
 ## 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()
 

tests/parabolic.R

@@ -44,27 +44,27 @@ 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
 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"))
 }
 # ------------------------------------------------------------------------------

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

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")
 ```
 

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()
+```