# 2021-11-10 Finite elements¶

## Last time¶

• Stable numerics

• Galerkin optimality

• Boundary conditions

• Properties of Galerkin methods

## Today¶

• Element restriction

• Finite element assembly

• Nonlinear problems

• LibCEED

using Plots
using LinearAlgebra
using SparseArrays

default(linewidth=4)

function my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

function vander_legendre_deriv(x, k=nothing)
if isnothing(k)
k = length(x) # Square by default
end
m = length(x)
Q = ones(m, k)
dQ = zeros(m, k)
Q[:, 2] = x
dQ[:, 2] .= 1
for n in 1:k-2
Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
dQ[:, n+2] = (2*n + 1) * Q[:,n+1] + dQ[:,n]
end
Q, dQ
end

function febasis(P, Q)
x, _ = gausslobatto(P)
q, w = gausslegendre(Q)
Pk, _ = vander_legendre_deriv(x)
Bp, Dp = vander_legendre_deriv(q, P)
B = Bp / Pk
D = Dp / Pk
x, q, w, B, D
end

┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] # General form¶

$\int_\Omega v\cdot f_0(u, \nabla u) + \nabla v\cdot f_1(u, \nabla u) = 0, \quad \forall v$

we discretize as

$\sum_e \mathcal E_e^T \Big( B^T W \left\lvert \frac{\partial x}{\partial X} \right\rvert f_0(\tilde u, \nabla \tilde u) + D^T \left(\frac{\partial X}{\partial x}\right)^{T} W \left\lvert \frac{\partial x}{\partial X} \right\rvert f_1(\tilde u, \nabla\tilde u) \Big) = 0$
where $$\tilde u = B \mathcal E_e u$$ and $$\nabla \tilde u = \frac{\partial X}{\partial x} D \mathcal E_e u$$ are the values and gradients evaluated at quadrature points.

## Isoparametric mapping¶

Given the reference coordinates $$X \in \hat K \subset R^n$$ and physical coordinates $$x(X)$$, an integral on the physical element can be written

$\int_{K = x(\hat K)} f(x) dx = \int_{\hat K} \underbrace{\left\lvert \frac{\partial x}{\partial X} \right\rvert}_{\text{determinant}} f(x(X)) dX .$

Notation

Meaning

$$x$$

physical coordinates

$$X$$

reference coordinates

$$\mathcal E_e$$

restriction from global vector to element $$e$$

$$B$$

values of nodal basis functions at quadrature ponits on reference element

$$D$$

gradients of nodal basis functions at quadrature points on reference element

$$W$$

diagonal matrix of quadrature weights on reference element

$$\frac{\partial x}{\partial X} = D \mathcal E_e x$$

gradient of physical coordinates with respect to reference coordinates

$$\left\lvert \frac{\partial x}{\partial X}\right\rvert$$

determinant of coordinate transformation at each quadrature point

$$\frac{\partial X}{\partial x} = \left(\frac{\partial x}{\partial X}\right)^{-1}$$

derivative of reference coordinates with respect to physical coordinates

# Finite element mesh and restriction¶

function fe1_mesh(P, nelem)
x = LinRange(-1, 1, nelem+1)
rows = Int[]
cols = Int[]
for i in 1:nelem
append!(rows, (i-1)*P+1:i*P)
append!(cols, (i-1)*(P-1)+1:i*(P-1)+1)
end
x, sparse(cols, rows, ones(nelem*P))'
end
P, nelem = 4, 3
x, E = fe1_mesh(P, nelem)
sparse(E)

12×10 SparseMatrixCSC{Float64, Int64} with 12 stored entries:
1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅    ⋅
⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅    ⋅
⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅
⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅    ⋅
⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅    ⋅
⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0   ⋅
⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   1.0

function xnodal(x, P)
nelem = length(x) - 1
xn = Float64[]
xref, _ = gausslobatto(P)
for i in 1:nelem
xL, xR = x[i:i+1]
append!(xn, (xL+xR)/2 .+ (xR-xL)/2 * xref[1+(i>1):end])
end
xn
end
xn = xnodal(x, 4)
scatter(xn, zero, marker=:auto) # Finite element building blocks¶

struct FESpace
P::Int
Q::Int
nelem::Int
x::Vector
xn::Vector
q::Vector
w::Vector
B::Matrix
D::Matrix
function FESpace(P, Q, nelem)
x, Et = fe1_mesh(P, nelem)
xn = xnodal(x, P)
_, q, w, B, D = febasis(P, Q)
new(P, Q, nelem, x, xn, Et, q, w, B, D)
end
end

function fe_element(fe, e)
xL, xR = fe.x[e:e+1]
q = (xL+xR)/2 .+ (xR-xL)/2*fe.q
w = (xR - xL)/2 * fe.w
E = fe.Et'[:, (e-1)*fe.P+1:e*fe.P]'
dXdx = ones(fe.Q) * 2 / (xR - xL)
q, w, E, dXdx
end

fe_element (generic function with 1 method)


# Finite element residual assembly¶

$v^T F(u) \sim \int_\Omega v\cdot f_0(u, \nabla u) + \nabla v\cdot f_1(u, \nabla u) = 0, \quad \forall v$
$\sum_e \mathcal E_e^T \Big( B^T W \left\lvert \frac{\partial x}{\partial X} \right\rvert f_0(\tilde u, \nabla \tilde u) + D^T \left(\frac{\partial X}{\partial x}\right)^{T} W \left\lvert \frac{\partial x}{\partial X} \right\rvert f_1(\tilde u, \nabla\tilde u) \Big) = 0$

where $$\tilde u = B \mathcal E_e u$$ and $$\nabla \tilde u = \frac{\partial X}{\partial x} D \mathcal E_e u$$ are the values and gradients evaluated at quadrature points.

kappa(x) = 0.6 .+ 0.4*sin(pi*x/2)
fq(q, u, Du) = 0*u .- 1, kappa.(q) .* Du
dfq(q, u, du, Du, Ddu) = 0*du, kappa.(q) .* Ddu

function fe_residual(u_in, fe, fq; bci=, bcv=[1.])
u = copy(u_in); v = zero(u)
u[bci] = bcv
for e in 1:fe.nelem
q, w, E, dXdx = fe_element(fe, e)
B, D = fe.B, fe.D
ue = E * u
uq = B * ue
Duq = dXdx .* (D * ue)
f0, f1 = fq(q, uq, Duq)
ve = B' * (w .* f0) + D' * (dXdx .* w .* f1)
v += E' * ve
end
v[bci] = u_in[bci] - u[bci]
v
end

fe_residual (generic function with 1 method)

import NLsolve: nlsolve

fe = FESpace(3, 3, 5)
u0 = zero(fe.xn)
sol = nlsolve(u -> fe_residual(u, fe, fq), u0; method=:newton)
plot(fe.xn, sol.zero, marker=:auto) # Finite element Jacobian assembly¶

$v^T F(u) \sim \int_\Omega v\cdot f_0(u, \nabla u) + \nabla v\cdot f_1(u, \nabla u) = 0, \quad \forall v$
$v^T J(u) du \sim \int_\Omega v\cdot df_0(u, du, \nabla u, \nabla du) + \nabla v\cdot df_1(u, du, \nabla u, \nabla du) = 0, \quad \forall v$
function fe_jacobian(u_in, fe, dfq; bci=, bcv=[1.])
u = copy(u_in); u[bci] = bcv
rows, cols, vals = Int[], Int[], Float64[]
for e in 1:fe.nelem
q, w, E, dXdx = fe_element(fe, e)
B, D, P = fe.B, fe.D, fe.P
ue = E * u
uq = B * ue; Duq = dXdx .* (D * ue)
K = zeros(P, P)
for j in 1:fe.P
du = B[:,j]
Ddu = dXdx .* D[:,j]
df0, df1 = dfq(q, uq, du, Duq, Ddu)
K[:,j] = B' * (w .* df0) + D' * (dXdx .* w .* df1)
end
inds = rowvals(E')
append!(rows, kron(ones(P), inds))
append!(cols, kron(inds, ones(P)))
append!(vals, vec(K))
end
A = sparse(rows, cols, vals)
A[bci, :] .= 0; A[:, bci] .= 0
A[bci,bci] = diagm(ones(length(bci)))
A
end

fe_jacobian (generic function with 1 method)

sol = nlsolve(u -> fe_residual(u, fe, fq),
u -> fe_jacobian(u, fe, dfq),
u0;
method=:newton)

Results of Nonlinear Solver Algorithm
* Algorithm: Newton with line-search
* Starting Point: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
* Zero: [1.0, 2.8410183049139395, 4.238906485042849, 5.1880169525357305, 5.8160506219732895, 6.227791209068442, 6.501764009933425, 6.681877765936257, 6.795822293071308, 6.8589074098787215, 6.879070690065479]
* Inf-norm of residuals: 0.000000
* Iterations: 1
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 2
* Jacobian Calls (df/dx): 2


# Spy on the Jacobian¶

#@show cond(Matrix(J))

fe = FESpace(6, 6, 4)
u0 = zero(fe.xn)
J = fe_jacobian(u0, fe, dfq)
my_spy(J) • What would the matrix structure look like for a finite difference method that is 6th order accurate?

# A nonlinear problem¶

$\int_\Omega \nabla v (1 + u^2) \nabla u - v = 0$
function fq(q, u, Du)
f0 = 0*u .- 1
f1 = (1 .+ u.^2) .* Du
f0, f1
end
function dfq(q, u, du, Du, Ddu)
df0 = 0*du
df1 = (1 .+ u.^2) .* Ddu + 2u .* du .* Du
df0, df1
end
sol = nlsolve(u -> fe_residual(u, fe, fq),
u -> fe_jacobian(u, fe, dfq),
u0; method=:newton)

Results of Nonlinear Solver Algorithm
* Algorithm: Newton with line-search
* Starting Point: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
* Zero: [1.0, 1.0562613160921506, 1.1576325990202827, 1.2590479765200593, 1.3315488441072163, 1.3634502374658097, 1.3932371719548644, 1.4481281730989712, 1.5041997855247982, 1.5445730865975467, 1.5622871531149463, 1.5787333640569756, 1.6085978073405816, 1.6380216706859818, 1.6579945342008353, 1.666258096607518, 1.6735545064531623, 1.6855249632772669, 1.6947698658243535, 1.6984415353929925, 1.698885490572898]
* Inf-norm of residuals: 0.000000
* Iterations: 6
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 7
* Jacobian Calls (df/dx): 7

plot(fe.xn, sol.zero, marker=:auto) # Advection-diffusion (time independent)¶

(45)\begin{align} \nabla\cdot\big[ \mathbf w u - \kappa \nabla u \big] &= 1 & \int_\Omega \nabla v \cdot \Big[-\mathbf w u + \kappa \nabla u \big] &= \int_\Omega v 1, \forall v \end{align}
wind = 50
fq(q, u, Du) = -one.(u), -wind .* u + 1 * Du
N = length(fe.xn)
sol = nlsolve(u -> fe_residual(u, fe, fq; bci=[1, N], bcv=[0, 0]), zero(fe.xn), method=:newton)
plot(fe.xn, sol.zero, marker=:auto) # LibCEED¶

using LibCEED

ceed = Ceed()

Ceed
Ceed Resource: /cpu/self/opt/blocked
Preferred MemType: host

function fe1_ceed(ceed, P, nelem)
offsets = CeedInt[]
# libceed uses 0-based indexing
for e in 0:nelem-1
append!(offsets, e*(P-1):(e+1)*(P-1))
end
E = create_elem_restriction(ceed, nelem, P, 1, 1, nelem*(P-1)+1, offsets)
end

P, Q, nelem = 5, 5, 4
E = fe1_ceed(ceed, P, nelem)

CeedElemRestriction from (17, 1) to 4 elements with 5 nodes each and component stride 1

basis = create_tensor_h1_lagrange_basis(ceed, 1, 1, P, Q, GAUSS)

CeedBasis: dim=1 P=5 Q=5
qref1d:	 -0.90617985	 -0.53846931	  0.00000000	  0.53846931	  0.90617985
qweight1d:	  0.23692689	  0.47862867	  0.56888889	  0.47862867	  0.23692689
interp1d:	  0.59337070	  0.51643520	 -0.16382364	  0.08322282	 -0.02920508
interp1d:	 -0.10048257	  0.93136610	  0.22966726	 -0.09069490	  0.03014411
interp1d:	  0.00000000	  0.00000000	  1.00000000	  0.00000000	  0.00000000
interp1d:	  0.03014411	 -0.09069490	  0.22966726	  0.93136610	 -0.10048257
interp1d:	 -0.02920508	  0.08322282	 -0.16382364	  0.51643520	  0.59337070
grad1d:	 -3.70533645	  4.33282117	 -0.90392454	  0.42067623	 -0.14423641
grad1d:	 -0.52871527	 -1.09765797	  2.13259378	 -0.74973857	  0.24351802
grad1d:	  0.37500000	 -1.33658458	  0.00000000	  1.33658458	 -0.37500000
grad1d:	 -0.24351802	  0.74973857	 -2.13259378	  1.09765797	  0.52871527
grad1d:	  0.14423641	 -0.42067623	  0.90392454	 -4.33282117	  3.70533645


# libCEED QFunction and Operator¶

$\begin{gather*} v^T F(u) \sim \int_\Omega v \cdot \color{olive}{f_0(u, \nabla u)} + \nabla v \!:\! \color{olive}{f_1(u, \nabla u)} \quad v^T J w \sim \int_\Omega \begin{bmatrix} v \\ \nabla v \end{bmatrix}^T \color{teal}{\begin{bmatrix} f_{0,0} & f_{0,1} \\ f_{1,0} & f_{1,1} \end{bmatrix}} \begin{bmatrix} w \\ \nabla w \end{bmatrix} \\ u = B_I \mathcal E_e u_L \qquad \nabla u = \frac{\partial X}{\partial x} B_{\nabla} \mathcal E_e u_L \\ J w = \sum_e \mathcal E_e^T \begin{bmatrix} B_I \\ B_{\nabla} \end{bmatrix}^T \underbrace{\begin{bmatrix} I & \\ & \left( \frac{\partial X}{\partial x}\right)^T \end{bmatrix} W_q \color{teal}{\begin{bmatrix} f_{0,0} & f_{0,1} \\ f_{1,0} & f_{1,1} \end{bmatrix}} \begin{bmatrix} I & \\ & \left( \frac{\partial X}{\partial x}\right) \end{bmatrix}}_{\text{coefficients at quadrature points}} \begin{bmatrix} B_I \\ B_{\nabla} \end{bmatrix} \mathcal E_e w_L \end{gather*}$
@interior_qf setup_geom = (
ceed, dim=1,
(dxdX, :in, EVAL_GRAD, 1),
(w, :in, EVAL_WEIGHT),
(qdata, :out, EVAL_NONE, 1, 2),
begin
qdata = w * abs(dxdX)
qdata = w / abs(dxdX)
end
)

0

@interior_qf apply_diff = (
ceed, dim=1,
(qdata, :in, EVAL_NONE, 2),
(dudX, :in, EVAL_GRAD, 1),
(v, :out, EVAL_INTERP, 1),
(dvdX, :out, EVAL_GRAD, 1),
begin
v .= qdata * -1
dvdX .= qdata * dudX
end
)

0


# Wire up restrictions and operators¶

qdata = CeedVector(ceed, nelem*Q*2)
x = CeedVector(ceed, nelem+1)
x[] = LinRange(-1, 1, nelem+1)
Ex = fe1_ceed(ceed, 2, nelem)
Eq = create_elem_restriction_strided(ceed, nelem, Q, 2, nelem*Q*2, STRIDES_BACKEND)
Bx = create_tensor_h1_lagrange_basis(ceed, 1, 1, 2, Q, GAUSS)

op_setup = Operator(ceed, qf=setup_geom, fields=[
(:dxdX, Ex, Bx, CeedVectorActive()),
(:w, ElemRestrictionNone(), Bx, CeedVectorNone()),
(:qdata, Eq, BasisCollocated(), CeedVectorActive()),
])

CeedOperator
3 Fields
2 Input Fields:
Input Field :
Name: "dxdX"
Active vector
Input Field :
Name: "w"
No vector
1 Output Field:
Output Field :
Name: "qdata"
Collocated basis
Active vector

apply!(op_setup, x, qdata)

op_diff = Operator(ceed, qf=apply_diff, fields=[
(:qdata, Eq, BasisCollocated(), qdata),
(:dudX, E, basis, CeedVectorActive()),
(:v, E, basis, CeedVectorActive()),
(:dvdX, E, basis, CeedVectorActive()),
])

CeedOperator
4 Fields
2 Input Fields:
Input Field :
Name: "qdata"
Collocated basis
Input Field :
Name: "dudX"
Active vector
2 Output Fields:
Output Field :
Name: "v"
Active vector
Output Field :
Name: "dvdX"
Active vector


# Solve¶

function ceed_residual(u_in, op_diff; bci=, bcv=[1.])
u = copy(u_in); v = zero(u)
u[bci] = bcv
uceed = CeedVector(ceed, length(u))
setarray!(uceed, MEM_HOST, USE_POINTER, u)
vceed = CeedVector(ceed, length(u))
apply!(op_diff, uceed, vceed)
v = witharray(Vector, vceed)
v[bci] = u_in[bci] - u[bci]
v
end

N = getlvectorsize(E)
sol = nlsolve(u -> ceed_residual(u, op_diff; bci=[1, N], bcv=[0, 0]), ones(N), method=:newton)

Results of Nonlinear Solver Algorithm
* Algorithm: Newton with line-search
* Starting Point: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
* Zero: [-3.9475089863572066e-12, 0.08260957960239312, 0.21874999997080646, 0.32810470604827624, 0.3749999999584841, 0.4144412884027462, 0.4687499999651483, 0.4962729972393789, 0.5000000000237609, 0.4962729973023331, 0.4687500000248659, 0.41444128846070394, 0.3750000000810918, 0.3281047062415061, 0.21875000011356172, 0.08260957964624749, -3.9475089863572066e-12]
* Inf-norm of residuals: 0.000000
* Iterations: 1
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 2
* Jacobian Calls (df/dx): 2

plot(sol.zero, marker=:auto) • Why is this figure weird?

• How are the x nodes spaced, and how should they be spaced?