# 2022-10-26 Finite Element Theory#

## Last time#

• Galerkin methods

• Weak forms

• $$L^2$$ Projection

• Galerkin “optimality”

## Today#

• Galerkin for Poisson

• Finite element constructions

using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays

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(x, k=nothing)
if isnothing(k)
k = length(x) # Square by default
end
m = length(x)
Q = ones(m, k)
Q[:, 2] = x
for n in 1:k-2
Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
end
Q
end

vander_legendre (generic function with 2 methods)


# Prime to nodal basis#

A nodal basis $$\{ \phi_j(x) \}$$ is one that satisfies $$n_i(\phi_j) = \delta_{ij} .$$

We write $$\phi_j$$ in the prime basis by solving with the generalized Vandermonde matrix $$V_{ij} = \langle n_i, p_j \rangle$$,

$\phi_j(x) = \sum_k p_k(x) (V^{-1})_{k,j} .$

k = 9
xn = LinRange(-1, 1, k)
V = vander_legendre(xn)
xx = LinRange(-1, 1, 50)
Pxx = vander_legendre(xx, k)
Pn = Pxx / V
plot(xx, Pn)
scatter!(xn, one.(xn))

using FastGaussQuadrature
xn, _ = gausslobatto(k)
V = vander_legendre(xn)
Pn = Pxx / V
plot(xx, Pn)
scatter!(xn, one.(xn))

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

vander_legendre_deriv (generic function with 2 methods)


# Galerkin method: $$L^2$$ projection#

A nice test problem that doesn’t require derivatives or boundary conditions: Find $$u \in \mathcal V_h$$ such that

$\int_{-1}^1 v(x) \big[ u(x) - f(x) \big] = 0, \quad \forall v \in \mathcal V_h$

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

function L2_galerkin(P, Q, f)
x, q, w, B, _ = febasis(P, Q)
A = B' * diagm(w) * B
rhs = B' * diagm(w) * f.(q)
u = A \ rhs
x, u
end

L2_galerkin (generic function with 1 method)

f(x) = tanh(5x)
P, Q = 16, 16
x, u = L2_galerkin(P, Q, f)
error = u - f.(x)
plot(x, u, marker=:circle, legend=:bottomright)
plot!(f, title="Galerkin error: $(norm(error))")  # Embed the solution in the continuous space# In thel last figure, the Galerkin solution was represented by a piecewise linear function, emphasizing its nodal values. Let’s make a nicer plot by multiplying by the corresponding nodal basis function. P, Q = 10,20 x, u = L2_galerkin(P, Q, f) error = u - f.(x) _, xx, _, BB, _ = febasis(P, 100) plot(f, label="exact", title="Galerkin error:$(norm(error))",
xlim=(-1, 1), legend=:bottomright)
plot!(xx, BB * u, label="Galerkin solution")
scatter!(x, u, marker=:circle, color=:black, label="Nodal values")

ns = 2 .^ (2:8)
function L2_error(n)
x, u = L2_galerkin(n, n, f)
norm(u - f.(x))
end
plot(ns, L2_error.(ns), marker=:circle, yscale=:log10,
title="Convergence", xlabel="number of points", ylabel="error")


# Galerkin error estimates#

We introduce the notation

$a(v, u) = \int_\Omega \nabla v(x) \cdot \nabla u(x)$
and note that $$a$$ is

• bilinear (linear in each of its arguments)

• symmetric: $$a(u, v) = a(v,u)$$

• positive definite: $$a(u, u) > 0$$ when $$u \ne 0$$ thus defines an inner product on the function space $$V$$.

We also introduce the $$L^2$$ inner product

$\langle u, v \rangle = \int_\Omega u(x) v(x)$
so that our continuous weak form is to find $$u \in V$$ such that
$a(v, u) = \langle v, f \rangle, \quad \forall v\in V.$
Our Galerkin discretization is to find $$u_h \in V_h \subset V$$ such that
$a(v_h, u_h) = \langle v_h, f \rangle, \quad \forall v_h \in V_h .$
Since $$V_h \subset V$$, we can subtract these two, yielding
$a(v_h, u_h - u) = 0, \quad \forall v_h \in V_h .$
This says that the error in the discrete solution $$u_h - u$$ is $$a$$-orthogonal to all test functions $$v_h$$.

# Galerkin optimality via energy norms#

We can also define the “energy norm” or $$a$$-norm,

$\lVert u \rVert_a = \sqrt{a(u,u)} .$
This norm satisfies the Cauchy-Schwarz inequality,
$\lvert a(u,v) \rvert \le \lVert u \rVert_a \lVert v \rVert_a .$
Now,

(49)#\begin{align} \lVert u_h - u \rVert_a^2 &= a(u_h - u, u_h - u) \\ &= a(u_h - v_h, u_h - u) + a(v_h - u, u_h - u) \\ &= a(v_h - u, u_h - u) \\ &\le \lVert v_h - u \rVert_a \lVert u_h - u \rVert_a . \end{align}

In other words,

$\lVert u_h - u \rVert_a \le \lVert v_h - u \rVert_a, \quad \forall v_h \in V_h .$
So the solution $$u_h$$ computed by the Galerkin discretization is optimal over the subspace $$V_h$$ as measured in the $$a$$-norm.

## Observations#

• The Galerkin method computes the exact solution any time it resides in the subspace $$V_h$$.

• The Galerkin method is automatically symmetric any time the weak form is symmetric.

• The Galerkin method can be spectrally accurate, similar to the Chebyshev finite difference methods.

• For a nonlinear problem, discretization and differentiation will commute.

# Galerkin method for Poisson#

$\int_{\Omega} \nabla v \cdot \kappa \cdot \nabla u = \int_{\Omega} v f, \forall v$
function poisson_galerkin(P, Q, kappa, f)
x, q, w, B, D = febasis(P, Q)
A = D' * diagm(w .* kappa.(q)) * D
rhs = B' * diagm(w) * f.(q)
# Boundary conditions
rhs[[1]] .= 0
A[1, :] .= 0
A[1, 1] = 1
#A[end, :] .= 0
#A[end, end] = 1
u = A \ rhs
x, u
end

poisson_galerkin (generic function with 1 method)

kappa(x) = 0.6 .+ 0.55 * sin.(pi*x)
P = 40
x, u = poisson_galerkin(P, P, one, x -> cos(x*pi/2))
uexact(x) = cos(x*pi/2) * (2/pi)^2
plot(uexact, title="Error: \$(norm(u - uexact.(x)))")
_, xx, _, BB, _ = febasis(P, 100)
plot!(xx, BB * u, label="FE solution", xlim=(-1, 1))
scatter!(x, u, label="Nodal solution", marker=:circle, color=:black,
legend=:bottomright)


# Galerkin preserves symmetry#

The bilinear operator is symmetric if $$a(v,u) = a(u,v)$$,

$a(v,u) = \int_{\Omega} \nabla v \cdot \kappa \cdot \nabla u$

This is symmetric even for

• variable coefficients

• any Neumann-type boundary conditions

• Dirichlet if eliminated or lifted symmetrically

• arbitrary grid spacing (and mesh topology in multiple dimensions)

x, q, w, B, D = febasis(7, 7)
A = D' * diagm(w .* kappa.(q)) * D
M = B' * diagm(w) * B
@show norm(A - A')
eigvals(A)

norm(A - A') = 6.869198314985155e-16

7-element Vector{Float64}:
-1.2257316920942458e-15
0.13096345516272143
1.2722189926037417
3.081972164593637
6.840904696876489
10.748215512493235
14.885725178270272


# 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_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$$

$$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(rows, cols, ones(nelem*P))
end
P, nelem = 4, 3
x, E = fe1_mesh(P, nelem)

(LinRange{Float64}(-1.0, 1.0, 4), sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], [1, 2, 3, 4, 4, 5, 6, 7, 7, 8, 9, 10], [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 12, 10))

function xnodal(x, P)
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)
end
xn
end
xn = xnodal(x, 5)
scatter(xn, zero, marker=:square)


# Finite element building blocks#

struct FESpace
P::Int
Q::Int
nelem::Int
x::Vector
xn::Vector
Et::SparseMatrixCSC{Float64, Int64}
q::Vector
w::Vector
B::Matrix
D::Matrix
function FESpace(P, Q, nelem)
x, Et = fe1_mesh(P, nelem)
@show Et
xn = xnodal(x, P)
_, q, w, B, D = febasis(P, Q)
new(P, Q, nelem, x, xn, Et, q, w, B, D)
end
end

# Extract out what we need for element e
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=[1], 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)

Et = sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], [1, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 9, 9, 10, 11], [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], 15, 11)

DimensionMismatch:

Stacktrace:
[1] mul!(C::Vector{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Vector{Float64}, α::Bool, β::Bool)
@ SparseArrays /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:31
[2] *
@ /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:51 [inlined]
[3] fe_residual(u_in::Vector{Float64}, fe::FESpace, fq::typeof(fq); bci::Vector{Int64}, bcv::Vector{Float64})
@ Main ./In[16]:11
[4] fe_residual(u_in::Vector{Float64}, fe::FESpace, fq::Function)
@ Main ./In[16]:5
[5] (::var"#4#5")(u::Vector{Float64})
@ Main ./In[17]:5
[6] (::NLSolversBase.var"#ff!#1"{var"#4#5"})(F::Vector{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/objective_types/inplace_factory.jl:11
[7] (::NLSolversBase.var"#fj_finitediff!#21"{NLSolversBase.var"#ff!#1"{var"#4#5"}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:central}(), Float64}})(F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/objective_types/oncedifferentiable.jl:138
[8] value_jacobian!!(obj::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/interface.jl:124
[9] value_jacobian!!
@ ~/.julia/packages/NLSolversBase/cfJrN/src/interface.jl:122 [inlined]
[10] newton_(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#27#29", cache::NLsolve.NewtonCache{Vector{Float64}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:48
[11] #newton#7
@ ~/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:146 [inlined]
[12] nlsolve(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:23
[13] nlsolve(f::Function, initial_x::Vector{Float64}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:52
[14] top-level scope
@ In[17]:5
[15] eval
@ ./boot.jl:368 [inlined]
[16] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


# 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=[1], 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)

DimensionMismatch:

Stacktrace:
[1] mul!(C::Vector{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Vector{Float64}, α::Bool, β::Bool)
@ SparseArrays /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:31
[2] *
@ /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:51 [inlined]
[3] fe_jacobian(u_in::Vector{Float64}, fe::FESpace, dfq::typeof(dfq); bci::Vector{Int64}, bcv::Vector{Float64})
@ Main ./In[18]:7
[4] fe_jacobian(u_in::Vector{Float64}, fe::FESpace, dfq::Function)
@ Main ./In[18]:1
[5] (::var"#8#10")(u::Vector{Float64})
@ Main ./In[19]:2
[6] f!
@ ~/.julia/packages/NLsolve/gJL1I/src/objectives/helpers.jl:6 [inlined]
[7] (::NLSolversBase.var"#fj!#7"{NLsolve.var"#f!#1"{var"#7#9"}, NLsolve.var"#f!#1"{var"#8#10"}})(fx::Vector{Float64}, jx::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/objective_types/abstract.jl:7
[8] value_jacobian!!(obj::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, F::Vector{Float64}, J::Matrix{Float64}, x::Vector{Float64})
@ NLSolversBase ~/.julia/packages/NLSolversBase/cfJrN/src/interface.jl:124
[9] value_jacobian!!
@ ~/.julia/packages/NLSolversBase/cfJrN/src/interface.jl:122 [inlined]
[10] newton_(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#27#29", cache::NLsolve.NewtonCache{Vector{Float64}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:48
[11] #newton#7
@ ~/.julia/packages/NLsolve/gJL1I/src/solvers/newton.jl:146 [inlined]
[12] nlsolve(df::NLSolversBase.OnceDifferentiable{Vector{Float64}, Matrix{Float64}, Vector{Float64}}, initial_x::Vector{Float64}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::LineSearches.Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64)
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:23
[13] nlsolve(f::Function, j::Function, initial_x::Vector{Float64}; inplace::Bool, kwargs::Base.Pairs{Symbol, Symbol, Tuple{Symbol}, NamedTuple{(:method,), Tuple{Symbol}}})
@ NLsolve ~/.julia/packages/NLsolve/gJL1I/src/nlsolve/nlsolve.jl:67
[14] top-level scope
@ In[19]:1
[15] eval
@ ./boot.jl:368 [inlined]
[16] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

fe = FESpace(3, 3, 5)
u0 = zero(fe.xn)
function mass(q, u, du, Du, Ddu)
df0 = du
df1 = 0*du
df0, df1
end
fe_jacobian(zero.(fe.xn), fe, mass, bci=[], bcv=[])

Et = sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], [1, 2, 3, 3, 4, 5, 5, 6, 7, 7, 8, 9, 9, 10, 11], [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], 15, 11)

DimensionMismatch:

Stacktrace:
[1] mul!(C::Vector{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Vector{Float64}, α::Bool, β::Bool)
@ SparseArrays /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:31
[2] *
@ /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:51 [inlined]
[3] fe_jacobian(u_in::Vector{Float64}, fe::FESpace, dfq::typeof(mass); bci::Vector{Any}, bcv::Vector{Any})
@ Main ./In[18]:7
[4] top-level scope
@ In[20]:8
[5] eval
@ ./boot.jl:368 [inlined]
[6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

B = diagm(11, 10, 0 => ones(10))
B[end,:] = B[1,:]
B

11×10 Matrix{Float64}:
1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

sigma = 0.2
exp.(-fe.xn[1:end-1].^2/sigma^2)

8-element Vector{Float64}:
1.388794386496407e-11
1.1253517471925912e-7
0.00012340980408667956
0.00012340980408667956
0.01831563888873418
0.36787944117144256
0.36787944117144256
1.0

exp(A)

7×7 Matrix{Float64}:
3.46915e8  -5.63207e8   3.53942e8  …   2.97179e8  -4.30059e8   2.58474e8
-5.63207e8   9.1506e8   -5.76891e8     -4.91434e8   7.14607e8  -4.30059e8
3.53942e8  -5.76891e8   3.68403e8      3.31971e8  -4.91434e8   2.97179e8
-2.63244e8   4.31923e8  -2.83171e8     -2.83171e8   4.31923e8  -2.63244e8
2.97179e8  -4.91434e8   3.31971e8      3.68403e8  -5.76891e8   3.53942e8
-4.30059e8   7.14607e8  -4.91434e8  …  -5.76891e8   9.1506e8   -5.63207e8
2.58474e8  -4.30059e8   2.97179e8      3.53942e8  -5.63207e8   3.46915e8

using LinearAlgebra


# 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)

Et = sparse([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], [1, 2, 3, 4, 5, 6, 6, 7, 8, 9, 10, 11, 11, 12, 13, 14, 15, 16, 16, 17, 18, 19, 20, 21], [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, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 24, 21)

DimensionMismatch:

Stacktrace:
[1] mul!(C::Vector{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Vector{Float64}, α::Bool, β::Bool)
@ SparseArrays /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:31
[2] *
@ /usr/share/julia/stdlib/v1.8/SparseArrays/src/linalg.jl:51 [inlined]
[3] fe_jacobian(u_in::Vector{Float64}, fe::FESpace, dfq::typeof(dfq); bci::Vector{Int64}, bcv::Vector{Float64})
@ Main ./In[18]:7
[4] fe_jacobian(u_in::Vector{Float64}, fe::FESpace, dfq::Function)
@ Main ./In[18]:1
[5] top-level scope
@ In[25]:5
[6] eval
@ ./boot.jl:368 [inlined]
[7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)