# The Finite Element Method in One Dimension {#sec-one-dimensional-fem}
```{=latex}
\chaptershorttitle{One-Dimensional FEM}
```
Chapter 4 derived the variational problem for a one-dimensional axial bar. We
now use that problem to develop the finite element method from beginning to
end. The one-dimensional setting is especially useful because the main ideas
can be seen without the geometric complications of two- and three-dimensional
meshes. We can construct the approximation functions explicitly, see where the
algebraic equations come from, follow an element contribution into the global
system, and recover the displacement, strain, stress, and reactions after the
solve.
The discussion is detailed because these ideas are reused throughout the rest
of the book. Chapter 6 will change the element geometry and introduce
multidimensional mappings, but it will not change the basic sequence developed
here.
## Galerkin approximation and the finite element idea {#sec-c5-galerkin}
Consider the axial bar on $(0,L)$ with a prescribed displacement at the left
endpoint, a distributed axial load $f(x)$, and an applied axial force
$\overline P$ at the right endpoint. We first take the prescribed displacement
to be zero. The admissible functions therefore belong to
$$
\mathcal{V}_0
=
\{v\in H^1(0,L):v(0)=0\}.
$$
The variational problem from Chapter 4 is
$$
\text{find }u\in\mathcal{V}_0\text{ such that}\qquad
a(u,v)=\ell(v)
\qquad\forall v\in\mathcal{V}_0,
$$ {#eq-c5-continuous-variational}
where
$$
a(u,v)=\int_0^L AE\,u'v'\,\dd x,
\qquad
\ell(v)=\int_0^L fv\,\dd x+\overline P\,v(L).
$$ {#eq-c5-bar-forms}
{#fig-c5-bar-loading width="86%" fig-align="center"}
The unknown in this problem is a function. In principle, the equation must hold
for every admissible test function $v\in\mathcal{V}_0$. A computer, however, cannot
search an infinite-dimensional collection of functions directly. We therefore
begin by replacing that collection with functions built from a finite number
of known functions.
Choose admissible functions
$$
\psi_1,\psi_2,\ldots,\psi_n\in\mathcal{V}_0.
$$
For any coefficients $V_1,\ldots,V_n$, the linear combination
$$
v_n(x)
=
V_1\psi_1(x)+\cdots+V_n\psi_n(x)
=
\sum_{i=1}^{n}V_i\psi_i(x)
$$ {#eq-c5-test-linear-combination}
also belongs to $\mathcal{V}_0$, because each $\psi_i$ satisfies the homogeneous
essential boundary condition. Instead of allowing the test function to range
over all of $\mathcal{V}_0$, we can require the variational equation only for
functions of the form @eq-c5-test-linear-combination.
The same known functions are used to approximate the displacement:
$$
u_n(x)
=
\sum_{j=1}^{n}U_j\psi_j(x),
$$ {#eq-c5-trial-linear-combination}
where the coefficients $U_j$ are unknown. The unknown function has therefore
been replaced by $n$ unknown numbers.
The collection of all linear combinations of
$\psi_1,\ldots,\psi_n$ is a vector space contained in $\mathcal{V}_0$. In the
terminology of Chapter 2,
$$
\mathcal{V}_n
=
\operatorname{span}\{\psi_1,\ldots,\psi_n\}
\subset\mathcal{V}_0.
$$ {#eq-c5-galerkin-space}
The word *span* means precisely the set of all linear combinations in
@eq-c5-test-linear-combination. If the $\psi_i$ are linearly independent, they
form a basis of $\mathcal{V}_n$, and $\dim \mathcal{V}_n=n$. Thus *finite dimensional* means that
every function in the approximation space is determined by finitely many
coefficients.
The Galerkin approximation restricts both the trial and test functions to this
space:
$$
\text{find }u_n\in \mathcal{V}_n\text{ such that}\qquad
a(u_n,v_n)=\ell(v_n)
\qquad\forall v_n\in \mathcal{V}_n.
$$ {#eq-c5-galerkin-problem}
Substituting @eq-c5-test-linear-combination and
@eq-c5-trial-linear-combination gives
$$
\sum_{i=1}^{n}
V_i
\left[
\sum_{j=1}^{n}a(\psi_j,\psi_i)U_j-\ell(\psi_i)
\right]
=0.
$$
The coefficients $V_i$ are arbitrary. We may therefore choose them one at a
time, or equivalently require each bracketed expression to vanish. This gives
$n$ scalar equations,
$$
\sum_{j=1}^{n}a(\psi_j,\psi_i)U_j
=
\ell(\psi_i),
\qquad
i=1,\ldots,n.
$$ {#eq-c5-galerkin-equations}
Define
$$
K_{ij}=a(\psi_j,\psi_i),
\qquad
F_i=\ell(\psi_i).
$$
Then
$$
\bK\bU=\bF.
$$ {#eq-c5-galerkin-system}
The passage from the variational problem to the matrix system is therefore not
specific to finite elements. Any finite-dimensional Galerkin space leads to a
system of this form.
What makes the finite element method different is the way the functions
$\psi_i$ are chosen. If they are global polynomials or trigonometric functions,
they are typically nonzero over most of $(0,L)$, so many pairs
$(\psi_i,\psi_j)$ interact in the bilinear form and the matrix is generally
dense. The finite element method instead constructs basis functions from small
pieces associated with a partition of the domain. Each basis function is then
nonzero only near a small number of neighboring points.
{#fig-c5-global-vs-local width="95%" fig-align="center"}
This local support is the structural idea we now develop. The basis functions
constructed below will play exactly the same role as the generic
$\psi_i$ above; to emphasize their connection to the mesh, we will denote them
by $\phi_i$.
::: {#exm-c5-why-finite-dimensional}
## Why finitely many basis functions lead to finitely many unknowns
Suppose
$\mathcal{V}_3=\operatorname{span}\{\psi_1,\psi_2,\psi_3\}$.
Every function in $\mathcal{V}_3$ then has the form
$$
v_3=V_1\psi_1+V_2\psi_2+V_3\psi_3.
$$
No matter how complicated the functions $\psi_i$ themselves may be, the unknown
function $v_3$ is completely determined once the three coefficients
$V_1$, $V_2$, and $V_3$ are known. This is the basic reason the Galerkin
approximation converts a variational problem for an unknown function into a
system for finitely many unknown numbers.
:::
::: {#exr-c5-global-galerkin}
## A Galerkin approximation with global functions
For the same bar, take
$$
\psi_1(x)=x,\qquad
\psi_2(x)=x^2,\qquad
\psi_3(x)=x^3.
$$
1. Verify that each function belongs to $\mathcal{V}_0$.
2. Write $u_3(x)=\sum_{j=1}^3U_j\psi_j(x)$.
3. Write the three Galerkin equations in terms of
$a(\psi_j,\psi_i)$ and $\ell(\psi_i)$.
4. Explain why the approximation is finite dimensional even though the basis
functions are not finite element basis functions.
:::
## Mesh, elements, nodes, and connectivity {#sec-c5-mesh}
The generic Galerkin construction tells us that we need a finite collection of
basis functions, but it does not tell us how those functions should be chosen.
The finite element method makes that choice by first introducing a geometric
partition of the domain. The basis functions will then be built from simple
polynomials defined on the pieces of that partition.
For the bar occupying $[0,L]$, choose ordered points
$$
0=x_1<x_2<\cdots<x_{N_n}=L.
$$ {#eq-c5-mesh-nodes}
These points are the **nodes**. For two-node linear elements,
$$
N_n=N_e+1,
$$
where $N_e$ is the number of elements and $N_n$ is the number of nodes.
Consecutive nodes define the **elements**
$$
\mathcal{K}_e=(x_e,x_{e+1}),
\qquad
e=1,\ldots,N_e.
$$ {#eq-c5-elements}
The collection of elements
$$
\Th=\{\mathcal{K}_1,\mathcal{K}_2,\ldots,\mathcal{K}_{N_e}\}
$$
is the **finite element mesh**, or simply the **mesh**.
{#fig-c5-mesh-partition width="90%" fig-align="center"}
The mesh has two roles. Geometrically, it partitions the bar into smaller
intervals. Computationally, it gives us places on which simple local
approximations can be defined and then connected. We will construct a
polynomial approximation separately on each element, but neighboring element
approximations will share information at their common node. This combination
of local approximation and shared degrees of freedom is what allows a global
finite element function to be assembled from small pieces.
The length of element $e$ is
$$
h_e=x_{e+1}-x_e,
$$
and the mesh size is
$$
h=\max_{1\le e\le N_e}h_e.
$$
A mesh is **uniform** if all element lengths are equal. It is **nonuniform** if
they are not. Nonuniform meshes are useful when smaller elements are desired
in selected parts of the domain, for example near a geometric feature or a
region where the solution varies rapidly. The convergence consequences of
mesh refinement are studied in Chapter 7; here we use the mesh to construct
the approximation itself.
The open element interiors are disjoint, while their closures cover the bar:
$$
[0,L]
=
\bigcup_{e=1}^{N_e}[x_e,x_{e+1}].
$$
This decomposition will later allow an integral over the whole bar to be
written as a sum of element integrals.
### Local and global numbering
A node can be described in two different ways. The **global node number**
identifies a point in the complete mesh. The **local node number** identifies
the position of that point within one element. Local numbers restart on each
element.
For a two-node linear element, local node 1 is the left endpoint and local
node 2 is the right endpoint. The **connectivity map**
$$
I^e:\{1,2\}\longrightarrow\{1,\ldots,N_n\}
$$
records which global nodes belong to element $e$. With the consecutive global
numbering used here,
$$
I^e(1)=e,
\qquad
I^e(2)=e+1.
$$ {#eq-c5-linear-connectivity}
{#fig-c5-local-global-numbering width="95%" fig-align="center"}
For example, the point $x_2$ is local node 2 of $\mathcal{K}_1$ and local node 1 of
$\mathcal{K}_2$, but it is only one global node. This distinction will matter during
assembly: both neighboring elements will contribute to equations associated
with that shared global node.
For a global node $i$ that belongs to element $e$, the inverse
$(I^e)^{-1}(i)$ gives its local number on that element. Thus
$$
(I^1)^{-1}(2)=2,
\qquad
(I^2)^{-1}(2)=1.
$$ {#eq-c5-inverse-connectivity}
The element label on $I^e$ is necessary because the same global node can have
a different local number in each incident element.
For a global node $i$, define the set of elements that contain that node by
$$
\mathcal{E}(i)
=
\{e\in\{1,\ldots,N_e\}:x_i\in\overline{\mathcal{K}}_e\}.
$$ {#eq-c5-incident-elements}
On the three-element mesh in @fig-c5-local-global-numbering,
$\mathcal{E}(1)=\{1\}$, $\mathcal{E}(2)=\{1,2\}$, and
$\mathcal{E}(3)=\{2,3\}$. This notation will be useful when element
contributions are collected at a global node.
::: {.remark}
## Numbering in theory and in Python
The mathematical development in this book numbers elements, nodes, basis
functions, and matrix entries starting from 1. Python arrays and the degree-of-freedom
maps used by FEniCSx start from 0. Thus mathematical node $i$ is stored at Python
index `i - 1`, and mathematical element $e$ is stored at Python index `e - 1`.
This translation is made only when the mathematical objects are represented in
code; the two conventions should not be mixed within a derivation.
:::
It is also useful to distinguish a **node**, a **degree of freedom**, and an
**unknown coefficient**. A node is a geometric point in the mesh. A degree of
freedom is a quantity used to determine the finite element function. For the
one-dimensional nodal elements considered first, the degree of freedom is the
displacement value at a node. The corresponding coefficient is the number that
will be determined by the algebraic system. In this simple setting these three
ideas are closely aligned, but they are not identical concepts and need not
remain so for other finite element families.
At this point the mesh has been defined, but no finite element basis functions
have yet been chosen. We now use the mesh to make the generic Galerkin basis
$\psi_i$ from @sec-c5-galerkin local.
## Finite element basis functions on the mesh {#sec-c5-basis}
The functions $\psi_i$ in the Galerkin discussion were deliberately generic.
They only had to belong to the admissible function space and be linearly
independent. For a finite element approximation, we now choose those basis
functions in a particular mesh-based way. To distinguish this construction,
we write the finite element basis as
$$
\phi_1,\phi_2,\ldots,\phi_{N_n}.
$$
Thus the symbols have changed, but their role has not:
the $\phi_i$ are the basis functions used in the Galerkin expansion. The new
feature is that each $\phi_i$ will be assembled from polynomial pieces defined
on the elements of $\Th$.
For the nodal $H^1$-conforming Lagrange elements used in this chapter, three
requirements guide the construction.
1. **Polynomial on each element.** On every element $\mathcal{K}_e$, the restriction of
a basis function belongs to a polynomial space
$\mathbb P_{p_e}(\mathcal{K}_e)$.
2. **Continuous across shared endpoints.** The polynomial pieces from
neighboring elements agree at their common node. Their derivatives are not
required to agree.
3. **Nodal interpolation.** The global basis function associated with node
$x_i$ satisfies
$$
\phi_i(x_j)=\delta_{ij}.
$$ {#eq-c5-global-nodal-property}
The first requirement makes the approximation simple on each element. The
second makes the global function conform to the $H^1$ space required by the
axial-bar variational problem. The third gives the coefficients a direct nodal
interpretation.
These requirements describe the Lagrange family used here; they are not the
most general definition of a finite element. We will return to the general
definition near the end of the chapter after the concrete construction is
familiar.
### Linear basis functions
Begin with $p_e=1$ on every element. Consider again the three-element mesh
$$
\mathcal{K}_1=[x_1,x_2],
\qquad
\mathcal{K}_2=[x_2,x_3],
\qquad
\mathcal{K}_3=[x_3,x_4].
$$
We want one global basis function for each global node. Consider
$\phi_2$, associated with the interior node $x_2$. The nodal interpolation
property requires
$$
\phi_2(x_1)=0,\qquad
\phi_2(x_2)=1,\qquad
\phi_2(x_3)=0,\qquad
\phi_2(x_4)=0.
$$
Now apply the elementwise polynomial requirement. On $\mathcal{K}_1$, $\phi_2$ must be
the linear polynomial that rises from zero at $x_1$ to one at $x_2$:
$$
\phi_2(x)
=
\frac{x-x_1}{x_2-x_1},
\qquad
x\in[x_1,x_2].
$$
On $\mathcal{K}_2$, it must fall from one at $x_2$ to zero at $x_3$:
$$
\phi_2(x)
=
\frac{x_3-x}{x_3-x_2},
\qquad
x\in[x_2,x_3].
$$
On $\mathcal{K}_3$, the required nodal values at both endpoints are zero. The only
linear polynomial with those endpoint values is the zero polynomial.
Consequently,
$$
\phi_2(x)
=
\begin{cases}
\dfrac{x-x_1}{x_2-x_1},
&x\in[x_1,x_2],\\[1em]
\dfrac{x_3-x}{x_3-x_2},
&x\in[x_2,x_3],\\[1em]
0,
&\text{otherwise}.
\end{cases}
$$ {#eq-c5-global-phi2}
The remaining global basis functions are obtained by applying the same three
requirements at the other nodes.
{#fig-c5-global-linear-basis width="95%" fig-align="center"}
This construction explains where **local support** comes from. On an element
that does not contain the node associated with $\phi_i$, the function is zero
at both endpoints. Since its restriction is linear, it must vanish throughout
that element. An interior linear basis function therefore has support on only
the two elements adjacent to its node; an endpoint basis function has support
on only one element.
More precisely, the support of $\phi_i$ is the closure of the set on which it
is nonzero. For the linear nodal basis,
$$
\operatorname{supp}\phi_i
=
\overline{\{x\in[0,L]:\phi_i(x)\ne0\}}
=
\bigcup_{e\in\mathcal{E}(i)}\overline{\mathcal{K}}_e.
$$ {#eq-c5-basis-support}
The same construction can be viewed locally. On one element
$\mathcal{K}_e=[x_e,x_{e+1}]$, denote the coordinates of its local nodes by
$$
x_1^e=x_e,
\qquad
x_2^e=x_{e+1}.
$$
Define two local shape functions
$\phi_1^e$ and $\phi_2^e$ by
$$
\phi_1^e(x_e)=1,\quad
\phi_1^e(x_{e+1})=0,
\qquad
\phi_2^e(x_e)=0,\quad
\phi_2^e(x_{e+1})=1.
$$
Since the local functions are linear,
$$
\phi_1^e(x)=\frac{x_{e+1}-x}{h_e},
\qquad
\phi_2^e(x)=\frac{x-x_e}{h_e}.
$$ {#eq-c5-linear-shape-functions}
They satisfy
$$
\phi_1^e(x)+\phi_2^e(x)=1.
$$
This partition-of-unity relation means that equal nodal values reproduce a
constant field exactly.
Let $U_a^e$ denote the displacement degree of freedom at local node $a$ of
element $e$:
$$
U_a^e=u_h(x_a^e),
\qquad a=1,2.
$$
The restriction of the finite element displacement to $\mathcal{K}_e$ is
$$
\uh|_{\mathcal{K}_e}
=
U_1^e\phi_1^e+U_2^e\phi_2^e.
$$ {#eq-c5-linear-interpolation}
{#fig-c5-linear-element-interpolation width="88%" fig-align="center"}
The interpolation property makes the meaning of the coefficients immediate:
evaluating @eq-c5-linear-interpolation at the left and right endpoints returns
$U_1^e$ and $U_2^e$. The connectivity map relates the local and global
coefficients:
$$
U_a^e=U_{I^e(a)},
\qquad a=1,2.
$$ {#eq-c5-local-global-dof}
For the consecutive numbering in @eq-c5-linear-connectivity,
$U_1^e=U_e$ and $U_2^e=U_{e+1}$. Differentiating the element interpolation
therefore gives
$$
\varepsilon_h|_{\mathcal{K}_e}
=
\uh'|_{\mathcal{K}_e}
=
\frac{U_2^e-U_1^e}{h_e}
=
\frac{U_{e+1}-U_e}{h_e}.
$$ {#eq-c5-linear-strain}
Thus the two-node bar element has a linear displacement field and a constant
strain field.
For later calculations, collect the element coefficients and shape functions
as
$$
\bU^e=
\begin{bmatrix}
U_1^e\\
U_2^e
\end{bmatrix},
\qquad
\bN^e=
\begin{bmatrix}
\phi_1^e&\phi_2^e
\end{bmatrix}.
$$
Then
$$
\uh|_{\mathcal{K}_e}
=
\bN^e\bU^e,
\qquad
\varepsilon_h|_{\mathcal{K}_e}
=
\bB^e\bU^e,
\qquad
\bB^e
=
\frac{\dd\bN^e}{\dd x}
=
\begin{bmatrix}
-1/h_e&1/h_e
\end{bmatrix}.
$$ {#eq-c5-linear-NB}
### From local shape functions to the global space
The local shape functions are not a second approximation. They are the pieces
from which the global Galerkin basis is assembled. If local node $a$ of element
$e$ corresponds to global node $I^e(a)$, then
$$
\phi_{I^e(a)}|_{\mathcal{K}_e}
=
\phi_a^e.
$$ {#eq-c5-local-global-relation}
For example, the global basis function associated with an interior node $x_i$
is $\phi_2^{\,i-1}$ on the element to the left, $\phi_1^{\,i}$ on the element
to the right, and zero elsewhere. Shared nodal values join the elementwise
pieces into one continuous global function.
The resulting piecewise-linear finite element space is
$$
\mathcal{V}_h
=
\left\{
v_h\in C([0,L]):
v_h|_{\mathcal{K}_e}\in\mathbb P_1(\mathcal{K}_e)
\text{ for }e=1,\ldots,N_e
\right\}.
$$ {#eq-c5-linear-space}
Its global nodal basis is
$\{\phi_1,\phi_2,\ldots,\phi_{N_n}\}$. Every $u_h\in\mathcal{V}_h$ has a unique expansion
$$
u_h(x)
=
\sum_{j=1}^{N_n}U_j\phi_j(x).
$$ {#eq-c5-global-expansion}
At node $x_i$,
$$
u_h(x_i)
=
\sum_{j=1}^{N_n}U_j\phi_j(x_i)
=
U_i.
$$ {#eq-c5-nodal-coefficient}
The global coefficients are therefore the nodal displacement values.
Connectivity also tells us how to recover the element description from the
global one. Restrict @eq-c5-global-expansion to one element $\mathcal{K}_e$.
By @eq-c5-basis-support, every basis function whose node does not belong to
$\mathcal{K}_e$ vanishes on that element, so only the terms with
$j=I^e(1)$ and $j=I^e(2)$ remain. Applying @eq-c5-local-global-relation to
those two terms gives
$$
u_h|_{\mathcal{K}_e}
=
\sum_{j=1}^{N_n}U_j\,\phi_j|_{\mathcal{K}_e}
=
U_{I^e(1)}\phi_1^e+U_{I^e(2)}\phi_2^e
=
\bN^e\bU^e,
\qquad
\bU^e=
\begin{bmatrix}
U_{I^e(1)}\\
U_{I^e(2)}
\end{bmatrix}.
$$ {#eq-c5-global-restriction}
This is the element field @eq-c5-linear-NB, now obtained from the global
expansion instead of being written directly on the element, and it confirms
the coefficient relation @eq-c5-local-global-dof. In a program, forming
$\bU^e$ from $\bU$ is called a **gather** operation: connectivity selects the
global coefficients that belong to element $e$. Assembly, developed in
@sec-c5-assembly, performs the reverse operation and adds element
contributions to the global positions selected by the same connectivity.
Because functions in $\mathcal{V}_h$ are continuous and have square-integrable
piecewise derivatives,
$$
\mathcal{V}_h\subset H^1(0,L).
$$
This is the conformity required by the axial-bar variational problem.
::: {#exm-c5-global-function}
## Constructing a finite element function
On the three-element mesh
$x_1=0$, $x_2=1$, $x_3=2$, $x_4=3$, assign
$$
U_1=0,\qquad
U_2=1,\qquad
U_3=-\frac12,\qquad
U_4=\frac12.
$$
Then
$$
u_h
=
U_1\phi_1+U_2\phi_2+U_3\phi_3+U_4\phi_4.
$$
On $\mathcal{K}_2=[1,2]$, only $\phi_2$ and $\phi_3$ are nonzero. Hence
$$
u_h|_{\mathcal{K}_2}
=
U_2\phi_2|_{\mathcal{K}_2}+U_3\phi_3|_{\mathcal{K}_2}
=
U_2\phi_1^2+U_3\phi_2^2.
$$
The global expansion and the local interpolation are two descriptions of the
same finite element function.
:::
::: {#exr-c5-construct-basis}
## Constructing a global basis function
For the nonuniform mesh
$x_1=0$, $x_2=1$, $x_3=3$, and $x_4=4$:
1. construct $\phi_3$ from elementwise linearity, continuity, and the nodal
interpolation property;
2. state its support; and
3. verify its nodal values.
:::
::: {#exr-c5-affine-reproduction}
## Exact reproduction of an affine displacement
Let $p(x)=c_0+c_1x$ with constants $c_0$ and $c_1$, and set
$U_1^e=p(x_e)$ and $U_2^e=p(x_{e+1})$.
1. Using @eq-c5-linear-shape-functions, show that
$$
p(x)=U_1^e\phi_1^e(x)+U_2^e\phi_2^e(x)
\qquad\text{for every }x\in\mathcal{K}_e.
$$
2. Use @eq-c5-linear-strain to show that the element strain equals $c_1$.
Thus a linear bar element reproduces every constant strain exactly.
3. Take $p(x)=x^2$ on $\mathcal{K}_e=[0,h_e]$ with the same nodal values.
Compute the largest value of $|p(x)-U_1^e\phi_1^e(x)-U_2^e\phi_2^e(x)|$ on
the element and state where it occurs.
:::
## From the finite element space to element equations {#sec-c5-discrete}
The approximation space has now been constructed. We next return to the
variational problem and use exactly the same Galerkin idea as in
@sec-c5-galerkin, but with the mesh-based basis
$\phi_i$ in place of the generic basis $\psi_i$.
For the homogeneous displacement condition at $x=0$, define
$$
\mathcal{V}_{h,0}
=
\mathcal{V}_h\cap\mathcal{V}_0
=
\{v_h\in\mathcal{V}_h:v_h(0)=0\}.
$$
The finite element problem is
$$
\text{find }u_h\in\mathcal{V}_{h,0}\text{ such that}\qquad
a(u_h,v_h)=\ell(v_h)
\qquad
\forall v_h\in\mathcal{V}_{h,0}.
$$ {#eq-c5-discrete-variational}
Write $u_h=\sum_{j=1}^{N_n}U_j\phi_j$. The coefficient $U_1$ is prescribed by
the essential boundary condition. Testing with each free nodal basis function
gives
$$
\sum_{j=1}^{N_n}K_{ij}U_j=F_i,
$$
with
$$
K_{ij}
=
a(\phi_j,\phi_i)
=
\int_0^L AE\,\phi_j'\phi_i'\,\dd x,
\qquad
F_i
=
\int_0^L f\phi_i\,\dd x+\overline P\,\phi_i(L).
$$ {#eq-c5-global-entries}
The formula for $K_{ij}$ also explains the sparsity of the global matrix. If
the supports of $\phi_i$ and $\phi_j$ do not overlap except on a set of measure
zero, then their derivatives are never simultaneously nonzero over an interval,
and
$$
K_{ij}=0.
$$ {#eq-c5-support-sparsity}
For linear elements, an interior basis function overlaps only with itself and
its immediate neighbors.
{#fig-c5-basis-overlap-sparsity width="60%" fig-align="center"}
@eq-c5-global-entries already defines the
discrete problem. It would be possible to evaluate each entry by integrating
over the whole interval. Doing so, however, would ignore the local structure
that motivated the finite element basis.
The mesh partitions the interval, so
$$
K_{ij}
=
\sum_{e=1}^{N_e}
\int_{\mathcal{K}_e}
AE\,\phi_j'\phi_i'\,\dd x.
$$ {#eq-c5-global-entry-sum}
Only elements lying in the common support of $\phi_i$ and $\phi_j$ contribute.
On such an element, suppose
$i=I^e(a)$ and $j=I^e(b)$. The local-global relation
@eq-c5-local-global-relation gives
$$
\phi_i|_{\mathcal{K}_e}=\phi_a^e,
\qquad
\phi_j|_{\mathcal{K}_e}=\phi_b^e.
$$
The corresponding contribution can therefore be computed entirely with local
shape functions:
$$
K_{ab}^e
=
\int_{\mathcal{K}_e}
AE
\frac{\dd\phi_b^e}{\dd x}
\frac{\dd\phi_a^e}{\dd x}
\,\dd x.
$$ {#eq-c5-element-stiffness-entry}
It remains to identify which elements contribute and which local entry each
one supplies. An element lies in the common support of $\phi_i$ and $\phi_j$
exactly when it contains both nodes, that is, when
$e\in\mathcal{E}(i)\cap\mathcal{E}(j)$. On such an element, the inverse
connectivity map @eq-c5-inverse-connectivity gives the local numbers
$a=(I^e)^{-1}(i)$ and $b=(I^e)^{-1}(j)$. Collecting the contributions of all such elements
in @eq-c5-global-entry-sum gives
$$
K_{ij}
=
\sum_{e\in\mathcal{E}(i)\cap\mathcal{E}(j)}
K_{(I^e)^{-1}(i),(I^e)^{-1}(j)}^e.
$$ {#eq-c5-assembled-stiffness-entry}
If $\mathcal{E}(i)\cap\mathcal{E}(j)$ is empty, the sum has no terms and
$K_{ij}=0$, in agreement with @eq-c5-support-sparsity.
The distributed part of a global load entry can be decomposed in the same way:
$$
F_i^{\mathrm{dist}}
=
\int_0^L f\phi_i\,\dd x
=
\sum_{e=1}^{N_e}\int_{\mathcal{K}_e}f\phi_i\,\dd x.
$$ {#eq-c5-global-load-sum}
For each element $e$, define its local load entries by
$$
F_a^e
=
\int_{\mathcal{K}_e}f\phi_a^e\,\dd x,
\qquad a=1,2.
$$ {#eq-c5-element-load-entry}
The argument is the same as for the stiffness entry, now with one global index
instead of two. If $e\notin\mathcal{E}(i)$, then $\phi_i$ vanishes on
$\mathcal{K}_e$, so that element contributes zero. If $e\in\mathcal{E}(i)$,
the local number of global node $i$ is $(I^e)^{-1}(i)$, and
$\phi_i|_{\mathcal{K}_e}=\phi_{(I^e)^{-1}(i)}^e$. The global distributed-load entry is
therefore
$$
F_i^{\mathrm{dist}}
=
\sum_{e\in\mathcal{E}(i)} F_{(I^e)^{-1}(i)}^e.
$$ {#eq-c5-assembled-load-entry}
The complete global entry also contains the prescribed end load:
$$
F_i=F_i^{\mathrm{dist}}+\overline P\phi_i(L).
$$ {#eq-c5-complete-load-entry}
Collecting the local entries gives
$$
\bK^e
=
\int_{\mathcal{K}_e}
(\bB^e)^T AE\,\bB^e\,\dd x,
\qquad
\bF^e
=
\int_{\mathcal{K}_e}
(\bN^e)^T f\,\dd x.
$$ {#eq-c5-element-matrix-vector}
For constant $A$, $E$, and $f$ on a linear element,
$$
\bK^e
=
\frac{A_eE_e}{h_e}
\begin{bmatrix}
1&-1\\
-1&1
\end{bmatrix},
\qquad
\bF^e
=
\frac{f_eh_e}{2}
\begin{bmatrix}
1\\
1
\end{bmatrix}.
$$ {#eq-c5-linear-element-matrix-vector}
{#fig-c5-typical-element width="92%" fig-align="center"}
The picture is worth pausing over. The geometry of one element
$\mathcal{K}_e$, the two nodal displacement values collected in $\bU^e$, the
interpolation formula for $u_h|_{\mathcal{K}_e}$, and the formulas for
$\bK^e$ and $\bF^e$ belong together. In implementation, one element at a time
supplies exactly these local objects, which are then sent to the global system
through assembly.
The load vector is called a **consistent load vector** because it is obtained
from the same variational equation and basis functions used for the
displacement approximation. It is not an arbitrary replacement of the
distributed load by nodal point forces.
For example, if the distributed load varies linearly from $f_e^-$ to $f_e^+$
across the element, then
$$
f(x)=f_e^-\phi_1^e(x)+f_e^+\phi_2^e(x),
$$
and direct integration gives
$$
\bF^e
=
\frac{h_e}{6}
\begin{bmatrix}
2f_e^-+f_e^+\\
f_e^-+2f_e^+
\end{bmatrix}.
$$ {#eq-c5-linear-varying-load}
The entries are generally unequal, but their sum is
$\int_{\mathcal{K}_e}f\,\dd x$. Thus the element vector preserves the total distributed
load.
The natural boundary term is already part of the global load functional. Since
$\phi_i(L)=\delta_{i,N_n}$, an applied force $\overline P$ at $x=L$ contributes
only to the right-end global load entry:
$$
F_{N_n}\mathrel{+}=\overline P.
$$
This is the discrete counterpart of the way a natural boundary condition enters
the continuous variational problem.
::: {#exm-c5-end-force}
## A natural boundary load at the right endpoint
Suppose the bar carries no distributed load, but the right endpoint is subject
to the axial force $\overline P$. Then every element load vector from the
distributed load is zero, yet the global load vector is not. The finite element
equations acquire the single contribution
$F_{N_n}\mathrel{+}=\overline P$ at the global degree of freedom associated with
$x=L$. This is why a natural boundary condition appears on the right-hand side
rather than in the definition of the admissible displacement space.
:::
::: {#rem-c5-direct-stiffness}
## Relation to the direct stiffness method
For a two-node linear axial bar, the direct stiffness method and the finite
element method give the same element stiffness matrix. The direct stiffness
method starts from an element force--displacement relation. The finite element
derivation starts from the variational problem and the finite-dimensional
basis. The latter construction extends without changing its basic logic to the
other boundary-value problems considered in this book.
:::
::: {#exr-c5-consistent-load}
## Consistent loading
On $\mathcal{K}_e=[0,h_e]$, let $f(x)=f_0+cx$.
1. Evaluate $\bF^e$ exactly.
2. Verify that the sum of its entries equals
$\int_{\mathcal{K}_e}f(x)\,\dd x$.
3. Compare the result with simply dividing the total load equally between the
two nodes.
:::
::: {#exr-c5-foundation-matrix}
## A bar attached to distributed axial springs
Suppose the bar is attached to distributed axial springs that exert the
resisting force $-c(x)u(x)$ per unit length, where $c(x)\ge 0$. The variational
problem then contains the additional term
$$
\int_0^L c\,u v\,\dd x.
$$
For a linear element $\mathcal K_e$ on which $c(x)=c_e$ is constant:
1. show that this term contributes the element matrix
$$
C_{ab}^e
=
\int_{\mathcal K_e}c_e\phi_a^e\phi_b^e\,\dd x,
\qquad a,b\in\{1,2\};
$$
2. evaluate the integrals to obtain
$$
\bC^e
=
\frac{c_eh_e}{6}
\begin{bmatrix}
2&1\\
1&2
\end{bmatrix};
$$
3. write the complete element matrix for the bar as
$\bK^e+\bC^e$; and
4. explain why $\bK^e[1,1]^T=\bm{0}$ but
$\bC^e[1,1]^T\ne\bm{0}$ when $c_e>0$.
The same matrix also appears in scalar reaction--diffusion equations
[@flaherty2000, Sec. 1.3].
:::
## Assembly, boundary conditions, and the computed solution {#sec-c5-assembly}
The element calculations produce small matrices and vectors in local numbering.
The global finite element equations require contributions in global numbering.
The connectivity map provides exactly the information needed to move from one
description to the other.
For each element $e$,
$$
F_{I^e(a)}\mathrel{+}=F_a^e,
\qquad
K_{I^e(a),I^e(b)}\mathrel{+}=K_{ab}^e.
$$ {#eq-c5-assembly-rule}
Consider again the three-element mesh with
$$
I^1=(1,2),
\qquad
I^2=(2,3),
\qquad
I^3=(3,4).
$$
Element 1 contributes only to global rows and columns 1 and 2. Element 2
contributes to rows and columns 2 and 3, and element 3 contributes to rows and
columns 3 and 4. At a shared global node, contributions from the neighboring
elements are added.
For example, the connectivity gives
$$
F_2^{\mathrm{dist}}=F_2^1+F_1^2,
\qquad
K_{22}=K_{22}^1+K_{11}^2,
\qquad
K_{23}=K_{12}^2,
\qquad
K_{14}=0.
$$ {#eq-c5-assembly-entry-examples}
The first two entries receive contributions from both elements incident on
node 2. The entry $K_{23}$ receives one contribution because nodes 2 and 3
share only element 2, while $K_{14}$ is zero because the supports of
$\phi_1$ and $\phi_4$ do not overlap. Each of these statements is an instance
of @eq-c5-assembled-load-entry or @eq-c5-assembled-stiffness-entry. For
$K_{22}$, the sum runs over $\mathcal{E}(2)=\{1,2\}$ with $J^1(2)=2$ and
$J^2(2)=1$; for $K_{14}$, the set $\mathcal{E}(1)\cap\mathcal{E}(4)$ is
empty.
If
$$
\bK^e=
\begin{bmatrix}
K_{11}^e&K_{12}^e\\
K_{21}^e&K_{22}^e
\end{bmatrix},
$$
the assembled matrix is
$$
\bK
=
\begin{bmatrix}
K_{11}^1
&
K_{12}^1
&
0
&
0\\
K_{21}^1
&
K_{22}^1+K_{11}^2
&
K_{12}^2
&
0\\
0
&
K_{21}^2
&
K_{22}^2+K_{11}^3
&
K_{12}^3\\
0
&
0
&
K_{21}^3
&
K_{22}^3
\end{bmatrix}.
$$ {#eq-c5-three-element-K}
The load vector can be built in the same way. Write the contribution of each
element as a vector of global length that is zero except in the two rows
selected by its connectivity:
$$
\begin{bmatrix}
F_1^1\\
F_2^1\\
0\\
0
\end{bmatrix}
+
\begin{bmatrix}
0\\
F_1^2\\
F_2^2\\
0
\end{bmatrix}
+
\begin{bmatrix}
0\\
0\\
F_1^3\\
F_2^3
\end{bmatrix}.
$$
Adding these three vectors gives the distributed-load vector
$$
\bF
=
\begin{bmatrix}
F_1^1\\
F_2^1+F_1^2\\
F_2^2+F_1^3\\
F_2^3
\end{bmatrix}
$$ {#eq-c5-three-element-F}
before any additional natural boundary contribution is added.
::: {.remark}
## Two equivalent views of assembly
@eq-c5-assembled-load-entry fixes a global node $i$ and collects the local
entry from every element in $\mathcal{E}(i)$; $(I^e)^{-1}$ identifies the
required local node on each element. In an implementation, one normally
uses the equivalent forward operation: compute $\bF^e$, loop over its local
entries $F_a^e$, and add each entry to the global position $I^e(a)$ according
to @eq-c5-assembly-rule. Both views produce the global vector in
@eq-c5-three-element-F. The stiffness matrix is treated in the same way:
@eq-c5-assembled-stiffness-entry gathers the contributions to one global entry
$K_{ij}$, whereas the forward rule loops over the local pairs $(a,b)$ of each
element and adds $K_{ab}^e$ to $K_{I^e(a),I^e(b)}$. The forward operation is
called **scatter-add** because each element scatters its local entries into
the global arrays and adds them to what is already there.
:::
{#fig-c5-assembly width="100%" fig-align="center"}
For equal element length $h$ and constant $AE$,
$$
\bK
=
\frac{AE}{h}
\begin{bmatrix}
1&-1&0&0\\
-1&2&-1&0\\
0&-1&2&-1\\
0&0&-1&1
\end{bmatrix}.
$$ {#eq-c5-unconstrained-K}
Before an essential displacement condition is imposed,
$[1,1,1,1]^T$ lies in the nullspace of this matrix. That vector represents a
rigid translation: all nodal displacements change by the same amount, so the
strain in every element remains zero.
### Essential boundary conditions and reactions
Suppose some nodal displacements are prescribed. Those coefficients are already
known and should not be solved for again. It is useful to see the algebra
before introducing block notation. If $U_1$ is prescribed, a free global
equation of the form
$$
K_{i1}U_1+K_{i2}U_2+\cdots+K_{i,N_n}U_{N_n}=F_i
$$
is rewritten by moving the known term to the right-hand side:
$$
K_{i2}U_2+\cdots+K_{i,N_n}U_{N_n}
=
F_i-K_{i1}U_1.
$$
A nonzero prescribed displacement therefore changes the right-hand side of the
free equations. The same operation can be written compactly by partitioning the
global vector into free and prescribed values,
$$
\bU=
\begin{bmatrix}
\bU_f\\
\bU_p
\end{bmatrix},
$$
and partitioning the assembled equations consistently:
$$
\begin{bmatrix}
\bK_{ff}&\bK_{fp}\\
\bK_{pf}&\bK_{pp}
\end{bmatrix}
\begin{bmatrix}
\bU_f\\
\bU_p
\end{bmatrix}
=
\begin{bmatrix}
\bF_f\\
\bF_p
\end{bmatrix}.
$$ {#eq-c5-partitioned-system}
The equations associated with the free degrees of freedom are
$$
\bK_{ff}\bU_f
=
\bF_f-\bK_{fp}\bU_p.
$$ {#eq-c5-reduced-system}
For a nonzero prescribed displacement, the term
$\bK_{fp}\bU_p$ is essential. It represents the contribution of the known
displacement values to the free equilibrium equations.
The equations associated with prescribed degrees of freedom are not discarded.
After the free displacements have been found, they provide the forces required
to enforce the prescribed values:
$$
\bm r_p
=
\bK_{pf}\bU_f+\bK_{pp}\bU_p-\bF_p.
$$ {#eq-c5-reactions}
In computation it is convenient to insert the prescribed values into $\bU$,
form the complete residual vector
$$
\bR=\bF-\bK\bU,
$$ {#eq-c5-residual-vector}
and recover the support reactions as $\bm r_p=-\bR_p$. This sign follows
because $\bR$ is external force minus internal force; the support supplies
the additional force needed for equilibrium. The free entries of
$\bR$ are $\bF_f-\bK_{ff}\bU_f-\bK_{fp}\bU_p$, which vanish, up to rounding
error, because they are exactly the equations solved in
@eq-c5-reduced-system. Checking them is therefore a simple test of the solve.
{#fig-c5-essential-bc width="100%" fig-align="center"}
The distinction from the continuous variational problem is worth keeping in
view. Natural boundary data contribute through the load functional. Essential
boundary data restrict the admissible functions and, after discretization,
identify coefficients whose values are prescribed.
The partition also completes the variational statement for a nonzero
prescribed displacement, which the beginning of this chapter set aside.
Suppose $u(0)=\overline u$. The finite element trial space is then
$$
\mathcal U_{h,\overline u}
=
\{w_h\in\mathcal{V}_h:w_h(0)=\overline u\},
$$
while the test space remains $\mathcal{V}_{h,0}$. As in Chapter 3,
$\mathcal U_{h,\overline u}$ is not a vector space when $\overline u\ne0$.
Since $\phi_1(0)=1$ and every other nodal basis function vanishes at $x=0$,
each $w_h\in\mathcal U_{h,\overline u}$ can be written uniquely as
$w_h=\overline u\,\phi_1+z_h$ with $z_h\in\mathcal{V}_{h,0}$. The finite
element problem is to find $u_h\in\mathcal U_{h,\overline u}$ such that
$a(u_h,v_h)=\ell(v_h)$ for every $v_h\in\mathcal{V}_{h,0}$. Substituting
$u_h=\overline u\,\phi_1+z_h$ and using bilinearity gives
$$
\text{find }z_h\in\mathcal{V}_{h,0}\text{ such that}\qquad
a(z_h,v_h)=\ell(v_h)-\overline u\,a(\phi_1,v_h)
\qquad\forall v_h\in\mathcal{V}_{h,0}.
$$ {#eq-c5-nonzero-prescribed}
The known function $\overline u\,\phi_1$ carries the prescribed value, and its
contribution moves to the right-hand side. Writing
$z_h=\sum_{j=2}^{N_n}U_j\phi_j$ and testing with $\phi_i$ for a free node $i$
gives
$$
\sum_{j=2}^{N_n}K_{ij}U_j=F_i-K_{i1}\overline u,
$$
because $a(\phi_1,\phi_i)=K_{i1}$. This is the free equation written at the
beginning of this subsection with $U_1=\overline u$, and hence
@eq-c5-reduced-system with $\bU_p=[\overline u]$.
The continuous problem @eq-model-bar-variational can be treated in the same
way. If $w$ is any function in $H^1(0,L)$ with $w(0)=\overline u$, write
$u=w+z$ with $z\in\mathcal V_0$. Then $z$ must satisfy
$a(z,v)=\ell(v)-a(w,v)$ for every $v\in\mathcal V_0$. The right-hand side is a
continuous linear functional on $\mathcal V_0$, so the Lax--Milgram theorem
@thm-lax-milgram gives a unique $z$ whenever $a$ is continuous and coercive on
$\mathcal V_0$. The resulting $u$ does not depend on the choice of $w$: the
difference of two solutions lies in $\mathcal V_0$, satisfies the homogeneous
variational equation, and is therefore zero by coercivity.
::: {#exm-c5-complete-bar}
## Complete analysis of a uniformly loaded bar
Consider a uniform bar of length $L$ and constant axial rigidity $AE$. Let the
left endpoint be fixed, the right endpoint be traction free, and let the
distributed load be the constant $f_0$.
Divide the bar into two equal elements. Then $h=L/2$ and
$$
I^1=(1,2),
\qquad
I^2=(2,3).
$$
Each element contributes
$$
\bK^e
=
\frac{AE}{h}
\begin{bmatrix}
1&-1\\
-1&1
\end{bmatrix},
\qquad
\bF^e
=
\frac{f_0h}{2}
\begin{bmatrix}
1\\
1
\end{bmatrix}.
$$
After assembly,
$$
\frac{AE}{h}
\begin{bmatrix}
1&-1&0\\
-1&2&-1\\
0&-1&1
\end{bmatrix}
\begin{bmatrix}
U_1\\
U_2\\
U_3
\end{bmatrix}
=
\frac{f_0h}{2}
\begin{bmatrix}
1\\
2\\
1
\end{bmatrix}.
$$ {#eq-c5-two-element-system}
The middle load entry contains one contribution from each neighboring element.
Since $U_1=0$, the free equations are
$$
\frac{AE}{h}
\begin{bmatrix}
2&-1\\
-1&1
\end{bmatrix}
\begin{bmatrix}
U_2\\
U_3
\end{bmatrix}
=
\frac{f_0h}{2}
\begin{bmatrix}
2\\
1
\end{bmatrix}.
$$
Solving and using $h=L/2$ gives
$$
U_2
=
\frac{3f_0L^2}{8AE},
\qquad
U_3
=
\frac{f_0L^2}{2AE}.
$$ {#eq-c5-two-element-solution}
The exact displacement is
$$
u(x)
=
\frac{f_0}{AE}
\left(
Lx-\frac{x^2}{2}
\right).
$$ {#eq-c5-uniform-load-exact}
For this particular problem, the finite element solution agrees with the exact
solution at the three nodes, but it remains linear within each element.
The corresponding element strains are
$$
\varepsilon_h^1
=
\frac{U_2-U_1}{h}
=
\frac{3f_0L}{4AE},
\qquad
\varepsilon_h^2
=
\frac{U_3-U_2}{h}
=
\frac{f_0L}{4AE}.
$$
The exact strain
$u'(x)=f_0(L-x)/(AE)$ varies linearly, whereas the finite element strain is
piecewise constant. The element stresses and axial forces are obtained from
$$
\sigma_h^e=E\varepsilon_h^e,
\qquad
N_h^e=A\sigma_h^e.
$$
Finally, the reaction at node 1 is obtained from the first row of the original
assembled system:
$$
r_1=-R_1
=
-\frac{AE}{h}U_2-\frac{f_0h}{2}
=
-f_0L.
$$
Because $\int_0^L f_0\,\dd x=f_0L$,
$$
r_1+\int_0^L f_0\,\dd x=0.
$$
The recovered reaction therefore satisfies global equilibrium.
:::
Once the global coefficient vector is known, postprocessing proceeds
element by element. Connectivity extracts the local coefficient vector
$\bU^e$ from $\bU$, and the formulas for $u_h$, strain, stress, and axial
force are evaluated on each element. Adjacent linear elements may give
different strains and stresses at a shared node. This does not violate
conformity: $\mathcal{V}_h$ requires displacement to be continuous, not its derivative.
Before accepting a result, check that the mesh and connectivity represent the
intended bar, all quantities use consistent units, prescribed and free degrees
of freedom are identified correctly, applied loads and recovered reactions
satisfy global equilibrium, and the displacement and axial-force trends are
consistent with the loading. Chapter 7 develops systematic error estimates and
convergence theory.
::: {#exr-c5-assembly-error}
## Diagnosing an assembly error
Two linear elements have
$k_1=(AE)_1/h_1$ and $k_2=(AE)_2/h_2$. A student assembles
$$
\bK
=
\begin{bmatrix}
k_1&-k_1&0\\
-k_1&k_1&-k_2\\
0&-k_2&k_2
\end{bmatrix}.
$$
Identify the missing contribution, correct the matrix, and explain which global
degree of freedom should receive the missing term.
:::
::: {#exm-c5-python-assembly}
## Assembly and mesh refinement in Python
We now repeat the uniformly loaded bar problem from
@exm-c5-complete-bar with an arbitrary number of uniform linear elements. Set
$L=AE=f_0=1$. The boundary-value problem and its exact solution are
$$
-u''=1\quad\text{in }(0,1),
\qquad
u(0)=0,
\qquad
u'(1)=0,
\qquad
u(x)=x-\frac{x^2}{2}.
$$
The following program exposes the element calculation, connectivity, assembly,
imposition of the prescribed displacement, solution of the reduced system,
and recovery of the reaction. The arrays use Python's zero-based indexing;
thus `dofs = [e, e + 1]` represents the theoretical connectivity
$I^{e+1}=(e+1,e+2)$. The two nested loops over the local indices `a` and `b`
transcribe the scatter-add rule @eq-c5-assembly-rule line by line, and the
reaction is the negative of the prescribed entry of the residual vector
@eq-c5-residual-vector, or equivalently the prescribed entry of `K @ U - F`.
```{python}
#| label: fig-c5-python-mesh-refinement
#| fig-cap: "Linear finite element solutions for a uniformly loaded bar. Mesh refinement reduces the error between nodes."
#| fig-width: 6.8
#| fig-height: 3.7
import numpy as np
import matplotlib.pyplot as plt
def solve_bar(num_elements):
x = np.linspace(0.0, 1.0, num_elements + 1)
K = np.zeros((num_elements + 1, num_elements + 1))
F = np.zeros(num_elements + 1)
for e in range(num_elements):
dofs = np.array([e, e + 1])
h = x[e + 1] - x[e]
K_e = np.array([[1.0, -1.0], [-1.0, 1.0]]) / h
F_e = 0.5 * h * np.ones(2)
for a in range(2):
F[dofs[a]] += F_e[a]
for b in range(2):
K[dofs[a], dofs[b]] += K_e[a, b]
prescribed = np.array([0])
free = np.arange(1, num_elements + 1)
U = np.zeros(num_elements + 1)
rhs = F[free] - K[np.ix_(free, prescribed)] @ U[prescribed]
U[free] = np.linalg.solve(K[np.ix_(free, free)], rhs)
R = F - K @ U
reaction = -R[prescribed]
return x, U, reaction
x_plot = np.linspace(0.0, 1.0, 501)
u_exact = x_plot - 0.5 * x_plot**2
for N in (2, 4, 8):
x, U, reaction = solve_bar(N)
u_h = np.interp(x_plot, x, U)
error = np.max(np.abs(u_h - u_exact))
plt.plot(x, U, "o-", label=fr"$N={N}$, max. error $={error:.4f}$")
plt.plot(x_plot, u_exact, "k--", linewidth=2, label="exact")
plt.xlabel("$x$")
plt.ylabel("$u(x)$")
plt.legend()
plt.tight_layout()
```
The finite element coefficients agree with the exact solution at the nodes for
this particular problem, but the linear finite element field does not reproduce
the quadratic variation between nodes. As $N$ increases, the element length
decreases and the piecewise-linear curves in @fig-c5-python-mesh-refinement
approach the exact solution. The recovered reaction is $-1$ for every mesh,
which balances the total distributed load.
The maximum error reported in the legend decreases by a factor of about four
each time the number of elements is doubled. Because the nodal values are
exact here, this error is entirely the error of linear interpolation between
the nodes. Chapter 7 explains why the nodal values are exact for this problem,
why the displacement error behaves like $h^2$, and why the strain error
decreases only like $h$.
:::
## Higher-order elements and the reference interval {#sec-c5-higher-order}
The construction above is not restricted to linear interpolation. We can keep
the same mesh-based idea and enrich the polynomial space on each element.
Consider a three-node quadratic element with local nodes
$$
x_1^e<x_2^e<x_3^e.
$$
The local basis functions belong to $\mathbb P_2(\mathcal{K}_e)$ and satisfy
$$
\phi_a^e(x_b^e)=\delta_{ab},
\qquad
a,b\in\{1,2,3\}.
$$ {#eq-c5-quadratic-nodal-property}
The three Lagrange polynomials are
$$
\begin{aligned}
\phi_1^e(x)
&=
\frac{(x-x_2^e)(x-x_3^e)}
{(x_1^e-x_2^e)(x_1^e-x_3^e)},\\
\phi_2^e(x)
&=
\frac{(x-x_1^e)(x-x_3^e)}
{(x_2^e-x_1^e)(x_2^e-x_3^e)},\\
\phi_3^e(x)
&=
\frac{(x-x_1^e)(x-x_2^e)}
{(x_3^e-x_1^e)(x_3^e-x_2^e)}.
\end{aligned}
$$ {#eq-c5-quadratic-basis}
The interpolation
$$
u_h|_{\mathcal{K}_e}
=
U_1^e\phi_1^e+
U_2^e\phi_2^e+
U_3^e\phi_3^e
$$
is quadratic, so its derivative is linear. A quadratic bar element can
therefore represent a linearly varying strain within one element.
On a mesh of quadratic elements, adjacent elements share endpoint nodes, while
each midpoint node belongs to a single element. One convenient connectivity
for three elements is
$$
I^1=(1,2,3),
\qquad
I^2=(3,4,5),
\qquad
I^3=(5,6,7).
$$ {#eq-c5-quadratic-connectivity}
{#fig-c5-linear-quadratic-connectivity width="95%" fig-align="center"}
A global quadratic basis function associated with a midpoint is supported on
one element. A basis function associated with a shared endpoint is obtained by
joining one quadratic piece from each neighboring element.
{#fig-c5-global-quadratic-basis width="95%" fig-align="center"}
The corresponding global space is
$$
\mathcal{V}_h^{(2)}
=
\left\{
v_h\in C([0,L]):
v_h|_{\mathcal{K}_e}\in\mathbb P_2(\mathcal{K}_e)
\text{ for every }e
\right\}.
$$
Its functions are continuous across element boundaries, but their derivatives
need not be continuous. Increasing the polynomial degree on a fixed mesh is
called **$p$-refinement**; subdividing the mesh so that the element size
decreases is called **$h$-refinement**.
The local matrix for a quadratic element is $3\times3$, but the assembly
principle does not change. For two elements with
$I^1=(1,2,3)$ and $I^2=(3,4,5)$, the shared global degree of freedom 3 receives
contributions from both elements. For example,
$$
K_{33}\mathrel{+}=K_{33}^1,
\qquad
K_{33}\mathrel{+}=K_{11}^2.
$$
The element has more local degrees of freedom, but connectivity still tells us
where each local entry belongs.
### Reference interval and general degree
Constructing the linear and quadratic basis functions directly on physical
elements makes their meaning clear. It also reveals a repeated calculation:
every element uses the same type of interpolation problem. In a finite element
implementation, it is more convenient to define the local basis once on a
fixed reference interval and map it to each physical element.
Let
$$
\widehat{\mathcal{K}}=[-1,1]
$$
with reference coordinate $\xi$. A coordinate map for element $e$ is a
one-to-one map
$$
x^e:\widehat{\mathcal K}\longrightarrow\mathcal K_e,
\qquad
x=x^e(\xi),
$$ {#eq-c5-abstract-element-map}
from the reference element to the physical element. The reference element and
its basis functions are the same for every element. The map $x^e$ contains the
element-specific geometry.
For the physical interval $\mathcal{K}_e=[x_e,x_{e+1}]$, first seek an affine
map $x^e(\xi)=m_e\xi+c_e$ satisfying
$x^e(-1)=x_e$ and $x^e(1)=x_{e+1}$. Solving for $m_e$ and $c_e$ gives
$$
x^e(\xi)
=
\frac{x_e+x_{e+1}}{2}
+
\frac{h_e}{2}\xi,
\qquad
J_e=\frac{\dd x}{\dd\xi}=\frac{h_e}{2}.
$$ {#eq-c5-reference-map}
Thus $x^e$ sends the two reference nodes to the two physical nodes. Its inverse
is
$$
(x^e)^{-1}(x)
=
\frac{2}{h_e}
\left(
x-\frac{x_e+x_{e+1}}{2}
\right).
$$
{#fig-c5-reference-map width="86%" fig-align="center"}
For a degree-$p$ Lagrange element, choose $p+1$ distinct reference nodes
$\xi_1,\ldots,\xi_{p+1}$. The reference basis function associated with
$\xi_a$ is
$$
\widehat\phi_a(\xi)
=
\prod_{\substack{b=1\\b\ne a}}^{p+1}
\frac{\xi-\xi_b}{\xi_a-\xi_b},
\qquad
a=1,\ldots,p+1.
$$ {#eq-c5-reference-lagrange}
By construction,
$\widehat\phi_a(\xi_b)=\delta_{ab}$. For $p=1$, with nodes $-1$ and $1$,
$$
\widehat\phi_1(\xi)=\frac{1-\xi}{2},
\qquad
\widehat\phi_2(\xi)=\frac{1+\xi}{2}.
$$
For $p=2$, with nodes $-1,0,1$,
$$
\widehat\phi_1(\xi)=\frac12\xi(\xi-1),
\qquad
\widehat\phi_2(\xi)=1-\xi^2,
\qquad
\widehat\phi_3(\xi)=\frac12\xi(\xi+1).
$$ {#eq-c5-quadratic-reference-basis}
If $U_a^e$ are the local coefficients, the finite element field is evaluated on
the reference interval as
$$
u_h(x^e(\xi))
=
\sum_{a=1}^{p+1}
U_a^e\widehat\phi_a(\xi).
$$ {#eq-c5-reference-field-interpolation}
Equivalently, the physical basis function is obtained by composition:
$$
\phi_a^e(x)
=
\widehat\phi_a((x^e)^{-1}(x)).
$$ {#eq-c5-reference-to-physical}
### Geometry and field interpolation {#sec-c5-parametric-interpolation}
The geometry map and the finite element field are two different
interpolations. In general, write the geometry map as
$$
x^e(\xi)
=
\sum_{A=1}^{n_g}
\widehat N_A^{\mathrm{geo}}(\xi)x_A^e,
$$ {#eq-c5-general-geometry-map}
where $\widehat N_A^{\mathrm{geo}}$ are geometry interpolation functions and
$x_A^e$ are physical geometry nodes. The field interpolation remains
@eq-c5-reference-field-interpolation, with basis functions
$\widehat\phi_a$ and field coefficients $U_a^e$. Keeping these two
interpolations distinct makes the standard terminology precise:
| description | relation between geometry and field interpolation |
|---|---|
| **subparametric** | the geometry interpolation has lower polynomial degree than the field interpolation |
| **isoparametric** | the geometry and field use the same interpolation functions |
| **superparametric** | the geometry interpolation has higher polynomial degree than the field interpolation |
These names compare the geometry interpolation with the approximation of the
unknown field [@jog, pp. 78--79].
{#fig-c5-parametric-types width="98%" fig-align="center"}
For an isoparametric degree-$p$ Lagrange element,
$n_g=p+1$ and $\widehat N_a^{\mathrm{geo}}=\widehat\phi_a$. Therefore,
$$
x^e(\xi)
=
\sum_{a=1}^{p+1}
\widehat\phi_a(\xi)x_a^e,
\qquad
J_e(\xi)
=
\frac{\dd x^e}{\dd\xi}.
$$ {#eq-c5-isoparametric-map}
For $p=1$, @eq-c5-isoparametric-map is exactly the affine map in
@eq-c5-reference-map. For $p=2$, the same three quadratic functions used to
interpolate $u_h$ also interpolate the geometry. If the middle geometry node
lies at the physical midpoint, the quadratic terms cancel and the map again
reduces to @eq-c5-reference-map. If the middle geometry node is not at the
midpoint, the parameterization is nonaffine and $J_e$ varies with $\xi$.
A valid orientation-preserving one-dimensional map must satisfy
$$
J_e(\xi)>0
\qquad
\text{for }-1\le\xi\le1.
$$
The chain rule gives
$$
\frac{\dd\phi_a^e}{\dd x}
=
\frac{1}{J_e}
\frac{\dd\widehat\phi_a}{\dd\xi}.
$$ {#eq-c5-derivative-transform}
The reference basis can therefore be reused on every element; element-specific
information enters through the map and connectivity.
::: {#exr-c5-parametric-classification}
## Classifying geometry and field interpolation
Consider a three-node quadratic field approximation on an interval.
1. Write the field interpolation using the three functions in
@eq-c5-quadratic-reference-basis.
2. Classify the description as subparametric when the geometry uses only the
two linear endpoint functions.
3. Write the isoparametric geometry map using all three quadratic functions.
4. Substitute
$x_2^e=(x_1^e+x_3^e)/2$ into the isoparametric map and show that it reduces
to an affine map.
5. Explain why moving $x_2^e$ away from the midpoint makes $J_e$ depend on
$\xi$, even though the physical element remains an interval.
:::
### Numerical integration on the reference interval
The same map standardizes element integration. For an integrable function
$g$,
$$
\int_{\mathcal{K}_e}g(x)\,\dd x
=
\int_{-1}^{1}g(x^e(\xi))J_e(\xi)\,\dd\xi.
$$ {#eq-c5-change-of-variables}
The element entries become
$$
K_{ab}^e
=
\int_{-1}^{1}
\frac{\dd\widehat\phi_b}{\dd\xi}
A(x^e(\xi))E(x^e(\xi))
\frac{\dd\widehat\phi_a}{\dd\xi}
\frac{1}{J_e(\xi)}
\,\dd\xi,
$$ {#eq-c5-reference-stiffness}
and
$$
F_a^e
=
\int_{-1}^{1}
\widehat\phi_a(\xi)
f(x^e(\xi))
J_e(\xi)
\,\dd\xi.
$$ {#eq-c5-reference-load}
These integrals have the same limits and use the same reference basis on every
element. For constant coefficients and simple polynomial data, some of them
can be evaluated analytically. A general finite element code, however, must
also handle spatially varying material properties, loads, and geometry. It is
therefore preferable to evaluate element integrals from values of the
integrand at a fixed set of points on the reference interval.
Gauss--Legendre quadrature provides such a rule. It approximates
$$
\int_{-1}^{1}g(\xi)\,\dd\xi
\approx
\sum_{r=1}^{n_q}w_r g(\xi_r).
$$ {#eq-c5-gauss-rule}
The first three rules are
| $n_q$ | quadrature points | weights | exact through degree |
|---:|---|---|---:|
| 1 | $0$ | $2$ | 1 |
| 2 | $-1/\sqrt{3},\,1/\sqrt{3}$ | $1,\,1$ | 3 |
| 3 | $-\sqrt{3/5},\,0,\,\sqrt{3/5}$ | $5/9,\,8/9,\,5/9$ | 5 |
Thus
$$
K_{ab}^e
\approx
\sum_{r=1}^{n_q}
w_r
\frac{\dd\widehat\phi_b}{\dd\xi}(\xi_r)
A(x_r)E(x_r)
\frac{\dd\widehat\phi_a}{\dd\xi}(\xi_r)
\frac{1}{J_e(\xi_r)},
$$
and
$$
F_a^e
\approx
\sum_{r=1}^{n_q}
w_r
\widehat\phi_a(\xi_r)
f(x_r)
J_e(\xi_r),
\qquad
x_r=x^e(\xi_r).
$$ {#eq-c5-quadrature-entries}
For affine geometry and constant $AE$, the linear-element stiffness integrand
is constant, so one Gauss point is sufficient. For a quadratic element the
basis derivatives are linear, so the stiffness integrand is quadratic and two
Gauss points integrate it exactly. If $A$, $E$, $f$, or the Jacobian varies,
the required quadrature order depends on the complete integrand.
For a three-node quadratic element with midpoint geometry and constant data,
two-point quadrature gives
$$
\bK^e
=
\frac{AE}{3h_e}
\begin{bmatrix}
7&-8&1\\
-8&16&-8\\
1&-8&7
\end{bmatrix},
\qquad
\bF^e
=
\frac{f_eh_e}{6}
\begin{bmatrix}
1\\
4\\
1
\end{bmatrix}.
$$ {#eq-c5-quadratic-element-matrix-vector}
The load entries add to $f_eh_e$, and
$\bK^e[1,1,1]^T=\bm{0}$ because equal nodal displacements produce zero strain.
::: {#exm-c5-quadratic-exact}
## One quadratic element for the uniformly loaded bar
Return to @exm-c5-complete-bar. The exact displacement
@eq-c5-uniform-load-exact is quadratic. Use one quadratic element with
interpolation nodes at $0$, $L/2$, and $L$.
After imposing $U_1=0$, the finite element equations give
$$
U_2
=
\frac{3f_0L^2}{8AE},
\qquad
U_3
=
\frac{f_0L^2}{2AE}.
$$
These are the exact displacement values at $x=L/2$ and $x=L$. Since the exact
solution and the finite element interpolant are both quadratic and agree at all
three interpolation nodes,
$$
u_h(x)=u(x)
\qquad\forall x\in[0,L].
$$
The agreement occurs because the exact solution belongs to the selected finite
element space; it is not a general property of quadratic elements.
:::
::: {#exr-c5-quadratic-connectivity}
## Quadratic connectivity and assembly
For two quadratic elements with
$I^1=(1,2,3)$ and $I^2=(3,4,5)$:
1. state the support of each global basis function $\phi_1,\ldots,\phi_5$;
2. identify every element contribution to the global entry $K_{33}$; and
3. explain why the assembly rule has the same form as for linear elements.
:::
::: {#exr-c5-quadrature}
## Selecting a quadrature rule
Determine the minimum Gauss--Legendre rule that integrates the stiffness matrix
exactly for:
1. a linear element with affine geometry and constant $AE$;
2. a quadratic element with affine geometry and constant $AE$; and
3. a quadratic element with affine geometry and $AE$ varying linearly with
$x$.
:::
## The finite element in general form {#sec-c5-general}
The linear and quadratic elements above were described through nodal values:
given the values of a polynomial at the interpolation nodes, the polynomial is
determined uniquely. The essential idea is more general than point values.
A local degree of freedom can be any linear rule that extracts a number from a
function in the local approximation space.
For example, point evaluation at $x_e$ is the rule
$$
N_1(w)=w(x_e).
$$
Such a rule is a **linear functional** because
$N_1(\alpha w+\beta z)=\alpha N_1(w)+\beta N_1(z)$. The collection of all
linear functionals on a vector space $\mathcal P$ is called its **dual
space** and is denoted by $\mathcal P^*$. With this terminology, the
construction can be stated without assuming that the degrees of freedom are
nodal values.
::: {#def-c5-finite-element}
## Finite element
A **finite element** is a triple $(K,\mathcal P,\mathcal N)$, where
1. $K$ is an element;
2. $\mathcal P$ is a finite-dimensional vector space of functions on $K$; and
3. $\mathcal N=\{N_1,\ldots,N_m\}$ is a basis for the dual space
$\mathcal P^*$.
Here $m=\dim\mathcal P$, and the functionals
$N_a:\mathcal P\to\R$ are the local degrees of freedom.
:::
For the two-node linear interval element,
$$
K=[x_e,x_{e+1}],
\qquad
\mathcal P=\mathbb P_1(K),
$$
and
$$
N_1(w)=w(x_e),
\qquad
N_2(w)=w(x_{e+1}).
$$
The degrees of freedom must determine every function in $\mathcal P$
uniquely. Equivalently,
$$
N_a(w)=0
\quad\forall a
\qquad\Longrightarrow\qquad
w=0.
$$
This property is called **unisolvence**. The nodal basis is characterized by
$$
N_a(\phi_b)=\delta_{ab}.
$$
For the Lagrange elements used in this chapter, the degrees of freedom are
point evaluations at interpolation nodes. Other finite element families use
different functionals.
The complete one-dimensional calculation can now be stated without assuming a
particular polynomial degree. Reference basis values and derivatives are
evaluated for the chosen element family and quadrature rule. Each element then
supplies its geometry, material data, loading, and connectivity:
```text
initialize the global matrix K and vector F
for each element e:
read its coordinates and global degree-of-freedom numbers
initialize K^e and F^e
for each quadrature point:
evaluate reference basis functions and derivatives
map the point to the physical element
compute the Jacobian and physical derivatives
accumulate the contribution to K^e and F^e
assemble K^e and F^e into K and F using connectivity
identify free and prescribed degrees of freedom
solve the reduced system
insert prescribed values into the global coefficient vector
recover reactions from the prescribed equations
postprocess displacement, strain, stress, and axial force element by element
```
Changing the polynomial degree modifies the local basis, the number of local
degrees of freedom, the size of the element arrays, the connectivity, and
usually the quadrature requirement. The Galerkin formulation,
local-to-global assembly, treatment of prescribed values, and overall solution
procedure remain the same.
::: {#exr-c5-trace-element}
## Tracing one element through the procedure
For the first element in @exm-c5-complete-bar, state
1. its connectivity and coordinate vector;
2. its element stiffness matrix and load vector;
3. every global matrix and vector entry to which it contributes;
4. the free and prescribed degrees of freedom after assembly; and
5. the quantities extracted during postprocessing.
:::
## Chapter summary {#sec-c5-summary}
- A Galerkin approximation replaces an infinite-dimensional variational
problem by a finite-dimensional problem. The finite element method constructs
the approximation space from a mesh and local polynomial spaces.
- Nodal $H^1$-conforming Lagrange basis functions are polynomial on each
element, continuous across shared endpoints, and satisfy the nodal
interpolation property.
- Each global basis function has local support. Only basis functions with
overlapping support contribute to the same matrix entry, which produces a
sparse global system.
- The local-to-global connectivity map identifies the restriction of a global
basis function with a local basis function. It gathers local coefficients
from $\bU$ and scatter-adds element matrices and vectors into $\bK$ and
$\bF$.
- Essential boundary conditions prescribe selected coefficients and produce a
reduced system for the free degrees of freedom. Natural boundary data enter
the load vector, and reactions are the negatives of the prescribed entries
of the residual after the solution is known.
- Quadratic and higher-order Lagrange elements follow the same construction
with richer local polynomial spaces. The reference element makes the basis,
coordinate mapping, and quadrature reusable across physical elements.
- Subparametric, isoparametric, and superparametric descriptions distinguish
the polynomial degrees used for geometry and field interpolation. The
triple $(K,\mathcal P,\mathcal N)$ describes the local finite element
independently of a particular degree.
