diff --git a/articles/Introduction.html b/articles/Introduction.html index 6c26321..2827d89 100644 --- a/articles/Introduction.html +++ b/articles/Introduction.html @@ -88,14 +88,14 @@ the following form:

\[ \begin{cases} --\mu \Delta f + \boldsymbol{b} \cdot \nabla f + a f = u \qquad & in +-\mu \Delta u + \boldsymbol{b} \cdot \nabla u + a u= f \qquad & in \ \Omega \\ boundary \ conditions \ (Dirichlet \ or \ Neumann) \qquad & on \ \partial \Omega, \end{cases} \] where \(\mu \in \mathbb{R}\), \(\boldsymbol{b} \in \mathbb{R}^2\), -\(c \in \mathbb{R}\) are the diffusion +\(a \in \mathbb{R}\) are the diffusion coefficient, the transport coefficient and the reaction coefficient, respectively. This vignette illustrates how to use femR to solve a partial differential equation. In particular, the following @@ -109,56 +109,61 @@

Solving a Poisson problem\(\Omega\):

\[ \begin{cases} --\Delta f = u &in \ \Omega \\ - \quad f = 0 &on \ \partial \Omega +-\Delta u = f &in \ \Omega \\ + \quad u = 0 &on \ \partial \Omega \end{cases} \]

-

where \(u(x,y) = 8\pi^2 \sin( 2\pi +

where \(f(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y)\) is the forcing term and \(\partial \Omega\) is the boundary of \(\Omega\) where we have prescribed homogeneous Dirichelet boundary condition. The exact solution of the previous problem is \(u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y)\). The following figure shows the exact solution of the problem, on the left hand side, and the triangulation of the unit square domain, on the right hand side.

-
-

The following window wraps all the steps needed to solve the +

+

The following window wraps all the steps needed to solve the problem.

-## 1. Defining domain
-data("unit_square", package="femR")
+## 1. Defining domain relying on RTriangle
+pslg <- 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)
 ## 2. Defining the solution of the PDE
-Vh <- FunctionSpace(mesh, fe_order=2) 
-f <- Function(Vh)
+Vh <- FunctionSpace(mesh, fe_order=1) 
+u <- Function(Vh)
 ## 3. Defining the differential operator
-L <- -laplace(f) 
-## 4. Defining the forcing term and the Dirichlet boundary conditions 
+Lu <- -laplace(u) 
+## 4. Defining the forcing term and Dirichlet boundary conditions
 ##    as standard R functions
 # forcing term
-u <- function(points){
+f <- function(points){
     return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
 }
-# Dirichlet BC
+# Dirichlet boundary conditions
 dirichletBC <- function(points){
   return(matrix(0,nrow=nrow(points), ncol=1))
 }
 ## 5. Building the PDE object
-pde <- Pde(L, u, dirichletBC)
-## 6. computing the discrete solution
+pde <- Pde(Lu, f)
+# setting boundary conditions
+pde$set_dirichletBC(dirichletBC)
+## 7. computing the discrete solution
 pde$solve()
-## 7. Plots
-plot(f)
-
-
-fig_h <- contour(f)
+
+## 8. Plots
+plot(u)
+
+
+fig_h <- contour(u)
 subplot(fig_exact %>%hide_colorbar(), fig_h, 
         margin = 0.05) %>%layout(annotations = list(
             list(x=0.155, y=1.05, text = "exact solution",    showarrow = F, 
                  xref='paper', yref='paper'),
             list(x=0.875,  y=1.05, text = "discrete solution", showarrow = F, 
                  xref='paper', yref='paper')))
-
-

The following sections explain more in detail the previous chunk of +

+

The following sections explain more in detail the previous chunk of code.

Defining the domain @@ -179,26 +184,29 @@

Defining the domain
 data("unit_square", package="femR")
-unit_square <- Mesh(unit_square)
-plot(unit_square)

-
-

femR package can read also triangulations provided by -RTriangle package, as the following chunk shows:

+mesh <- Mesh(unit_square) +plot(mesh) +
+

femR package can read also Delaunay triangulations +provided by RTriangle package, as the following chunk +shows:

 library(RTriangle)
 
-## create an object with a concavity 
-pslg <- pslg(P=rbind(c(0, 0), c(0, 1), c(0.5, 0.5), c(1, 1), c(1, 0)),
-          S=rbind(c(1, 2), c(2, 3), c(3, 4), c(4, 5), c(5, 1)))
+## create a Planar Straigh Line Graph object (RTriangle)  
+pslg <- 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)))
 
 ## triangulate it using RTriangle 
-domain <- triangulate(pslg, a=0.01, q=20)
+# a = maximum triangle area 
+# q = minimum triangle angle in degrees
+unit_square <- triangulate(pslg, a=0.00125, q=30)
 
 ## create mesh object
-mesh <- Mesh(domain)
+mesh <- Mesh(unit_square)
 plot(mesh)
-
-

Moreover, femR package provides an utility that reads +

+

Moreover, femR package provides an utility that reads files storing meshes provided by third-party software such as FreeFem++. Indeed, read_mesh(filename) reads a .mesh file from FreeFem++ and returns a named list that can @@ -221,35 +229,40 @@

Defining function space Vh.

 ## solution of the PDE
-Vh <- FunctionSpace(unit_square, fe_order=2)
-f <- Function(Vh)
-

Then, we define the differential operator in its strong +Vh <- FunctionSpace(mesh, fe_order=1) +u <- Function(Vh) +

The following figure shows a basis function of the space +Vh, i.e. the space of functions which are globally +continuous and are polynomials of degree 1 once restricted to each +element of the triangulation.

+
+

Then, we define the differential operator in its strong formulation.

-L <- -laplace(f) # L <- -div(I*grad(f))
-                 # In general:
-                 # L <- -mu*laplace(f) + dot(b,grad(f)) + a*f
-

Note that, according to the well-known identity \(div(\nabla f) = \Delta f\), we could have -written L<--div(I*grad(f)), where I is the +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 -\(Lf= -div(K\nabla f) + \mathbf{b}\cdot \nabla -f + a f\) as follows:

+\(Lu= -div(K\nabla u) + \mathbf{b}\cdot \nabla +u + a u\) as follows:

Although, the interface let the user define a generic second-order -linear differential operator: -div(K*grad(f)) + ..., at +linear differential operator: -div(K*grad(u)) + ..., at this stage the package has been tested only in the isotropic case, i.e K=I, where I is the two-dimensional identity matrix.

@@ -258,16 +271,16 @@

Defining

Building the PDE object

We define the forcing term of the differential problem and the -Dirichlet boundary condition as simple R functions. Such -functions will be provided as parameters to the object which +Dirichlet boundary condition as simple R functions. The +forcing term will be provided as parameter to the object which encapsulates the concept of PDEs.

 ## forcing term
-u <- function(points){
+f <- function(points){
     return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
 }
 
-## Dirichlet BC
+## dirichletBC
 dirichletBC <- function(points){
   return(matrix(0,nrow=nrow(points), ncol=1))
 }
@@ -276,12 +289,14 @@

Building the PDE objectComputing the FE solution

Finally, we build a pde object passing the differential -operator L, as first parameter, the forcing term -u, as second parameter, the Dirichlet boundary condition -dirichletBC, as third parameter:

+operator L, as first parameter and the forcing term +u, as second parameter:

 ## create pde
-pde <- Pde(L, u, dirichletBC)
+pde <- Pde(Lu, f) + +## set Dirichlet boundary conditions +pde$set_dirichletBC(dirichletBC)

We can compute the discrete solution of the Poisson problem calling the solve method:

@@ -297,7 +312,7 @@ 
Computing the FE solutionu_ex <- as.matrix(exact_solution(pde$get_dofs_coordinates()))
 error.L2 <- sqrt(sum(pde$get_mass()%*%(u_ex-pde$solution())^2))
 cat("L2 error = ", error.L2, "\n")
-#> L2 error =  0.0004136347
+#> L2 error = 0.001903506

Graphics @@ -306,20 +321,20 @@

Graphics -
-
+plot(u)
+
+
 # contour 
-contour(f)
-
-

Note that, the base::plot function and +contour(u)

+
+

Note that, the base::plot function and graphics::contourfunction have both been overloaded and return a plotly object that can be modified as one wishes. 

-plot(f, colorscale="Jet") %>% layout(scene=list(aspectmode="cube"))
-
- +plot(u, colorscale="Jet") %>% layout(scene=list(aspectmode="cube")) +
+