The Finite Element Method (FEM) is a mathematical tool using for solving differential equations. It has been popularly used (also, initially developed) for solving elastostatics problems which is to find the displacement $u$ such that,
$$
\sigma_{ij,j} + f_i = 0 \text{ in the domain }\Omega
$$
To solve this, the domain $\Omega$ is split into finite sub domains $\Omega_e$ in which the solution $u$ is (usually) approximated by a polynomial of certain order which is governed by the regularity requirements of the weak form of above PDE.
The division into finite elements $\Omega_e$ is usually done as
Lines for 1D problem (obviously)
Quadrilaterals in 2D problems
Hexahedrons in 3D problems
These finite elements are then mapped to a “Natural co-ordinate” system where they accordingly become
Lines with $\xi_1 \in (-1,1)$
Square with $\xi_1, \xi_2 \in (-1,1)$
Cubes with $\xi_1, \xi_2, \xi_3 \in (-1,1)$
where $\xi_1, \xi_2, \xi_3$ are the $x,y,z$ co-ordinates in this bi-unit domain (called so since these values here range from $-1$ to $1$).
While the FEM is solved for discrete values of $u$, the solution is presented as a continuous function over the domain as well as within the elements. This equation of the solution is obtained by scaling the Lagrangian polynomials in this space by the nodal values of $u$ (it is these nodal values that FEM finally solves for). These Lagrangian polynomials are (usually) the basis functions of the finite elements.
1-D Elements
See
The above diagram shows the nodes for 1D elements. There are two main points to be noted:
There are polynomials for each node of an element. And these polynomials are chosen such a way that these polynomials have value $1$ at that node and $0$ at all other nodes.
Order of the polynomial $=\text{Number of nodes}-1$
The final solution is the sum of all these polynomials which can be given as
$$
u(\xi) = \sum_{A=1}^{N_{\text{nodes}}}N_A(\xi)\bf u_A
$$
where,
$N_A$ is the basis function at node A.
$\bf u_A$ is nodal degree of freedom (value of displacement).
$N_\text{nodes}$ is the total number of nodes per element.
The general Lagrangian basis function $N_A$ is given as,
$$
N_A(\xi) = \prod_{B=1 \\ B\not=A}^{N_\text{nodes}}\frac {\xi - \xi_B}{\xi_A - \xi_B}
$$
where $B$ iterates through all the nodes except A. It can be observed here that the value of basis function is always zero at all the nodes (except A).
1
2
3
4
5
6
7
8
9
10
11
12
# The generic Lagrangian basis function isfunctionN(ξ,order,node)ξ_B=range(-1,1,length=order+1)|>collectξ_A=ξ_B[node]N_A=1foriinrange(1,order+1)ifi≠nodeN_A*=(ξ-ξ_B[i])/(ξ_A-ξ_B[i])endendN_Aend
Note that, we won’t be using the generic function for the ease of comparison of various orders. Now that we are moving into plotting of basis functions, let me initialize the prerequisites in Julia.
The superscript 1 denotes that it is a linear element. The same basis functions can be obtained by calling the above function as N(ξ, 1, node) for each of the respective node. The corresponding plots are,
For the case of 2D elements, we will consider the quadratic elements which can be obtained as a tensor product of the above 1D basis functions. This can be represented as,
Generally, an element has the same order of the basis functions along both $\xi_1$ and $\xi_2$ directions though it need not be the case.
Thus, for a quadilateral element, the basis function in the natural co-ordinates (where it is transformed to square elements) over each node can be obtained by tensor product as,
Here is the plot for the $N_{11}$. The $z$ axis represents the nodal values of the shape function while the $xy$ plane represents the vector field of $\xi_1$ and $ \xi_2$.
One interesting point to observe here is that, though the function is linear, the plot is curved along the diagonal. That is because the shape function is actually a bilinear function of $\xi_1$ and $\xi_2$. Here is the plot for each shape functions.
2D Quadratic Elements
The quadratic elements are obtained through the tensor product as,
Thus a general representation of the Lagrangian basis functions have been achieved for the 1D and 2D elements. I do wish to implement the same for the 3D elements but I seem to be running out of spaces! (Actually, it is possible and will implement and update with link soon). This had been a very good learning experience for me and am looking forward to dwelving more into this realm soon!
