32.4.1 Localized Support in One Dimension: Hat Functions

However, FEM takes a different approach. Here, we use functions with localized support. Localized-support functions are functions that are only defined in a very small interval and are zero everywhere else. They are therefore piecewise-defined. The most common functions for one-dimensional FEM are hat functions (see Fig. 32.3a). The name stems from the fact that they look like a very pointed hat. As you can see from Fig. 32.3a the domain is first subdivided into intervals that are not necessarily of the same size. In the interval j a piecewise-linear hat function ψ_j is defined such that it reaches the center value ${\tilde{g}}^{\left(j\right)}$ of the interval and thus “locally supports” the solution function $\tilde{g}$ .

As you can see, the function ψ_j is not defined in a way that it is zero at the boundaries of the interval but at the location of the neighboring intervals’ center values. One may ask why these functions “overlap.” The reason for this is simple: The piecewise-linear profile has two subintervals. The first subinterval in which the function increases and the second subinterval in which the function decreases. Both intervals require a gradient, i.e., we have two unknowns. Therefore it is not sufficient to consider a single interval alone, which only gives us one equation of the form Eq. 29.10. We have already encountered this problem before when introducing the piecewise-linear reconstruction in FVM (see section 31.5.4). We introduced stencil methods as a means of achieving a similar “overlapping” of the reconstruction function. We have to do the same thing here.

32.4.2 Localized Support in Two Dimensions: Pyramid Functions

The basis functions in two dimensions are simply a two-dimensional version of a hat function. These functions are referred to as pyramid functions (see Fig. 32.3b). Depending on the number of triangles that are adjacent in a specific point (in this case point 1) the function involves several neighboring triangles. In all cases, the function ψ_j is defined such that it supports the solution function at the central value ${\tilde{g}}^{\left(j,1\right)}$ and is zero at all other points. Each triangle amounts to one equation of the form of Eq. 29.10. If there are four triangles adjoining in point 1 (as shown in Fig. 32.3b) we have a total of four unknowns to solve for the unknown inclination values of each of the four triangles.

32.4.3 FEM-Galerkin Method

We can now proceed by writing out the N equations for the FEM-Galerkin method by Eq. 29.10 as

${\displaystyle \int {\psi }_j}\left(\overrightarrow{x}\right)L\left(\overrightarrow{x},g,\frac{\mathit{dg}}{dx},\frac{d^{2}g}{d{x}^2}\right)d\Omega =0$

We have modified the notation to account for multiple dependent variables x₁, x₂,…, x_N, and the domain was rewritten to Ω. In this sense, Ω can refer to a one-dimensional domain (linear), a two-dimensional domain (plane), or a three-dimensional domain (volume). We apply this equation to each of the j nodes in our domain. At each node we include all adjacent triangles from which we form the basis function ψ_j. We apply this function to Eq. 32.17, obtaining one equation per node. Please note that in this case the index j does not refer to the individual triangle but to the pyramid made from the triangles that are adjacent in a given point. This is important to keep in mind as we are used to using fixed indices to refer to volumes, intervals, etc. This also implies that the value of the point 1 in Fig. 32.3b when calculating the pyramid j is nonzero whereas it is zero when calculating all other pyramids in the domain Ω. This is due to the nature of localized-support functions.

We will now turn to investigating potential terms that may arise from the linear differential operator L. There is only a limited choice of potential terms that we will need to discuss. We already did a similar discussion of such contributions when deriving FVM (see section 31.3.3). However, in this case, the mathematics is slightly different.

Constants: Sources and Sinks. Source and sink terms yield constants. Obviously, constants are simply integrated over the domain. As already stated, these values are often given as normalized values, e.g., in the form of the volume forces $\overrightarrow{k}$ of the Navier-Stokes equation (see Eq. 11.40), which are simply multiplied with the domain. As stated, the domain may not be a volume. In this case, these values may be given as the “quantity of the dependent variable per area” (two-dimensional) or the “quantity of the dependent variable per length” (one-dimensional). Thus the integration simply is

${\displaystyle \int {\psi }_{j}cd\Omega =c{\displaystyle \int {\psi }_{j}d\Omega }}$

which only amounts to the integration of ψ_j.

Nabla Operators: Divergence. In general, a divergence of the dependent variable g gives rise to the integral

${\displaystyle {\int }_{\Omega }{\psi }_j\left(\overrightarrow{\nabla }\cdot g\right)d\Omega }$

which cannot be simplified further.

Laplace Operators: Diffusion. The third case we will discuss are diffusion effects of the dependent variable’s quantity, which is described by Laplace operators.

$\begin{array}{lll}{\displaystyle {\int }_{\Omega }\omega \Delta gd\Omega }& =& {\displaystyle {\int }_{\Omega }{\psi }_j\left(\frac{{\partial }^{2}g}{\partial x^2}+\frac{{\partial }^{2}g}{\partial y^2}+\frac{{\partial }^{2}g}{\partial z^2}\right)}d\Omega \\ & =& {\displaystyle {\int }_{\Omega }\left(\frac{\partial }{\partial x}\left({\psi }_j\frac{\partial g}{\partial x}\right)-\frac{\partial {\psi }_j}{\partial x}\frac{\partial g}{\partial x}+\frac{\partial }{\partial y}\left({\psi }_j\frac{\partial g}{\partial y}\right)-\frac{\partial {\psi }_j}{\partial y}\frac{\partial g}{\partial y}+\frac{\partial }{\partial z}\left({\psi }_j\frac{\partial g}{\partial z}\right)-\frac{\partial {\psi }_j}{\partial z}\frac{\partial g}{\partial z}\right)d\Omega }\\ & =& {\displaystyle {\int }_{\Omega }\frac{\partial }{\partial x}\left({\psi }_j\frac{\partial g}{\partial x}\right)+\frac{\partial }{\partial y}\left({\psi }_j\frac{\partial g}{\partial y}\right)}+\frac{\partial }{\partial z}\left({\psi }_j\frac{\partial g}{\partial z}\right)d\Omega -\left({\displaystyle {\int }_{\Omega }\frac{\partial {\psi }_j}{\partial x}\frac{\partial g}{\partial x}+\frac{\partial {\psi }_j}{\partial y}\frac{\partial g}{\partial y}+\frac{\partial {\psi }_j}{\partial z}\frac{\partial g}{\partial z}}d\Omega \right)\\ & =& {\displaystyle {\int }_{\Omega }\overrightarrow{\nabla }\cdot \left({\psi }_j\overrightarrow{\nabla }\cdot g\right)d\Omega -{\displaystyle {\int }_{\Omega }\left(\overrightarrow{\nabla }\cdot {\psi }_j\right)\left(\overrightarrow{\nabla }\cdot g\right)d\Omega }}\end{array}$

where we apply Gauss’s theorem (see Eq. 7.10) on the first integral for which we find

${\displaystyle {\int }_{\Omega }\overrightarrow{\nabla }\cdot \left({\psi }_j\overrightarrow{\nabla }\cdot g\right)d\Omega ={\displaystyle {\oint }_{\partial \Omega }\left({\psi }_j\overrightarrow{\nabla }\cdot g\right)\cdot \overrightarrow{n}d\partial \Omega =0}}$

which is zero because our localized-support functions ψ_j are designed such that they are zero on the boundary of the domain. In this case Eq. 32.18 simplifies to

${\displaystyle {\int }_{\Omega }{\psi }_j\mathit{\Delta gd\Omega }=-{\displaystyle {\int }_{\Omega }\left(\overrightarrow{\nabla }\cdot {\psi }_j\right)\left(\overrightarrow{\nabla }\cdot g\right)}}d\Omega$

Different Reconstruction Functions. During our derivation we have used linear interpolation within the individual triangles. Obviously, this is not the only reconstruction scheme that can be used. It is equally possible to use quadratic approximations in the same way we have used these for FVM (see section 31.5.5). These approximations give rise to slightly more complex integrals but often to (significantly) more exact solutions.

32.5 One-Dimensional Example: Flow in Infinitesimally Extended Channels

32.5.1 Problem Definition

We return to the example of the infinitesimally extended channel (see section 16.2) for which we found the following ODE (see vEq. 16.9)

$\frac{d^2{v}_x}{d{z}^2}=\frac{1}{\eta }\frac{dp}{dx}$

We use the same example we have used during the illustration of the Chebyshev polynomials in spectral methods using the values

$\begin{array}{lll}\eta & =& 1\mathrm{m}\mathrm{P}\mathrm{a}\mathrm{s}\\ \frac{dp}{dx}& =& -0.001\mathrm{mbar}{\mathrm{mm}}^{-1}\\ \frac{1}{\eta }\frac{dp}{dx}& =& -100{\mathrm{mm}}^{-1}{\mathrm{s}}^{-1}\end{array}$

which allows us to rewrite Eq. 32.20 to

$\frac{d^2{v}_x}{d{z}^2}+100=0$

where the domain is − 1 ≤ z ≤ 1 (in mm) and the result is obtained in mms−1. The coordinate system is shifted into the center of the channel.

32.5.2 Hat Functions

We will solve this example using a one-dimensional FEM approach with hat functions as basis functions. For this, $\tilde{g}$ will be a function of z only, which is the only independent variable. We will split the domain into intervals of equal size according to Fig. 32.3a. Within the interval j our approximation function $\tilde{g}$ (j) takes the form

${\psi }_j\left(z\right)={\tilde{g}}^{\left(j\right)}\left\{\begin{array}{cc}\hfill \frac{z-z^{\left(j-1\right)}}{z^{\left(j\right)}-z^{\left(j-1\right)}}\hfill & \hfill z^{\left(j-1\right)}\le z\le z^{\left(j\right)}\hfill \\ \hfill \frac{z^{\left(j+1\right)}-z}{z^{\left(j+1\right)}-z^{\left(j\right)}}\hfill & \hfill z^{\left(j\right)}\le z\le z^{\left(j+1\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

Our approximation function $\tilde{g}$ on the domain is thus given by

$\tilde{g}\left(z\right)={\displaystyle \sum _{j=1}^N{\tilde{g}}^{\left(j\right)}{\psi }_j\left(z\right)}$

32.5.3 Integration

Referring to section 4.3 we now set up the integration equation given by Eq. 32.17. The differential operator only has a single entry and contains a constant and a diffusion term. According to the derivation in section 4.3 we can write the integral as

$\begin{array}{lll}{\displaystyle \underset{{\Omega }_j}{\int }{\psi }_j\left(\frac{d^2\tilde{g}}{d{z}^2}+100\right)}d\Omega & =& 0\\ -{\displaystyle \underset{{\Omega }_j}{\int }\left(\overrightarrow{\nabla }\cdot {\psi }_j\right)}\left(\overrightarrow{\nabla }\cdot \tilde{g}\right)d\Omega +100{\displaystyle \underset{{\Omega }_j}{\int }{\psi }_{j}d\Omega }& =& 0\\ {\displaystyle \underset{{\Omega }_j}{\int }\frac{\partial {\psi }_j}{\partial z}\frac{\partial \tilde{g}}{dz}d\Omega }& =& 100{\displaystyle \underset{{\Omega }_j}{\int }{\psi }_{j}d\Omega }\end{array}$

which amounts to one equation per interval j. We now require the first derivative of our approximation function $\tilde{g}$ , which is given by

$\frac{\partial \tilde{g}}{\partial z}={\displaystyle \sum _{i=1}^N{\tilde{g}}^{\left(i\right)}\frac{\partial {\psi }_i}{\partial z}}$

Please note that we have to use the index i because the derivative within the interval j will also include contributions from ${\tilde{g}}^{\left(j-1\right)}$ and ${\tilde{g}}^{\left(j+1\right)}$ , which is why we need to use a separate index variable. Using the definition of ${\tilde{g}}^{\left(j\right)}$ we can now derive the partial differential and the integral given by

$\frac{\partial {\psi }_j}{\partial z}=\left\{\begin{array}{cc}\hfill \frac{1}{z^{\left(j\right)}-z^{\left(j-1\right)}}\hfill & \hfill z^{\left(j-1\right)}\le z\le z^{\left(j\right)}\hfill \\ \hfill -\frac{1}{z^{\left(j+1\right)}-z^{\left(j\right)}}\hfill & \hfill z^{\left(j\right)}\le z\le z^{\left(j+1\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\begin{array}{l}{\displaystyle \underset{{\Omega }_j}{\int }{\psi }_{j}d\Omega }={\displaystyle \underset{z^{\left(j-1\right)}}{\overset{z^{\left(j\right)}}{\int }}\frac{z-z^{\left(j-1\right)}}{z^{\left(j\right)}-z^{\left(j-1\right)}}dz+}{\displaystyle \underset{z^{\left(j-1\right)}}{\overset{z^{\left(j\right)}}{\int }}\frac{z-z^{\left(j-1\right)}}{z^{\left(j\right)}-z^{\left(j-1\right)}}dz}\\ ={\displaystyle {\left[\frac{\frac{z^2}{2}-z{z}^{\left(j-1\right)}}{z^{\left(j\right)}-z^{\left(j-1\right)}}\right]}_{z^{\left(j-1\right)}}^{z^{\left(j\right)}}}+{\displaystyle {\left[\frac{z{z}^{\left(j-1\right)}-\frac{z^2}{2}}{z^{\left(j-1\right)}-z^{\left(j\right)}}\right]}_{z^{\left(j\right)}}^{z^{\left(j+1\right)}}}\\ =\frac{\frac{{\left(z^{\left(j\right)}\right)}^2-{\left(z^{\left(j-1\right)}\right)}^2}{2}-\left(z^{\left(j\right)}-z^{\left(j-1\right)}\right)z^{\left(j-1\right)}}{z^{\left(j\right)}-z^{\left(j-1\right)}}+\frac{\left(z^{\left(j+1\right)}-z^{\left(j\right)}\right)z^{\left(j+1\right)}-\frac{{\left(z^{\left(j+1\right)}\right)}^2-{\left(z^{\left(j\right)}\right)}^2}{2}}{z^{\left(j+1\right)}-z^{\left(j\right)}}\end{array}$

As you can see, these values depend solely on the geometry of the chosen mesh and can thus be precalculated. In this example, we assume a total of 10 intervals of equal size. Our boundary − 1 ≤ z ≤ 1 is thus subdivided into intervals of size Δz = 0.2. Using this simplification we can simplify Eq. 32.22, Eq. 32.26, and Eq. 32.27 to

${\tilde{g}}^{\left(j\right)}={\tilde{g}}^{\left(j\right)}\left\{\begin{array}{cc}\hfill \frac{z-z^{\left(j-1\right)}}{\Delta z},\hfill & \hfill z^{\left(j-1\right)}\le z\le z^{\left(j\right)}\hfill \\ \hfill \frac{z^{\left(j+1\right)}-z}{\Delta z},\hfill & \hfill z^{\left(j\right)}\le z\le z^{\left(j+1\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\frac{\partial {\psi }_j}{\partial z}=\left\{\begin{array}{cc}\hfill \frac{1}{\Delta z},\hfill & \hfill z^{\left(j-1\right)}\le z\le z^{\left(j\right)}\hfill \\ \hfill -\frac{1}{\Delta z},\hfill & \hfill z^{\left(j\right)}\le z\le z^{\left(j+1\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\begin{array}{l}{\displaystyle \underset{{\Omega }_j}{\int }\begin{array}{l}{\psi }_{j}d\Omega =\frac{\frac{{\left(z^{\left(j\right)}\right)}^2-{\left(z^{\left(j-1\right)}\right)}^2}{2}-\Delta z{z}^{\left(j-1\right)}+\Delta z{z}^{\left(j+1\right)}-\frac{{\left(z^{\left(j+1\right)}\right)}^2-{\left(z^{\left(j\right)}\right)}^2}{2}}{\Delta z}\\ =\frac{\frac{{\left(z^{\left(j\right)}\right)}^2-{\left(z^{\left(j\right)}-\Delta x\right)}^2}{2}-\Delta z\left(z^{\left(j\right)}-\Delta z\right)+\Delta z\left(z^{\left(j\right)}+\Delta z\right)-\frac{{\left(z^{\left(j\right)}+\Delta z\right)}^2-{\left(z^{\left(j\right)}\right)}^2}{2}}{\Delta z}\end{array}}\\ =\frac{2\frac{z^{\left(j\right)}\Delta z-{\left(\Delta x\right)}^2}{2}-\Delta z\left(z^{\left(j\right)}-\Delta z\right)+\Delta z\left(z^{\left(j\right)}+\Delta z\right)-\frac{2z^{\left(j\right)}\Delta z+{\left(\Delta z\right)}^2}{2}}{\Delta z}\\ =\Delta z\end{array}$

Obviously, these are some very neat simplifications. In most cases, we cannot make these types of simplifications because the mesh is not regular and thus the width of the individual interval is different.

32.5.4 Matrix Notation

We now combine Eq. 32.24 and Eq. 32.30 finding

${\displaystyle \underset{{\Omega }_j}{\int }\frac{\partial {\psi }_j}{\partial z}{\displaystyle \sum _{i=1}^N{\tilde{g}}^{\left(i\right)}\frac{\partial {\psi }_i}{\partial z}}}d\Omega =100{\displaystyle \underset{{\Omega }_j}{\int }{\psi }_{j}d\Omega }$

We can now insert Eq. 32.30 in the right-hand side of Eq. 32.31 finding

${\displaystyle \underset{{\Omega }_j}{\int }\frac{\partial {\psi }_j}{\partial z}{\displaystyle \sum _{i=1}^N{\tilde{g}}^{\left(i\right)}\frac{\partial {\psi }_i}{\partial z}}}d\Omega =100\Delta z$

We need a little more caution when evaluating the left-hand side of Eq. 32.32. At interval j we only need to take into account the terms ${\displaystyle \sum _{i=1}^N{\tilde{g}}^{\left(i\right)}\frac{\partial {\psi }_i}{\partial z}}$ with i = j − 1, i = j, and i = j + 1. Remember that we designed the local support functions such that they overlap by one interval. We therefore can rewrite Eq. 32.32 to

${\displaystyle \underset{{\Omega }_j}{\int }\frac{\partial {\psi }_j}{\partial z}}\left({\tilde{g}}^{\left(j-1\right)}\frac{\partial {\psi }_{j-1}}{\partial z}+{\tilde{g}}^{\left(j\right)}\frac{\partial {\psi }_j}{\partial z}+{\tilde{g}}^{\left(j+1\right)}\frac{\partial {\psi }_{j+1}}{\partial z}\right)d\Omega =100\Delta z$

We note that the first integral of Eq. 32.33 must only be evaluated on the range z^(j) – Δz ≤ z ≤ z^(j) because ${\psi }_{j-1}$ is zero for z > z^(j). We therefore find

$\begin{array}{l}{\tilde{g}}^{\left(j-1\right)}{\displaystyle \underset{z^{\left(j\right)}-\Delta z}{\overset{z^{\left(j\right)}}{\int }}\frac{\partial {\psi }_j}{\partial z}}\frac{\partial {\psi }_{j-1}}{\partial z}dz={\tilde{g}}^{\left(j-1\right)}{\displaystyle \underset{z^{\left(j\right)}-\Delta z}{\overset{z^{\left(j\right)}}{\int }}\frac{1}{{\left(\Delta z\right)}^2}}\\ =-{\tilde{g}}^{\left(j-1\right)}{\displaystyle {\left[\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(j\right)}-\Delta z}^{z^{\left(j\right)}}}\\ =-\frac{{\tilde{g}}^{\left(j-1\right)}}{\Delta z}\end{array}$

Likewise, the third integral of Eq. 32.33 must only be evaluated on the range z^(j) ≤ z ≤ z^(j) + Δz because ψ_j is zero for z < z(j). We therefore find

$\begin{array}{l}{\tilde{g}}^{\left(j+1\right)}{\displaystyle \underset{z^{\left(j\right)}}{\overset{z^{\left(j\right)}+\Delta z}{\int }}\frac{\partial {\psi }_j}{\partial z}}\frac{\partial {\psi }_{j+1}}{\partial z}dz={\tilde{g}}^{\left(j+1\right)}{\displaystyle \underset{z^{\left(j\right)}}{\overset{z^{\left(j\right)}+\Delta z}{\int }}\frac{1}{{\left(\Delta z\right)}^2}}\\ =-{\tilde{g}}^{\left(j+1\right)}{\displaystyle {\left[\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(j\right)}}^{z^{\left(j\right)}+\Delta z}}\\ =-\frac{{\tilde{g}}^{\left(j+1\right)}}{\Delta z}\end{array}$

The middle integral in Eq. 32.33 must be evaluated from z^(j) − Δz ≤ z ≤ z^(j) + Δz because ψ_j is nonzero in this range. We therefore find

$\begin{array}{l}{\tilde{g}}^{\left(j\right)}{\displaystyle \underset{z^{\left(j\right)}-\Delta z}{\overset{z^{\left(j\right)}+\Delta z}{\int }}\frac{\partial {\psi }_j}{\partial z}}\frac{\partial {\psi }_j}{\partial z}dz={\tilde{g}}^{\left(j\right)}{\displaystyle \underset{z^{\left(j\right)}-\Delta z}{\overset{z^{\left(j\right)}+\Delta z}{\int }}\frac{1}{{\left(\Delta z\right)}^2}}\\ =-{\tilde{g}}^{\left(j\right)}{\displaystyle {\left[\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(j\right)}-\Delta z}^{z^{\left(j\right)}+\Delta z}}\\ =2\frac{{\tilde{g}}^{\left(j\right)}}{\Delta z}\end{array}$

Inserting Eq. 32.34, Eq. 32.35, and Eq. 32.36 into Eq. 32.33 yields

$\begin{array}{l}\frac{1}{\Delta z}\left(-{\tilde{g}}^{\left(j-1\right)}+2{\tilde{g}}^{\left(j\right)}-{\tilde{g}}^{\left(j+1\right)}\right)=100\Delta z\\ -{\tilde{g}}^{\left(j-1\right)}+2{\tilde{g}}^{\left(j\right)}-{\tilde{g}}^{\left(j+1\right)}=100{\left(\Delta z\right)}^2\end{array}$

where we have a linear equation for the interval j that involves the unknown values ${\tilde{g}}^{\left(j-1\right)}$ , ${\tilde{g}}^{\left(j\right)}$ , ${\tilde{g}}^{\left(j+1\right)}$ .

Finite-Difference Scheme. If you look at Eq. 32.37 you can see that this is a finite difference approximation for the second derivative according to Eq. 4.10. This is due to the fact that we used a regular mesh for this FEM approach. We would also obtain this equation from an FDM approach. As an example, refer to Eq. 30.2 where we have obtained this equation in two-dimensional form. We could use the spreadsheet prepared for Fig. 30.5 directly to solve this problem using an FDM approach.

Left-Hand Boundary Values. If you refer to Fig. 32.3a you can see that we have to evaluate a slightly different integral for the first and the Last intervals in the domain in order to account for the boundary values. This is due to the fact that the basis functions ψ₀ and ψ_N are only defined on half of the interval. The value $\stackrel{\sim \left(0\right)}{g}$ is the boundary value of the left-hand side of the domain. For this we find Microfluidics: Modeling, Mechanics, and Mathematics

${\tilde{g}}^{\left(0\right)}={\tilde{g}}^{\left(0\right)}\left\{\begin{array}{cc}\hfill \frac{z^{\left(1\right)}-z}{\Delta z}\hfill & \hfill z^{\left(0\right)}\le z\le z^{\left(1\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\frac{\partial {\psi }_0}{\partial z}=\left\{\begin{array}{cc}\hfill -\frac{1}{\Delta z}\hfill & \hfill z^{\left(0\right)}\le z\le z^{\left(1\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_0}{\psi }_0}={\left[\frac{z{z}^{\left(1\right)}-\frac{z^2}{2}}{\Delta z}\right]}_{z\left(0\right)}^{z\left(1\right)}=\frac{\Delta z{z}^{\left(1\right)}-\frac{\left({\left(z^{\left(1\right)}\right)}^2-{\left(z^{\left(1\right)}-\Delta z\right)}^2\right)}{2}}{\Delta z}\\ =\frac{\Delta z}{2}\end{array}$

which we require to calculate Eq. 32.33 for the first interval j = 0. Since the domain Ω0 is ${\tilde{g}}^{\left(0\right)}\le z\le {\tilde{g}}^{\left(1\right)}$ Eq. 32.33 simplifies to

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_0}\frac{\partial {\psi }_0}{\partial z}}\left({\tilde{g}}^{\left(0\right)}\frac{\partial {\psi }_0}{\partial z}+{\tilde{g}}^{\left(1\right)}\frac{\partial {\psi }_1}{\partial z}\right)d\Omega =100{\displaystyle {\int }_{{\Omega }_0}{\psi }_{0}d\Omega }\\ {\tilde{g}}^{\left(0\right)}{\displaystyle {\int }_{{\Omega }_0}\frac{\partial {\psi }_0}{\partial z}}\frac{\partial {\psi }_0}{\partial z}d\Omega +{\tilde{g}}^{\left(1\right)}{\displaystyle {\int }_{{\Omega }_0}\frac{\partial {\psi }_0}{\partial z}}\frac{\partial {\psi }_1}{\partial z}d\Omega =100\frac{\Delta z}{2}\end{array}$

$\begin{array}{l}{\tilde{g}}^{\left(0\right)}{\left[\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(0\right)}}^{z^{\left(1\right)}}+{\tilde{g}}^{\left(1\right)}{\left[-\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(0\right)}}^{z^{\left(1\right)}}=50\Delta z\\ \frac{{\tilde{g}}^{\left(0\right)}}{\Delta z}-\frac{{\tilde{g}}^{\left(1\right)}}{\Delta z}=50\Delta z\\ {\tilde{g}}^{\left(0\right)}-{\tilde{g}}^{\left(1\right)}=50{\left(\Delta z\right)}^2\end{array}$

where we have the equation of the left-hand side boundary of the domain.

Right-Hand Boundary. We repeat this calculation for the right-hand side boundary of the domain. The last interval N in our domain is only defined for z^(N-1) ≤ z ≤ z^(N), where ${\tilde{g}}^{\left(N\right)}$ is the boundary value. As stated, ψ_N is only defined on a half-interval. For this we find

${\tilde{g}}^{\left(N\right)}={\tilde{g}}^{\left(N\right)}\left\{\begin{array}{cc}\hfill \frac{z-z^{\left(N-1\right)}}{\Delta z}\hfill & \hfill z^{\left(N-1\right)}\le z\le z^{\left(N\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\frac{\partial {\psi }_N}{\partial z}=\left\{\begin{array}{cc}\hfill -\frac{1}{\Delta z}\hfill & \hfill z^{\left(N-1\right)}\le z\le z^{\left(N\right)}\hfill \\ \hfill 0\hfill & \hfill \mathrm{otherwise}\hfill \end{array}\right.$

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_N}{\psi }_N}={\left[\frac{\frac{z^2}{2}}{\Delta z}-z{z}^{\left(N-1\right)}\right]}_{z\left(N-1\right)}^{z\left(N\right)}=\frac{\frac{{\left(z^{\left(N-1\right)}+\Delta z\right)}^2-{\left(z^{\left(N-1\right)}\right)}^2}{2}-\Delta z{z}^{\left(N-1\right)}}{\Delta z}\\ =\frac{\Delta z}{2}\end{array}$

which we use to calculate the integral from Eq. 32.33 of j = N. Since the domain ΩN is ${\tilde{g}}^{\left(N-1\right)}\le z\le {\tilde{g}}^{\left(N\right)}$ Eq. 32.33 simplifies to

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_N}\frac{\partial {\psi }_N}{\partial z}}\left({\tilde{g}}^{\left(N-1\right)}\frac{\partial {\psi }_{N-1}}{\partial z}+{\tilde{g}}^{\left(N\right)}\frac{\partial {\psi }_N}{\partial z}\right)d\Omega =100{\displaystyle {\int }_{{\Omega }_N}{\psi }_{0}d\Omega }\\ {\tilde{g}}^{\left(N-1\right)}{\displaystyle {\int }_{{\Omega }_N}\frac{\partial {\psi }_N}{\partial z}}\frac{\partial {\psi }_{N-1}}{\partial z}d\Omega +{\tilde{g}}^{\left(N\right)}{\displaystyle {\int }_{{\Omega }_N}\frac{\partial {\psi }_N}{\partial z}}\frac{\partial {\psi }_N}{\partial z}d\Omega =100\frac{\Delta z}{2}\end{array}$

$\begin{array}{l}{\tilde{g}}^{\left(N-1\right)}{\left[-\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(N-1\right)}}^{z^{\left(N\right)}}+{\tilde{g}}^{\left(N\right)}{\left[-\frac{z}{{\left(\Delta z\right)}^2}\right]}_{z^{\left(N-1\right)}}^{z^{\left(N\right)}}=50\Delta z\\ -\frac{{\tilde{g}}^{\left(N-1\right)}}{\Delta z}+\frac{{\tilde{g}}^{\left(N\right)}}{\Delta z}=50\Delta z\\ -{\tilde{g}}^{\left(N-1\right)}+{\tilde{g}}^{\left(N\right)}=50{\left(\Delta z\right)}^2\end{array}$

Where we have the equation of the right-hand side boundary of the domain.

Stiffness Matrix and Load Vector. With Eq. 32.42 for the left-hand boundary of the domain, Eq. 32.47 for the right-hand boundary, and Eq. 32.37 for all intervals between we can now write out the individual equations for all intervals. Inserting the value Δz = 0.2, i.e., accounting for a total of 10 intervals on the domain, we find

$\begin{array}{ccccc}\hfill {\tilde{g}}^{\left(0\right)}-{\tilde{g}}^{\left(1\right)}\hfill & \hfill =\hfill & \hfill 2\hfill & \hfill j=0\hfill & \hfill \left(\mathrm{left}‐\mathrm{hand}\mathrm{boundary}\right)\hfill \\ \hfill -{\tilde{g}}^{\left(j-1\right)}+2{\tilde{g}}^{\left(j\right)}-{\tilde{g}}^{\left(j+1\right)}\hfill & \hfill =\hfill & \hfill 4\hfill & \hfill 1\le j\le N-1\hfill & \hfill \left(\mathrm{all}\mathrm{nonboundary}\mathrm{intervals}\right)\hfill \\ \hfill -{\tilde{g}}^{\left(N-1\right)}+{\tilde{g}}^{\left(N\right)}\hfill & \hfill =\hfill & \hfill 2\hfill & \hfill j=N\hfill & \hfill \left(\mathrm{right}‐\mathrm{hand}\mathrm{boundary}\right)\hfill \end{array}$

If we arrange the 10 intervals we obtain a total of 12 linear equations. The two additional equations are due to the boundary conditions. Written out in a matrix we obtain the following form:

$\left(\begin{array}{cccccccccccc}\hfill -1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\tilde{g}}^{\left(0\right)}\hfill \\ \hfill {\tilde{g}}^{\left(1\right)}\hfill \\ \hfill {\tilde{g}}^{\left(2\right)}\hfill \\ \hfill {\tilde{g}}^{\left(3\right)}\hfill \\ \hfill {\tilde{g}}^{\left(4\right)}\hfill \\ \hfill {\tilde{g}}^{\left(5\right)}\hfill \\ \hfill {\tilde{g}}^{\left(6\right)}\hfill \\ \hfill {\tilde{g}}^{\left(7\right)}\hfill \\ \hfill {\tilde{g}}^{\left(8\right)}\hfill \\ \hfill {\tilde{g}}^{\left(9\right)}\hfill \\ \hfill {\tilde{g}}^{\left(10\right)}\hfill \\ \hfill {\tilde{g}}^{\left(11\right)}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill 2\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 2\hfill \end{array}\right)$

Eq. 32.48 in a regular linear system of type A$\overrightarrow{x}=\overrightarrow{b}$ as we have studied them in section 25.3. In FEM the matrix A is referred to as the stiffness matrix. The vector $\overrightarrow{b}$ is referred to as the load vector.

The matrix of Eq. 32.48 contains a total of 10 intervals, but we have written a total of 12 equations in order to account for the boundary conditions, where ${\tilde{g}}^{\left(0\right)}$ is the value on the left-hand side boundary and ${\tilde{g}}^{\left(11\right)}$ is the value at the right-hand side boundary. If the boundary condition is Neumann type, we do not know these values, only their relation to the neighboring values. As an example, if the boundary values are expected to have a zero-gradient we could change the system such that ${\tilde{g}}^{\left(0\right)}={\tilde{g}}^{\left(1\right)}$ and ${\tilde{g}}^{\left(11\right)}={\tilde{g}}^{\left(10\right)}$ .

However, in our example, the boundary values are given. As the boundary condition is the no-slip condition we know that ${\tilde{g}}^{\left(0\right)}={\tilde{g}}^{\left(11\right)}=0$ . Therefore, our system of equation reduces to 10 equations, as we can drop the first and the last equations brought in by the boundary conditions. This simplifies the system to

$\left(\begin{array}{cccccccccc}\hfill -2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 2\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\tilde{g}}^{\left(1\right)}\hfill \\ \hfill {\tilde{g}}^{\left(2\right)}\hfill \\ \hfill {\tilde{g}}^{\left(3\right)}\hfill \\ \hfill {\tilde{g}}^{\left(4\right)}\hfill \\ \hfill {\tilde{g}}^{\left(5\right)}\hfill \\ \hfill {\tilde{g}}^{\left(6\right)}\hfill \\ \hfill {\tilde{g}}^{\left(7\right)}\hfill \\ \hfill {\tilde{g}}^{\left(8\right)}\hfill \\ \hfill {\tilde{g}}^{\left(9\right)}\hfill \\ \hfill {\tilde{g}}^{\left(10\right)}\hfill \end{array}\right)=\left(\begin{array}{c}\hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \\ \hfill 4\hfill \end{array}\right)$

The only thing that is left to do is to solve Eq. 32.49 using any suitable method for solving linear systems of equations.

32.5.5 Maple Worksheet

For solving Eq. 32.49 we will make use of the Thomas algorithm (see section 25.3.8) and our implemented function DoThomas from listing 25.12. The listing is shown in listing 32.1. As you can see, it is fairly short. It repeats the example we have already used in listing 25.12, where we actually anticipated this example. After restarting the

 1 restart:read "Core.txt":read "CoreFunctions.txt":
 2
 3 A:=Matrix ([[2 ,−1,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0],
 4             [− 1,2 , − 1,0 ,0 ,0 ,0 ,0 ,0 ,0],
 5             [ 0, − 1,2 , − 1,0 ,0 ,0 ,0 ,0 ,0],
 6             [ 0, 0, − 1, 2, − 1,0 ,0 ,0 ,0 ,0],
 7             [ 0, 0, 0, − 1,2 , − 1,0 ,0 ,0 ,0],
 8             [ 0, 0, 0, 0, − 1,2 , − 1,0 ,0 ,0],
 9             [ 0, 0, 0, 0,0 , − 1,2 , − 1,0 ,0],
10             [ 0, 0, 0, 0,0 ,0 , − 1,2 , − 1,0],
11             [ 0, 0, 0, 0,0 ,0 ,0 , − 1,2 , − 1],
12             [ 0, 0, 0, 0,0 ,0 ,0 , 0, − 1,2]]);
13 b:=Vector( [4,4,4,4,4,4,4,4,4,4] );
14 vx:=DoThomas ( A, b);

Maple server and including the files Core.txt and CoreFunctions.txt we define the matrix A in line 3 and the vector $\overrightarrow{b}$ in line 13. We then call the functions DoThomas, obtaining the solution

${\overrightarrow{v}}_x=\left(\begin{array}{c}\hfill 20\hfill \\ \hfill 36\hfill \\ \hfill 48\hfill \\ \hfill 56\hfill \\ \hfill 60\hfill \\ \hfill 60\hfill \\ \hfill 56\hfill \\ \hfill 48\hfill \\ \hfill 36\hfill \\ \hfill 20\hfill \end{array}\right)$

We already know from Eq. 32.28 that the correct solution is

v_x (z) = 50 (1 − z²)

If we plot the results from line 14 against the analytical solution we can see that the results are slightly off (see Fig. 32.4). This is due to the fact that we chose a very rough mesh with a spacing of 200 μm, which is too coarse for obtaining good solutions. However, if the mesh is refined, i.e., Δz is reduced, the solution becomes more and more exact. For Δz = 100 μm the solution is already significantly closer to the analytical solution. This requires a total of only 20 intervals.

Plotting the Hat Functions. In order to illustrate the hat functions, i.e., our basis functions ψj we plot the functions using the resulting values ${\tilde{g}}^{\left(j\right)}$ obtained from listing 32.1. The result is shown in Fig. 32.5. As you can see, the hat functions locally support the solution $\tilde{g}$ (x) which is linearly interpolated between the support points.

32.5.6 General Notation

As we have seen, it is not complicated to apply FEM to solving PDEs in one dimension. Ultimately, we are left with solving a linear system of equations of the form given in Eq. 32.49. We have derived the equations for the stiffness matrix A and the load vector $\overrightarrow{b}$ analytically in our example. We will now give the general notation for them depending on the type of boundary conditions.

Known Boundary Values. In this case we do not require the calculation of the boundary values as they are given. In this case, the entries a_j,k of the matrix A can be calculated according to Eq. 32.33 as

${\displaystyle a_{j,k}}=\left\{\begin{array}{cc}\hfill {\displaystyle {\int }_{{\Omega }_j}\frac{\partial {\psi }_j}{\partial z}{\tilde{g}}^{\left(j-1\right)}\frac{\partial {{\psi }_j}_{-1}}{\partial z}d\Omega }\hfill & \hfill k=j-1\ge 1\hfill \\ \hfill {\displaystyle {\int }_{{\Omega }_j}\frac{\partial {\psi }_j}{\partial z}{\tilde{g}}^{\left(j\right)}\frac{\partial {\psi }_j}{\partial z}d\Omega }\hfill & \hfill k=j\hfill \\ \hfill {\displaystyle {\int }_{{\Omega }_j}\frac{\partial {\psi }_j}{\partial z}{\tilde{g}}^{\left(j\mathit{+}\mathit{1}\right)}\frac{\partial {{\psi }_j}_{\mathit{+}\mathit{1}}}{\partial z}d\Omega }\hfill & \hfill k=j+1\le N\hfill \\ \hfill 0\hfill & \hfill \text{otherwise}\hfill \end{array}\right.$

where j is the row and k is the column index. Likewise the entries of the load vector $\overrightarrow{b}$ can be calculated according to Eq. 32.31 as

${\displaystyle b_j}=\begin{array}{cc}\hfill {\displaystyle {\int }_{{\Omega }_j}\frac{\partial {\psi }_j}{\partial z}\left({\tilde{g}}^{\left(j-1\right)}\frac{\partial {\psi }_{j-1}}{\partial z}+{\tilde{g}}^{\left(j\right)}\frac{\partial {\psi }_j}{\partial z}+{\tilde{g}}^{\left(j+1\right)}\frac{\partial {\psi }_{j+1}}{\partial z}\right)d\text{Ω},}\hfill & \hfill 1\le j\le N\hfill \end{array}$

where the value ${\tilde{g}}^{\left(0\right)}$ is the left-hand side boundary value and the value ${\tilde{g}}^{\left(N\right)}$ is the right-hand side boundary value.

Unknown Boundary Values. In case the boundary values are not known they must be calculated as well. In this case, the first (j = 0) and the last (j = N + 1) rows of the stiffness matrix are calculated differently. We find the following from Eq. 32.41 and Eq. 32.46:

${\displaystyle a_{0,k}}=\left\{\begin{array}{cc}\hfill {\displaystyle {\int }_{{\Omega }_j}\frac{\partial {\psi }_j}{\partial z}{\tilde{g}}^{\left(j\right)}\frac{\partial {\psi }_j}{\partial z}d\Omega }\hfill & \hfill k=j\hfill \\ \hfill {\displaystyle {\int }_{{\Omega }_j}\frac{\partial {\psi }_j}{\partial z}{\tilde{g}}^{\left(j+1\right)}\frac{\partial {{\psi }_{j+}}_1}{\partial z}d\Omega }\hfill & \hfill k=j+1\le N\hfill \\ \hfill 0\hfill & \hfill \text{otherwise}\hfill \end{array}\right.$