FEM Introduction¶

Let

$\Omega$ be the domain, we aim at solving the Laplace equation (in 1D): \[-\Delta u(x) = -\frac{\partial^2}{\partial x^2}u(x) =0\]

Subject to boundary conditions

$g$ \[u|_{\partial\Omega} = g,\]

where $u|_{\partial\Omega}$ is the boundary of $\Omega$.

In 1D, the boundaries are the 2 endpoints $\Omega_0$ and $\Omega_1$, so the boundary conditions are:

\[ u(\Omega_0) = u_0\qquad\mathrm{and}\qquad u(\Omega_1) = u_1. \]

Weak Form¶

Instread of solving \[-\Delta u =0\]

we multiply it by a test function $v$ and integrate over the domain

\[-\int_\Omega\Delta u v =0, \qquad \forall v.\]

This equation is called weak form of the original PDE.

If it holds for any $v$, then $u$ is also a solution of the original PDE (strong form).

We use integration by parts to simplify the weak form

\[-\int_\Omega\Delta u v = \]

\[ =-\Bigg(\bigg[\nabla u\, v \bigg]_{\Omega_0}^{\Omega_1}-\int_\Omega\nabla u \cdot \nabla v\Bigg) = \]

\[=\int_\Omega\nabla u \cdot \nabla v = 0, \qquad \forall v\]

Discretization¶

We express the unknown function

$u$ in terms of a dicrete basis $\phi_i$, $i=0,\dots,n$.

\[u(x)=\sum_{i=0}^n u_i \phi_i(x)\]

We insert this definition in the weak form

\[ \int_\Omega\nabla \sum_{i=0}^n(u_i \phi_i) \cdot \nabla v = \sum_{i=0}^n u_i \int_\Omega\nabla \phi_i \cdot \nabla v = 0, \qquad \forall v\]

We use the same bases $\phi_j$ for $v$ and plug it in the weak form

\[\sum_{i=0}^n u_i \int_\Omega\nabla \phi_i \cdot \nabla \phi_j = 0, \qquad \forall j=0,\dots,n\]

This expression can be rewritten in matrix form \[ L \vec{u} =0,\qquad\mathrm{where}\qquad L_{i,j} = \int_\Omega\nabla \phi_i \cdot \nabla \phi_j \] and $\vec{u}$ is the vector containing the $u_i$.

Example¶

import numpy as np
import scipy.sparse as spr
from scipy.sparse.linalg import spsolve

import plotly.offline as plotly
import plotly.graph_objects as go
import plotly.figure_factory as ff

# For static visualization
import plotly.io as pio
pio.renderers.default = "svg"

# For interactive visualization
#pio.renderers.defaut = "notebook"

The domain $\Omega = [0, 1]$, which we discretize with $n_{el}$ segments (or elements) $s_i$

#domain
omega = np.array([0, 1])

#number of bases and elements
n_el = 10
n_bases = n_el + 1

#segments
s = np.linspace(omega[0], omega[1], num=n_el+1)
# s = np.cumsum(np.random.rand(n_elements+1))
s = (s-s[0])/(s[-1]-s[0])

#plot
go.Figure(data=[go.Scatter(x=s, y=np.zeros(s.shape), mode='lines+markers')])

We can split the integrals in the weak form

\[\sum_{i=0}^n u_i \int_\Omega\nabla \phi_i \cdot \nabla \phi_j = 0, \qquad \forall j=0,\dots,n\] from the whole domain $\Omega$ to each individual element $s_e$.

\[ %\sum_{i=0}^n u_i \int_\Omega\nabla \phi_i \cdot \nabla \phi_j = \sum_{i=0}^n u_i \sum_{e=0}^{n_{el}} \int_{s_e}\nabla \phi_i \cdot \nabla \phi_j = 0, \qquad \forall j=0,\dots,n\]

We want to have interpolatory bases, so let’s have value 1 over the node and 0 everywhere else

phis = []

for i in range(n_bases):
    phi = np.zeros(s.shape)
    phi[i] = 1
    phis.append(go.Scatter(x=s, y=phi, mode='lines+markers', name="$\phi_{{{}}}$".format(i)))

go.Figure(data=phis)

Local bases¶

We want to localize the bases, since every piece looks similar.

For simplicity we define the reference element $\hat s= [0, 1]$, a segment of unit length.

On each element we have only 2 non-zero local bases. We define their “pieces” on $\hat s$.

#definition of bases
def hat_phi0(x):
    return 1-x
def hat_phi1(x):
    return x

We can now plot the two bases

x = np.linspace(0, 1)
go.Figure(data=[
    go.Scatter(x=x, y=hat_phi0(x), mode='lines', name="$\hat\phi_0$"),
    go.Scatter(x=x, y=hat_phi1(x), mode='lines', name="$\hat\phi_1$")
])

We now need to map the reference element $\hat s$ to each individual segment $s_e$.

This is called geometric mapping and maps the local segment $\hat s$ to each global segment $s_e$:

\[g_e(\hat x) = s_{e,0} + \hat x (s_{e,1} - s_{e,0})\]

where $s_{e,0}$ and $s_{e,1}$ are the start and end points of $s_e$.

The integrals are over the global elements $s_e$.

\[ \sum_{i=0}^n u_i \sum_{e=0}^{n_{el}} \int_{s_e}\nabla \phi_i(x) \cdot \nabla \phi_j(x)\, \mathrm{d} x =0 , \qquad \forall j \]

We perform a change of variable in the integral to integrate over the reference element $\hat s$.

\[ \sum_{i=0}^n u_i \sum_{e=0}^{n_{el}} \int_{\hat s}\nabla \phi_i(g_e(\hat x)) \cdot \nabla \phi_j(g_e(\hat x))\, g_e^{\prime}(\hat x)\, \mathrm{d}\hat x =0 , \qquad \forall j \]

Note that for 2D and 3D the jacobian of the geometric mapping appears, as explained here.

We can do more!

We note:

\[ \phi_i(g_e(\hat x)) = \hat \phi_i(\hat x) \]

\[ g_e^{\prime}(x)=s_{e, 1} - s_{e, 0} \]

and \[ \nabla_x \phi_i(g_e(\hat x)) = \frac{\nabla_{\hat x} \hat \phi_i(\hat x)}{s_{e, 1} - s_{e, 0}} \]

The weak form

\[ \sum_{i=0}^n u_i \sum_{e=0}^{n_{el}} \int_{\hat s}\nabla \phi_i(g_e(\hat x)) \cdot \nabla \phi_j(g_e(\hat x))\, g_e^{\prime}(\hat x)\, \mathrm{d}\hat x =0 , \qquad \forall j \]

simplifies to

\[ \sum_{i=0}^n u_i \sum_{e=0}^{n_{el}} \int_{\hat s}\frac{\nabla \hat\phi_i(\hat x)}{s_{e, 1} - s_{e, 0}} \cdot \frac{\nabla \hat\phi_j(\hat x)}{s_{e, 1} - s_{e, 0}}\,(s_{e, 1} - s_{e, 0}) \, \mathrm{d}\hat x = \]

\[ \sum_{i=0}^n u_i \sum_{e=0}^{n_{el}} \int_{\hat s}\frac{\nabla \hat\phi_i(\hat x) \cdot \nabla \hat\phi_j(\hat x)}{s_{e, 1} - s_{e, 0}} \, \mathrm{d}\hat x = 0 , \qquad \forall j \]

This localization forces us to keep track of the mapping between the 2 local nodes and their respective global indices.

This mapping is called local to global.

In other words, the local to global mapping $g_e^i$ maps the local indices $i=0,1$ of element $e$ to its corresponding global indices.

Note that most of the terms are zero since only 2 bases are not zero.

Using this note and the local to global mapping we can further simplify

\[ \sum_{i=0}^n u_i \sum_{e=0}^{n_{el} } \int_{\hat s} \frac{\nabla \hat\phi_i \cdot \nabla \hat\phi_j}{s_{e, 1} - s_{e, 0} } = \sum_{e=0}^{n_{el} }\sum_{i=0}^1 u_{g_e^i} \int_{\hat s} \frac{\nabla \hat\phi_i \cdot \nabla \hat\phi_j}{s_{e, 1} - s_{e, 0}} , \qquad \forall j=0,1 \]

We need the gradients of the local bases

\[ \hat \phi_0 = 1-x \]

and

\[ \hat \phi_1 = x \]

def grad_hat_phi0(x):
    return -np.ones(x.shape)
def grad_hat_phi1(x):
    return np.ones(x.shape)

Basis construction¶

We now construct an array of elements, one for each

$s_e$ .

Each element $e$ contains: - the number of non-zero bases (always 2), - the 2 functions and their 2 gradients, - the local to global mapping $g_e^i$ , - the geometric mapping and its gradient.

Note that in this case the local to global mapping is trivial:

\[g_e^i = e+i.\]

elements = []
for e in range(n_el):
    el = {}

    el["n_bases"] = 2

    #2 bases
    el["phi"] = [hat_phi0, hat_phi1]
    el["grad_phi"] = [grad_hat_phi0, grad_hat_phi1]

    #local to global mapping
    el["loc_2_glob"] = [e, e+1]

    #geometric mapping
    el["gmapping"] = lambda x, e=e : s[e] + x*(s[e+1]-s[e])
    el["grad_gmapping"] = lambda x : (s[e+1]-s[e])

    elements.append(el)

We define a function to interpolate the vector $\vec{u}$ using the local to global, geometric mapping, and local bases to interpolate the data.

We first define a vector of $\hat x$ values (xhat = np.linspace(0, 1)) where we want to evaluate the local bases,

then we iterate over all elements and compute

\[ \sum_{i=0}^1 u_{g_e^i} \hat \phi_i(\hat x). \]

We also map $\hat x$ to its global position with

\[ x = g_e(\hat x). \]

def interpolate(u):
    uinterp = np.array([])
    x = np.array([])

    #create the reference evaluation points
    xhat = np.linspace(0, 1)

    #loop over bases
    for e in range(n_el):
        # pick an element
        el = elements[e]

        # we want to sum, we initialize to zero
        uloc = np.zeros(xhat.shape)

        # sum over the local non-zero bases (2)
        for i in range(el["n_bases"]):
            # g_e^i
            glob_node = el["loc_2_glob"][i]
            # \phi_{g_e^i}
            loc_base = el["phi"][i]

            uloc += u[glob_node] * loc_base(xhat)

        uinterp = np.append(uinterp, uloc)
        # g_e(\hat x)
        x = np.append(x, el["gmapping"](xhat))

    return x, uinterp

We can generate a random vector $\vec{u}$ and use the previous function.

u = np.random.rand(n_bases)

x, uinterp = interpolate(u)

go.Figure(data=[
    go.Scatter(x=x, y=uinterp, mode='lines'),
    go.Scatter(x=s, y=u, mode='markers'),
])

Assembly¶

We are now ready the assemble the global stiffness matrix.

The local entries are \[ L^e_{i,j} = %\int_\Omega\nabla \phi_i \cdot \nabla \phi_j = \int_{\hat s} \frac{\nabla \hat\phi_{i} \cdot \nabla \hat\phi_{j}}{s_{e, 1} - s_{e, 0}} \]

which are then mapped to the global entries $g_e^i, g_e^j$.

Note that the integrals are performed with quadpy.

import quadpy
# some quadrature rule
scheme = quadpy.line_segment.gauss_patterson(5)

#triplets for the matrix
rows = []
cols = []
vals = []

for e in range(n_el):
    el = elements[e]

    for i in range(el["n_bases"]):
        for j in range(el["n_bases"]):
            # evaluation of the integral:
            # \int_{\hat s} \frac{\nabla \hat\phi_{i} \cdot \nabla \hat\phi_{j}}{s_{e, 1} - s_{e, 0}}
            val = scheme.integrate(
                lambda x:
                el["grad_phi"][i](x) * el["grad_phi"][j](x) / el["grad_gmapping"](x),
                [0.0, 1.0])

            # the local entry val at i, j goes to g_e^i, g_e^j
            rows.append(el["loc_2_glob"][i])
            cols.append(el["loc_2_glob"][j])
            vals.append(val)


rows = np.array(rows)
cols = np.array(cols)
vals = np.array(vals)

L = spr.coo_matrix((vals, (rows, cols)))
L = spr.csr_matrix(L)

L.toarray()
#this looks exacly like FD!

array([[ 10., -10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [-10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -10.,  20., -10.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -10.,  10.]])

We set the row 0 and n_el to identity for the boundary conditions

for bc in [0, n_el]:
    _, nnz = L[bc,:].nonzero()
    for j in nnz:
        if j != bc:
            L[bc, j] = 0.0
    L[bc, bc] = 1.0

L.A

array([[  1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [-10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0., -10.,  20., -10.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0., -10.,  20., -10.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   1.]])

We set the right-hand side to zero, and set the two boundary conditions to 1 and 4.

f = np.zeros((n_bases, 1))
f[0] = 1
f[-1] = 4

We now solve $L\vec{u}=f$ for $\vec{u}$.

u = spsolve(L, f)

We now plot the solution $\vec{u}$. Will we get a line?

x, uinterp = interpolate(u)


go.Figure(data=[
    go.Scatter(x=x, y=uinterp, mode='lines', name="solution"),
    go.Scatter(x=s, y=u, mode='markers', name="$u$"),
])

Mass Matrix¶

We change the pde from Laplace to Poisson

\[ -\Delta u = f \]

If we assume that $f$ is also expressed in terms of $\phi_i$ we can rewrite the weak form as

\[\sum_{i=0}^n u_i \int_\Omega\nabla \phi_i \cdot \nabla \phi_j = \sum_{i=0}^n f_i \int_\Omega\phi_i \phi_j, \qquad \forall j=0,\dots,n\]

Which can be represented in matrix form

\[ L \vec{u} = M \vec{f}, \]

where $\vec{f}$ is the vector of $f_i$ and

\[ M_{i,j} = \int_\Omega\phi_i \phi_j \]

is the mass matrix.

As for the stiffness matrix, the mass matrix can be localized

\[ M^e_{i,j} = \int_{\hat s_j} \hat\phi_i \cdot \hat\phi_j\,(s_{j, 1} - s_{j, 0}). \]

Note that there is no division since there are no gradients!

rows = []
cols = []
vals = []


#same as above but now we use phi instead of grad_phi and no division
for e in range(n_el):
    el = elements[e]

    for i in range(el["n_bases"]):
        for j in range(el["n_bases"]):
            # \int_{\hat s_j} \hat\phi_i \cdot \hat\phi_j\,(s_{j, 1} - s_{j, 0})
            val = scheme.integrate(
                lambda x:
                el["phi"][i](x) * el["phi"][j](x) * el["grad_gmapping"](x),
                [0.0, 1.0])

            rows.append(el["loc_2_glob"][i])
            cols.append(el["loc_2_glob"][j])
            vals.append(val)


rows = np.array(rows)
cols = np.array(cols)
vals = np.array(vals)

M = spr.coo_matrix((vals, (rows, cols)))
M = spr.csr_matrix(M)

Now we set $\vec{f}=4$ and zero boundary conditions.

f = 4*np.ones((n_bases, 1))
f = M*f

f[0] = 0
f[-1] = 0

We now solve $L\vec{u}=M\vec{f}$ for $\vec{u}$.

u = spsolve(L, f)

x, uinterp = interpolate(u)

go.Figure(data=[
    go.Scatter(x=x, y=uinterp, mode='lines', name="solution"),
    go.Scatter(x=s, y=u, mode='markers', name="$u$"),
])