{
"cells": [
{
"cell_type": "markdown",
"id": "3471b94b",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2022-10-26 Finite Element Theory"
]
},
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"## Last time\n",
"\n",
"* Questions about projects\n",
"* Galerkin methods\n",
" * Weak forms\n",
" * $L^2$ Projection\n",
" * Galerkin \"optimality\" \n",
" \n",
"## Today\n",
"* Galerkin for Poisson\n",
"* Finite element constructions"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "e82d841d",
"metadata": {
"hideOutput": true,
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"vander_legendre (generic function with 2 methods)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"default(linewidth=3)\n",
"using LinearAlgebra\n",
"using SparseArrays\n",
"\n",
"function my_spy(A)\n",
" cmax = norm(vec(A), Inf)\n",
" s = max(1, ceil(120 / size(A, 1)))\n",
" spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))\n",
"end\n",
"\n",
"\n",
"function vander_legendre(x, k=nothing)\n",
" if isnothing(k)\n",
" k = length(x) # Square by default\n",
" end\n",
" m = length(x)\n",
" Q = ones(m, k)\n",
" Q[:, 2] = x\n",
" for n in 1:k-2\n",
" Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)\n",
" end\n",
" Q\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "1657e91e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Prime to nodal basis\n",
"\n",
"A **nodal basis** $\\{ \\phi_j(x) \\}$ is one that satisfies\n",
"$ n_i(\\phi_j) = \\delta_{ij} . $\n",
"\n",
"We write $\\phi_j$ in the prime basis by solving with the generalized Vandermonde matrix $V_{ij} = \\langle n_i, p_j \\rangle$,\n",
"$$ \\phi_j(x) = \\sum_k p_k(x) (V^{-1})_{k,j} . $$"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "3fb34a25",
"metadata": {
"cell_style": "split",
"scrolled": true,
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 9\n",
"xn = LinRange(-1, 1, k) \n",
"V = vander_legendre(xn)\n",
"xx = LinRange(-1, 1, 50)\n",
"Pxx = vander_legendre(xx, k)\n",
"Pn = Pxx / V\n",
"plot(xx, Pn)\n",
"scatter!(xn, one.(xn))"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "1efc623f",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using FastGaussQuadrature\n",
"xn, _ = gausslobatto(k)\n",
"V = vander_legendre(xn)\n",
"Pn = Pxx / V\n",
"plot(xx, Pn)\n",
"scatter!(xn, one.(xn))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "230e8866",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"vander_legendre_deriv (generic function with 2 methods)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function vander_legendre_deriv(x, k=nothing)\n",
" if isnothing(k)\n",
" k = length(x) # Square by default\n",
" end\n",
" m = length(x)\n",
" Q = ones(m, k)\n",
" dQ = zeros(m, k)\n",
" Q[:, 2] = x\n",
" dQ[:, 2] .= 1\n",
" for n in 1:k-2\n",
" Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)\n",
" dQ[:, n+2] = (2*n + 1) * Q[:,n+1] + dQ[:,n]\n",
" end\n",
" Q, dQ\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "ed8f0dc4",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Galerkin method: $L^2$ projection\n",
"\n",
"A nice test problem that doesn't require derivatives or boundary conditions:\n",
"Find $u \\in \\mathcal V_h$ such that\n",
"$$ \\int_{-1}^1 v(x) \\big[ u(x) - f(x) \\big] = 0, \\quad \\forall v \\in \\mathcal V_h$$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "96ec301b",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"L2_galerkin (generic function with 1 method)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function febasis(P, Q)\n",
" x, _ = gausslobatto(P)\n",
" q, w = gausslegendre(Q)\n",
" Pk, _ = vander_legendre_deriv(x)\n",
" Bp, Dp = vander_legendre_deriv(q, P)\n",
" B = Bp / Pk\n",
" D = Dp / Pk\n",
" x, q, w, B, D\n",
"end\n",
"\n",
"function L2_galerkin(P, Q, f)\n",
" x, q, w, B, _ = febasis(P, Q)\n",
" A = B' * diagm(w) * B\n",
" rhs = B' * diagm(w) * f.(q)\n",
" u = A \\ rhs\n",
" x, u\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 48,
"id": "139878d2",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(x) = tanh(5x)\n",
"P, Q = 16, 16\n",
"x, u = L2_galerkin(P, Q, f)\n",
"error = u - f.(x)\n",
"plot(x, u, marker=:circle, legend=:bottomright)\n",
"plot!(f, title=\"Galerkin error: $(norm(error))\")"
]
},
{
"cell_type": "markdown",
"id": "63fecb29",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Embed the solution in the continuous space\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 51,
"id": "7bf49b06",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P, Q = 10,20\n",
"x, u = L2_galerkin(P, Q, f)\n",
"error = u - f.(x)\n",
"_, xx, _, BB, _ = febasis(P, 100)\n",
"plot(f, label=\"exact\", title=\"Galerkin error: $(norm(error))\",\n",
" xlim=(-1, 1), legend=:bottomright)\n",
"plot!(xx, BB * u, label=\"Galerkin solution\")\n",
"scatter!(x, u, marker=:circle, color=:black, label=\"Nodal values\")"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "78f13bca",
"metadata": {
"cell_style": "split",
"hideOutput": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ns = 2 .^ (2:8)\n",
"function L2_error(n)\n",
" x, u = L2_galerkin(n, n, f)\n",
" norm(u - f.(x))\n",
"end\n",
"plot(ns, L2_error.(ns), marker=:circle, yscale=:log10,\n",
" title=\"Convergence\", xlabel=\"number of points\", ylabel=\"error\")"
]
},
{
"cell_type": "markdown",
"id": "a27da104",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Galerkin error estimates"
]
},
{
"cell_type": "markdown",
"id": "ff2b2754",
"metadata": {
"cell_style": "split"
},
"source": [
"We introduce the notation\n",
"$$ a(v, u) = \\int_\\Omega \\nabla v(x) \\cdot \\nabla u(x) $$\n",
"and note that $a$ is\n",
"* bilinear (linear in each of its arguments)\n",
"* symmetric: $a(u, v) = a(v,u)$\n",
"* positive definite: $a(u, u) > 0$ when $u \\ne 0$\n",
"thus defines an inner product on the function space $V$."
]
},
{
"cell_type": "markdown",
"id": "4966b4c5",
"metadata": {
"cell_style": "split"
},
"source": [
"We also introduce the $L^2$ inner product\n",
"$$ \\langle u, v \\rangle = \\int_\\Omega u(x) v(x) $$\n",
"so that our continuous weak form is to find $u \\in V$ such that\n",
"$$ a(v, u) = \\langle v, f \\rangle, \\quad \\forall v\\in V. $$\n",
"Our Galerkin discretization is to find $u_h \\in V_h \\subset V$ such that\n",
"$$ a(v_h, u_h) = \\langle v_h, f \\rangle, \\quad \\forall v_h \\in V_h . $$\n",
"Since $V_h \\subset V$, we can subtract these two, yielding\n",
"$$ a(v_h, u_h - u) = 0, \\quad \\forall v_h \\in V_h .$$\n",
"This says that the error in the discrete solution $u_h - u$ is $a$-orthogonal to all test functions $v_h$."
]
},
{
"cell_type": "markdown",
"id": "842f1ca8",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Galerkin optimality via energy norms"
]
},
{
"cell_type": "markdown",
"id": "fba20ff3",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"We can also define the \"energy norm\" or $a$-norm,\n",
"$$ \\lVert u \\rVert_a = \\sqrt{a(u,u)} . $$\n",
"This norm satisfies the Cauchy-Schwarz inequality,\n",
"$$ \\lvert a(u,v) \\rvert \\le \\lVert u \\rVert_a \\lVert v \\rVert_a . $$\n",
"Now,\n",
"\\begin{align}\n",
"\\lVert u_h - u \\rVert_a^2 &= a(u_h - u, u_h - u) \\\\\n",
"&= a(u_h - v_h, u_h - u) + a(v_h - u, u_h - u) \\\\\n",
"&= a(v_h - u, u_h - u) \\\\\n",
"&\\le \\lVert v_h - u \\rVert_a \\lVert u_h - u \\rVert_a .\n",
"\\end{align}\n",
"In other words,\n",
"$$\\lVert u_h - u \\rVert_a \\le \\lVert v_h - u \\rVert_a, \\quad \\forall v_h \\in V_h .$$\n",
"So the solution $u_h$ computed by the Galerkin discretization is optimal over the subspace $V_h$ as measured in the $a$-norm."
]
},
{
"cell_type": "markdown",
"id": "58a462a3",
"metadata": {
"cell_style": "split"
},
"source": [
"## Observations\n",
"* The Galerkin method computes the exact solution any time it resides in the subspace $V_h$.\n",
"* The Galerkin method is automatically symmetric any time the weak form is symmetric.\n",
"* The Galerkin method can be spectrally accurate, similar to the Chebyshev finite difference methods.\n",
"* For a nonlinear problem, discretization and differentiation will commute."
]
},
{
"cell_type": "markdown",
"id": "2da97560",
"metadata": {
"hidePrompt": true,
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Galerkin method for Poisson\n",
"\n",
"$$ \\int_{\\Omega} \\nabla v \\cdot \\kappa \\cdot \\nabla u = \\int_{\\Omega} v f, \\forall v$$"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "5250e6e0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"poisson_galerkin (generic function with 1 method)"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function poisson_galerkin(P, Q, kappa, f)\n",
" x, q, w, B, D = febasis(P, Q)\n",
" A = D' * diagm(w .* kappa.(q)) * D\n",
" rhs = B' * diagm(w) * f.(q)\n",
" # Boundary conditions\n",
" rhs[[1]] .= 0\n",
" A[1, :] .= 0\n",
" A[1, 1] = 1\n",
" #A[end, :] .= 0\n",
" #A[end, end] = 1\n",
" u = A \\ rhs\n",
" x, u\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 65,
"id": "bc10831f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kappa(x) = 0.6 .+ 0.55 * sin.(pi*x)\n",
"P = 40\n",
"x, u = poisson_galerkin(P, P, one, x -> cos(x*pi/2))\n",
"uexact(x) = cos(x*pi/2) * (2/pi)^2\n",
"plot(uexact, title=\"Error: $(norm(u - uexact.(x)))\")\n",
"_, xx, _, BB, _ = febasis(P, 100)\n",
"plot!(xx, BB * u, label=\"FE solution\", xlim=(-1, 1))\n",
"scatter!(x, u, label=\"Nodal solution\", marker=:circle, color=:black,\n",
" legend=:bottomright)"
]
},
{
"cell_type": "markdown",
"id": "0403d1a2",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## What if we only impose the boundary condition on one side?"
]
},
{
"cell_type": "markdown",
"id": "ec253cd3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Galerkin preserves symmetry"
]
},
{
"cell_type": "markdown",
"id": "dd93aadd",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"The bilinear operator is symmetric if $a(v,u) = a(u,v)$,\n",
"\n",
"$$ a(v,u) = \\int_{\\Omega} \\nabla v \\cdot \\kappa \\cdot \\nabla u $$\n",
"\n",
"This is symmetric even for\n",
"* variable coefficients\n",
"* any Neumann-type boundary conditions\n",
"* Dirichlet if eliminated or lifted symmetrically\n",
"* arbitrary grid spacing (and mesh topology in multiple dimensions)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"id": "0f2a45e6",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm(A - A') = 6.869198314985155e-16\n"
]
},
{
"data": {
"text/plain": [
"7-element Vector{Float64}:\n",
" -1.2257316920942458e-15\n",
" 0.13096345516272143\n",
" 1.2722189926037417\n",
" 3.081972164593637\n",
" 6.840904696876489\n",
" 10.748215512493235\n",
" 14.885725178270272"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, q, w, B, D = febasis(7, 7)\n",
"A = D' * diagm(w .* kappa.(q)) * D\n",
"M = B' * diagm(w) * B\n",
"@show norm(A - A')\n",
"eigvals(A)"
]
},
{
"cell_type": "markdown",
"id": "51500e5c",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## What if we use too few quadrature points?"
]
},
{
"cell_type": "markdown",
"id": "4e012df5",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"![](../img/libCEED-2-trim.svg)"
]
},
{
"cell_type": "markdown",
"id": "8061be0c",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# General form\n",
"$$ \\int_\\Omega v\\cdot f_0(u, \\nabla u) + \\nabla v\\cdot f_1(u, \\nabla u) = 0, \\quad \\forall v$$\n",
"we discretize as\n",
"$$ \\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 $$\n",
"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."
]
},
{
"cell_type": "markdown",
"id": "a04a2180",
"metadata": {
"cell_style": "split"
},
"source": [
"## Isoparametric mapping\n",
"\n",
"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\n",
"$$ \\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 .$$"
]
},
{
"cell_type": "markdown",
"id": "f81921c7",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"| Notation | Meaning |\n",
"|---------|:-------------|\n",
"| $x$ | physical coordinates |\n",
"| $X$ | reference coordinates |\n",
"| $\\mathcal E_e$ | restriction from global vector to element $e$ |\n",
"| $B$ | values of nodal basis functions at quadrature ponits on reference element |\n",
"| $D$ | gradients of nodal basis functions at quadrature points on reference element|\n",
"| $W$ | diagonal matrix of quadrature weights on reference element |\n",
"| $\\frac{\\partial x}{\\partial X} = D \\mathcal E_e x $ | gradient of physical coordinates with respect to reference coordinates |\n",
"| $\\left\\lvert \\frac{\\partial x}{\\partial X}\\right\\rvert$ | determinant of coordinate transformation at each quadrature point |\n",
"| $\\frac{\\partial X}{\\partial x} = \\left(\\frac{\\partial x}{\\partial X}\\right)^{-1}$ | derivative of reference coordinates with respect to physical coordinates |"
]
},
{
"cell_type": "markdown",
"id": "a140a545",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Finite element mesh and restriction"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "109306ab",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"(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))"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function fe1_mesh(P, nelem)\n",
" x = LinRange(-1, 1, nelem+1)\n",
" rows = Int[]\n",
" cols = Int[]\n",
" for i in 1:nelem\n",
" append!(rows, (i-1)*P+1:i*P)\n",
" append!(cols, (i-1)*(P-1)+1:i*(P-1)+1)\n",
" end\n",
" x, sparse(rows, cols, ones(nelem*P))\n",
"end\n",
"P, nelem = 4, 3\n",
"x, E = fe1_mesh(P, nelem)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "9c5ffff5",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function xnodal(x, P)\n",
" xn = Float64[]\n",
" xref, _ = gausslobatto(P)\n",
" for i in 1:nelem\n",
" xL, xR = x[i:i+1]\n",
" append!(xn, (xL+xR)/2 .+ (xR-xL)/2 * xref)\n",
" end\n",
" xn\n",
"end\n",
"xn = xnodal(x, 5)\n",
"scatter(xn, zero, marker=:square)"
]
},
{
"cell_type": "markdown",
"id": "72c23be8",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Finite element building blocks"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "c7a39135",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [],
"source": [
"struct FESpace\n",
" P::Int\n",
" Q::Int\n",
" nelem::Int\n",
" x::Vector\n",
" xn::Vector\n",
" Et::SparseMatrixCSC{Float64, Int64}\n",
" q::Vector\n",
" w::Vector\n",
" B::Matrix\n",
" D::Matrix\n",
" function FESpace(P, Q, nelem)\n",
" x, Et = fe1_mesh(P, nelem)\n",
" @show Et\n",
" xn = xnodal(x, P)\n",
" _, q, w, B, D = febasis(P, Q)\n",
" new(P, Q, nelem, x, xn, Et, q, w, B, D)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "44ef27d9",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"fe_element (generic function with 1 method)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Extract out what we need for element e\n",
"function fe_element(fe, e)\n",
" xL, xR = fe.x[e:e+1]\n",
" q = (xL+xR)/2 .+ (xR-xL)/2*fe.q\n",
" w = (xR - xL)/2 * fe.w\n",
" E = fe.Et[(e-1)*fe.P+1:e*fe.P, :]\n",
" dXdx = ones(fe.Q) * 2 / (xR - xL)\n",
" q, w, E, dXdx\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "11bc0cda",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Finite element residual assembly\n",
"$$ 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$$\n",
"$$ \\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 $$\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "6066e3a7",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"fe_residual (generic function with 1 method)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kappa(x) = 0.6 .+ 0.4*sin(pi*x/2)\n",
"fq(q, u, Du) = 0*u .- 1, kappa.(q) .* Du\n",
"dfq(q, u, du, Du, Ddu) = 0*du, kappa.(q) .* Ddu\n",
"\n",
"function fe_residual(u_in, fe, fq; bci=[1], bcv=[1.])\n",
" u = copy(u_in); v = zero(u)\n",
" u[bci] = bcv\n",
" for e in 1:fe.nelem\n",
" q, w, E, dXdx = fe_element(fe, e)\n",
" B, D = fe.B, fe.D\n",
" ue = E * u\n",
" uq = B * ue\n",
" Duq = dXdx .* (D * ue)\n",
" f0, f1 = fq(q, uq, Duq)\n",
" ve = B' * (w .* f0) + D' * (dXdx .* w .* f1)\n",
" v += E' * ve\n",
" end\n",
" v[bci] = u_in[bci] - u[bci]\n",
" v\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "cf04d909",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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)\n"
]
},
{
"ename": "LoadError",
"evalue": "DimensionMismatch: ",
"output_type": "error",
"traceback": [
"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)",
" @ Base ./loading.jl:1428"
]
}
],
"source": [
"import NLsolve: nlsolve\n",
"\n",
"fe = FESpace(3, 3, 5)\n",
"u0 = zero(fe.xn)\n",
"sol = nlsolve(u -> fe_residual(u, fe, fq), u0; method=:newton)\n",
"plot(fe.xn, sol.zero, marker=:auto)"
]
},
{
"cell_type": "markdown",
"id": "4fa8f8e9",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Finite element Jacobian assembly\n",
"\n",
"$$ 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$$\n",
"$$ 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$$"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "de0b0f51",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"fe_jacobian (generic function with 1 method)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function fe_jacobian(u_in, fe, dfq; bci=[1], bcv=[1.])\n",
" u = copy(u_in); u[bci] = bcv\n",
" rows, cols, vals = Int[], Int[], Float64[]\n",
" for e in 1:fe.nelem\n",
" q, w, E, dXdx = fe_element(fe, e)\n",
" B, D, P = fe.B, fe.D, fe.P\n",
" ue = E * u\n",
" uq = B * ue; Duq = dXdx .* (D * ue)\n",
" K = zeros(P, P)\n",
" for j in 1:fe.P\n",
" du = B[:,j]\n",
" Ddu = dXdx .* D[:,j]\n",
" df0, df1 = dfq(q, uq, du, Duq, Ddu)\n",
" K[:,j] = B' * (w .* df0) + D' * (dXdx .* w .* df1)\n",
" end\n",
" inds = rowvals(E')\n",
" append!(rows, kron(ones(P), inds))\n",
" append!(cols, kron(inds, ones(P)))\n",
" append!(vals, vec(K))\n",
" end\n",
" A = sparse(rows, cols, vals)\n",
" A[bci, :] .= 0; A[:, bci] .= 0\n",
" A[bci,bci] = diagm(ones(length(bci)))\n",
" A\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "d585a0be",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"ename": "LoadError",
"evalue": "DimensionMismatch: ",
"output_type": "error",
"traceback": [
"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)",
" @ Base ./loading.jl:1428"
]
}
],
"source": [
"sol = nlsolve(u -> fe_residual(u, fe, fq),\n",
" u -> fe_jacobian(u, fe, dfq),\n",
" u0;\n",
" method=:newton)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "362d82e1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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)\n"
]
},
{
"ename": "LoadError",
"evalue": "DimensionMismatch: ",
"output_type": "error",
"traceback": [
"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)",
" @ Base ./loading.jl:1428"
]
}
],
"source": [
"fe = FESpace(3, 3, 5)\n",
"u0 = zero(fe.xn)\n",
"function mass(q, u, du, Du, Ddu)\n",
" df0 = du\n",
" df1 = 0*du\n",
" df0, df1\n",
"end\n",
"fe_jacobian(zero.(fe.xn), fe, mass, bci=[], bcv=[])"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "7e19a574",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"11×10 Matrix{Float64}:\n",
" 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0\n",
" 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = diagm(11, 10, 0 => ones(10))\n",
"B[end,:] = B[1,:]\n",
"B"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "5aa9269d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8-element Vector{Float64}:\n",
" 1.388794386496407e-11\n",
" 1.1253517471925912e-7\n",
" 0.00012340980408667956\n",
" 0.00012340980408667956\n",
" 0.01831563888873418\n",
" 0.36787944117144256\n",
" 0.36787944117144256\n",
" 1.0"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sigma = 0.2\n",
"exp.(-fe.xn[1:end-1].^2/sigma^2)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "2f6484fd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"7×7 Matrix{Float64}:\n",
" 3.46915e8 -5.63207e8 3.53942e8 … 2.97179e8 -4.30059e8 2.58474e8\n",
" -5.63207e8 9.1506e8 -5.76891e8 -4.91434e8 7.14607e8 -4.30059e8\n",
" 3.53942e8 -5.76891e8 3.68403e8 3.31971e8 -4.91434e8 2.97179e8\n",
" -2.63244e8 4.31923e8 -2.83171e8 -2.83171e8 4.31923e8 -2.63244e8\n",
" 2.97179e8 -4.91434e8 3.31971e8 3.68403e8 -5.76891e8 3.53942e8\n",
" -4.30059e8 7.14607e8 -4.91434e8 … -5.76891e8 9.1506e8 -5.63207e8\n",
" 2.58474e8 -4.30059e8 2.97179e8 3.53942e8 -5.63207e8 3.46915e8"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"exp(A)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "e35322bc",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "markdown",
"id": "3e43ea8c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Spy on the Jacobian"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "8e87d631",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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)\n"
]
},
{
"ename": "LoadError",
"evalue": "DimensionMismatch: ",
"output_type": "error",
"traceback": [
"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)",
" @ Base ./loading.jl:1428"
]
}
],
"source": [
"#@show cond(Matrix(J))\n",
"\n",
"fe = FESpace(6, 6, 4)\n",
"u0 = zero(fe.xn)\n",
"J = fe_jacobian(u0, fe, dfq)\n",
"my_spy(J)"
]
},
{
"cell_type": "markdown",
"id": "b6c7a911",
"metadata": {
"cell_style": "split"
},
"source": [
"* What is interesting about this matrix structure?\n",
" * What would the matrix structure look like for a finite difference method that is 6th order accurate?"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "21d9247c",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.7.2",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.2"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}