From 5a01604d2782d07ad834d75fe4b5f708881b53ce Mon Sep 17 00:00:00 2001 From: aldoclemente Date: Fri, 10 Nov 2023 18:13:58 +0100 Subject: [PATCH] Update vignettes. Added example concerning mixed boundary conditions. --- tests/mixed_bc.R | 61 ++++++++++++++++++++++++++++++++++++++ tests/parabolic.R | 17 +++++++++-- vignettes/Introduction.Rmd | 20 +++++++++---- 3 files changed, 90 insertions(+), 8 deletions(-) create mode 100644 tests/mixed_bc.R diff --git a/tests/mixed_bc.R b/tests/mixed_bc.R new file mode 100644 index 0000000..5392516 --- /dev/null +++ b/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) diff --git a/tests/parabolic.R b/tests/parabolic.R index a235fbf..38579e7 100644 --- a/tests/parabolic.R +++ b/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) # diff --git a/vignettes/Introduction.Rmd b/vignettes/Introduction.Rmd index 00da4e2..47e8e60 100644 --- a/vignettes/Introduction.Rmd +++ b/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.