Introduction • femR

femR is an R package that solves PDEs on R relying on the Finite Element Method (FEM). Let \(\Omega \subset \mathbb{R}^2\) be the domain of our interest, femR package solves advection-reaction-diffusion problem of the following form:

\[ \begin{cases} -\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\), \(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 section explains step by step the solution of a simple Poisson problem.

Solving a Poisson problem

Let \(\Omega\) be the unit square \([0,1]^2\). Consider the following problem defined over \(\Omega\):

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

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 problem.

## 1. Defining domain relying on RTriangle
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)
## 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(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
}
# Dirichlet boundary conditions
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_boundary_condition(type="dirichlet",
                           boundary_condition=g)
## 7. computing the discrete solution
pde$solve()

## 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 code.

Defining the domain

First, we create a proper domain object exploiting the function Mesh which takes a R named list as parameter. Such named list should contain geometrical information of the domain:

nodes, a 2-columns matrix containing the coordinates of the nodes of the mesh.
elements, a 3-columns matrix that for each row \(i\) contains the indexes of the mesh nodes that are vertices of the \(i\)-th element.
boundary, a 1-column matrix that indicates the nodes on the boundary of the domain.

femR contains the unit_square named list that can be used to perform some tests.

data("unit_square", package="femR")
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 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 
# a = maximum triangle area 
# q = minimum triangle angle in degrees
unit_square <- triangulate(pslg, a=0.00125, q=30)

## create mesh object
mesh <- Mesh(unit_square)
plot(mesh)

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 be passed to Mesh function. The following chunk, that is not executed, shows how read_mesh works.

filename <- "path/to/domain.mesh" # domain.mesh stored by the FreeFem++ 
                                  # utility savemesh() 
domain <- read_mesh(filename)
mesh <- Mesh(domain)

Defining the solution of the PDE and the differential operator

First, we initialize a FunctionSpace object passing the mesh that we previously build and the finite element order. Note that, the package provide first order and second order finite elements. First order finite element is the default parameter. Then, we define a Function, solution of the PDEs, belonging to the function space Vh.

## solution of the PDE
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.

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 u\), we could have written L<--div(grad(u)) or 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.
Advection: dot(b,grad(u)), where b is a two-dimensional vector, implements the second term of differential operator.
Reaction: a*u, where a is a scalar, implements the last term of the differential operator.

Although, the interface let the user define a generic second-order 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.

Building the PDE object

We define the forcing term of the differential problem and the 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
f <- function(points){
    return(8.*pi^2* sin(2.*pi*points[,1])*sin(2.*pi*points[,2])) 
}

## dirichletBC
g <- function(points){
  return(matrix(0,nrow=nrow(points), ncol=1))
}

Computing the FE solution

Finally, we build a pde object passing the differential operator L, as first parameter and the forcing term f, as second parameter:

## create pde
pde <- Pde(Lu, f)

## set Dirichlet boundary conditions
pde$set_boundary_condition(type="dirichlet",
                           boundary_condition=g)

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

## solve problem
pde$solve()

Since, we are dealing with a simple differential problem and we know the analytic form of the exact solution, we can compute the \(L^2\) norm of the error exploiting on the mass matrix provided by the pde object:

exact_solution <- function(points){
    return( sin(2.*pi*points[,1])*sin(2.*pi*points[,2]))
}
u_ex <- as.matrix(exact_solution(pde$dofs_coordinates()))
error.L2 <- sqrt(sum(pde$mass()%*%(u_ex-u$coefficients())^2))
cat("L2 error = ", error.L2, "\n")
#> L2 error =  0.001903506

Graphics

The R package plotly is exploited to plot the computed solution.

## plot solution 
plot(u)

# contour 
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(u, colorscale="Jet") %>% layout(scene=list(aspectmode="cube"))

Solving a parabolic problem

Let \(\Omega\) be the unit square \([0,1]^2\). Consider the following problem defined over \(\Omega \times [0,T]\):

\[ \begin{cases} \frac{\partial u}{\partial t}-\Delta u = f &in \ \Omega \times [0,T] \\ \quad u = u_0 &in \ \Omega, \ t=0 \\ \quad u = 0 &on \ \partial \Omega,\ t>0 \end{cases} \]

where \(u_0 = \sin( 2\pi x)\sin(2\pi y)\) is the initial condition, \(f(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y) e^{-t}\) is the forcing term and \(\partial \Omega\) is the boundary of \(\Omega\) where we have prescribed homogeneous Dirichelet boundary condition for all \(t\) grater than \(0\). The exact solution of the previous problem is \(u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y) e^{-t}\). The following window wraps all the steps needed to solve the problem.

## spatio-temporal domain
nodes <- rbind(c(0,0), c(1,0), c(1,1), c(0,1))
edges <- cbind(1:nrow(nodes), c(2:nrow(nodes),1))
domain <- Domain(list(nodes=nodes, edges=edges))
domain <- domain %X% c(0,1)

## build mesh
mesh <- build_mesh(domain, maximum_area = 0.00125, minimum_angle = 20)
# time step
dT = 1e-2
mesh$set_time_step(dT)
## create Functional Space
Vh <- FunctionSpace(mesh)
## define differential operator in its strong formulation
u <- Function(Vh)
## define differential operator in its strong formulation
Lu <- dt(u) - laplace(u)
## forcing term
f <- function(points,times){
  res <- matrix(0, nrow=nrow(points), ncol=length(times))
  for( t in 1:length(times)){
    res[,t] = (8*pi^2 -1) * sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]) * exp(-times[t])
  }
  return(res) 
}
## initial condition 
u0 <- function(points){
 return(sin( 2.* pi * points[,1]) * sin(2.*pi* points[,2]))
}
## create pde
pde <- Pde(Lu, f)
# set initial condition and dirichlet BC
pde$set_initial_condition(u0)
pde$set_boundary_condition(boundary_condition= 0.0, type="dirichlet")
## solve problem
pde$solve()
## plot solution
image(u)