Update vignettes. Added example concerning mixed boundary conditions.

tests/mixed_bc.R

@@ -0,0 +1,61 @@
+library(RTriangle)
+library(femR)
+
+p <- pslg(P=rbind(c(0, 0),  # first node
+                  c(1, 0),  # ...   node
+                  c(1, 1),  # ...   node
+                  c(0, 1)), # ...   node
+          S=rbind(c(1, 2),  # first physical edge 
+                  c(2, 3),  # ...   physical edge
+                  c(3, 4),  # ...   physical edge
+                  c(4,1)),  # ...   physical edge
+    SB = cbind(c(1,2,3,4))) # use SB to specify "physical boundaries"
+                            # 1 bottom, 2 right, 3 up, 4 left 
+
+plot(p)
+unit_square <- triangulate(p, a = 0.00125, q=30)
+# select nodes which do not belong to 1 and 4 physical edges
+mask <- unit_square$PB != 1 & unit_square$PB != 4  
+
+# Homogeneous Neumann BC will be imposed on nodes belonging to 2 and 3
+unit_square$PB[mask,] = 0 
+
+# Dirichlet BC will be imposed on nodes belonging to 1 and 4
+unit_square$PB[!mask,] = 1
+
+plot(unit_square)
+points(unit_square$P[unit_square$PB==1,],pch=16,col="red") # dirichlet
+
+# create list to be passed to Mesh()
+domain <- list( elements = unit_square$T, 
+                nodes = unit_square$P, boundary = unit_square$PB)
+
+## 1. building mesh
+mesh <- Mesh(domain)
+plot(mesh)
+
+## 2. Defining the solution of the PDE
+Vh <- FunctionSpace(mesh, fe_order=1) 
+u <- Function(Vh)
+## 3. Defining the differential operator
+Lu <- -laplace(u) 
+## 4. Defining the forcing term and Dirichlet boundary conditions
+##    as standard R functions
+# forcing term
+f <- function(points){
+  return(4*((points[,1]^2+points[,2]^2) * sin(points[,1]^2 + points[,2]^2 - 1) 
+            - cos(points[,1]^2 + points[,2]^2 - 1))) 
+}
+# Dirichlet boundary conditions 
+dirichletBC <- 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)
+## 7. computing the discrete solution
+pde$solve()
+
+## 8. Plots
+plot(u)

tests/parabolic.R

@@ -94,7 +94,18 @@ for( t in 1:length(times)){
 # 
 # ## plot solution 
 # options(warn=-1)
-# plot(f) %>% hide_colorbar()
-#  
-# contour(f)
+plot(u) %>% hide_colorbar()
+
+# scene = list(camera = list(eye = list(x = 0, y = 0, z = 1.25)),
+#              xaxis = list(autorange="reversed"),
+#              yaxis = list(autorange="reversed"))
+# 
+# plot_tmp <- plot(u) %>% 
+#   layout(dragmode=FALSE, scene=scene,aspectmode="cube") %>%
+#   hide_colorbar() %>% 
+#   config(displayModeBar =FALSE)
+# 
+# plot_tmp
+# TO DO  
+# contour(u)
 # 

vignettes/Introduction.Rmd

@@ -83,7 +83,7 @@ The following window wraps all the steps needed to solve the problem.
 
 ```{r code block}
 ## 1. Defining domain relying on RTriangle
-pslg <- pslg(P=rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1)),
+p <- pslg(P=rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1)),
           S=rbind(c(1, 2), c(2, 3), c(3, 4), c(4,1)))
 unit_square <- triangulate(p, a = 0.00125, q=30)
 mesh <- Mesh(unit_square)
@@ -198,11 +198,21 @@ coeff_aux[5] = 1
 
 hat_function$set_coeff(as.matrix(coeff_aux))
 
+# scene = list(camera = list(eye = list(x = 0, y = 0, z = 1.25)),
+#              yaxis = list(autorange="reversed"),
+#              projection=list(type="orthographic"))
+# 
+# plot_hat <- plot(hat_function) %>% layout(scene=scene,aspectmode="cube",
+#                                           dragmode=FALSE) %>%hide_colorbar() %>% config(displayModeBar =FALSE, scrollZoom = FALSE)
+
 scene = list(camera = list(eye = list(x = 0, y = 0, z = 1.25)),
-             yaxis = list(autorange="reversed"),
-             projection=list(type="orthographic"))
+             xaxis = list(autorange="reversed"),
+             yaxis = list(autorange="reversed"))
 
-plot_hat <- plot(hat_function) %>% layout(scene=scene,aspectmode="cube") %>%hide_colorbar()
+plot_hat <- plot(hat_function) %>% 
+            layout(dragmode=FALSE, scene=scene,aspectmode="cube") %>%
+            hide_colorbar() %>% 
+            config(displayModeBar =FALSE)
 
 plot_hat
 ```
@@ -213,7 +223,7 @@ Lu <- -laplace(u) # L <- -div(I*grad(u))
                   # In general:
                   # L <- -mu*laplace(u) + dot(b,grad(u)) + a*u
 ```
-Note that, according to the well-known identity $div(\nabla u) = \Delta f$, we could have written `L<--div(I*grad(u))`, where `I` is the 2-dimensional identity matrix. Indeed, `femR` package let you define every term of a generic second-order linear differential operator $Lu= -div(K\nabla u) + \mathbf{b}\cdot \nabla u + a u$ as follows:
+Note that, according to the well-known identity $div(\nabla u) = \Delta u$, we could have written `L<--div(I*grad(u))`, where `I` is the 2-dimensional identity matrix. Indeed, `femR` package let you define every term of a generic second-order linear differential operator $Lu= -div(K\nabla u) + \mathbf{b}\cdot \nabla u + a u$ as follows:
 
   * Diffusion: `-div(K*grad(u))` or `-mu*laplacian(u)`, where either `K` or `mu` represents the diffusion matrix in an anisotropic problem or the diffusion coefficient in a isotropic problem, implements the first term of the differential operator according to the problem at hand.