Skip to content

FEM Higher Order

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"

As in the linear case, the domain is \(\Omega = [0, 1]\), which we discretize with \(n_{el}\) segments (or elements) \(s_i\).

Now we also create the high-order nodes.

Note that we append them at the end. 

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

#number of bases and elements
n_el = 10


#Regular nodes, as before
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])

# we now pick the order
order = 3

nodes = s

#more bases
n_bases = n_el + 1 + n_el*(order-1)

#create the nodes for the plots
for e in range(n_el):
    #For every segment, we create order + 1 new high-order nodes
    tmp = np.linspace(s[e], s[e+1], num=order+1)
    #exclude the first and last since they already exists
    tmp = tmp[1:-1]

    #and append at the end of nodes
    nodes = np.append(nodes, tmp)

Plot, in orange the high-order nodes, in blue the linear nodes. 

go.Figure(data=[
    go.Scatter(x=s, y=np.zeros(s.shape), mode='lines+markers'),
    go.Scatter(x=nodes[n_el+1:], y=np.zeros(nodes.shape), mode='markers')
])
00.20.40.60.81−1−0.500.51trace 0trace 1

Local bases

As in the linear case, we define the reference element \(\hat s= [0, 1]\), a segment of unit length.

On each element we now have order+1 (e.g., 2 for linear, 3 for quadratic) non-zero local bases.

We define their “piece” on \(\hat s\).

It is important to respect the order of the nodes: the first 2 bases are always for the endpoints, and the others are ordered left to right.

Definition of linear bases, same as before 

def hat_phi_1_0(x):
    return 1-x
def hat_phi_1_1(x):
    return x

Definition of quadratic bases 

def hat_phi_2_0(x):
    return 2*(x-0.5)*(x-1)
def hat_phi_2_1(x):
    return 2*(x-0)*(x-0.5)
def hat_phi_2_2(x):
    return -4*(x-0.5)**2+1

Definition of cubic bases 

def hat_phi_3_0(x):
    return -9/2*(x-1/3)*(x-2/3)*(x-1)
def hat_phi_3_1(x):
    return 9/2*(x-0)*(x-1/3)*(x-2/3)
def hat_phi_3_2(x):
    return 27/2*(x-0)*(x-2/3)*(x-1)
def hat_phi_3_3(x):
    return -27/2*(x-0)*(x-1/3)*(x-1)

Utility function to return the list of functions 

def hat_phis(order):
    if order == 1:
        return [hat_phi_1_0, hat_phi_1_1]
    elif order == 2:
        return [hat_phi_2_0, hat_phi_2_1, hat_phi_2_2]
    elif order == 3:
        return [hat_phi_3_0, hat_phi_3_1, hat_phi_3_2, hat_phi_3_3]

We can now plot the order+1 local bases, same code as before.

Note that the first two bases correspond to the end-points, and the others are ordered. 

x = np.linspace(0, 1)

data = []

tmp = hat_phis(order)

for o in range(order+1):
    data.append(go.Scatter(x=x, y=tmp[o](x), mode='lines', name="$\hat\phi_{}$".format(o)))

go.Figure(data=data)
00.20.40.60.81−0.200.20.40.60.81

We use sympy to compute the gradients of the local bases. 

import sympy as sp

xsym = sp.Symbol('x')

def grad_hat_phis(order):
    #For linear we need to get the correct size
    if order == 1:
        return [lambda x : -np.ones(x.shape), lambda x : np.ones(x.shape)]

    res = []
    tmp = hat_phis(order)

    for fun in tmp:
        res.append(sp.lambdify(xsym, fun(xsym).diff(xsym)))
    return res

Plotting gradients 

x = np.linspace(0, 1)

data = []

tmp = grad_hat_phis(order)

for o in range(order+1):
    data.append(go.Scatter(x=x, y=tmp[o](x), mode='lines', name="$\hat\phi_{}$".format(o)))

go.Figure(data=data)
00.20.40.60.81−10−50510

Basis construction

This code is exacly as before.

The only difficulty is the local to global mapping: - the first 2 nodes are always the same

\[ g_e^0 = e \qquad\mathrm{and}\qquad g_e^1=g+1 \] - the others are

\[ g_e^i = n_{el} + e (\mathrm{order}-1) + i. \] 

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

    el["n_bases"] = order+1

    #2 bases
    el["phi"] = hat_phis(order)
    el["grad_phi"] = grad_hat_phis(order)

    #local to global mapping
    high_order_nodes = list(range(n_el + 1 + e*(order-1), n_el + e*(order-1) + order))
    el["loc_2_glob"] = [e, e+1] + high_order_nodes

    #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, as before. 

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

    xhat = np.linspace(0, 1)


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

        uloc = np.zeros(xhat.shape)

        for i in range(el["n_bases"]):
            glob_node = el["loc_2_glob"][i]
            loc_base = el["phi"][i]

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

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

    return x, uinterp

We can generate a random vector \(\vec{u}\) and use the previous function. This will interpolate all nodes. 

u = np.random.rand(n_bases)

x, uinterp = interpolate(u)


go.Figure(data=[
    go.Scatter(x=x, y=uinterp, mode='lines'),
    go.Scatter(x=nodes, y=u, mode='markers'),
])
00.20.40.60.8100.20.40.60.81trace 0trace 1

Assembly

We are now ready the assemble the global stiffness matrix, which is exacly as before. 

import quadpy
scheme = quadpy.line_segment.gauss_patterson(5)

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"]):
            val = scheme.integrate(
                lambda x:
                el["grad_phi"][i](x) * el["grad_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)

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

We set the rows 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

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[n_el] = 4

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

u = spsolve(L, f)

We now plot the solution. We expect a line, independently of order! 

x, uinterp = interpolate(u)

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

Mass Matrix

This is exactly as before! 

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"]):
            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[n_el] = 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=nodes, y=u, mode='markers', name="$u$"),
])
00.20.40.60.8100.10.20.30.40.5solution