FEATool has been designed with flexibility and user customization in mind allowing for solving custom equations, calling external solvers, grid generators, and post processing tools, and allowing for both MATLAB m-script as well as GUI based modeling. This post shows how to easily define custom finite element shape functions (also called basis functions or ansatz functions) which can be of interest for researchers, teachers, and advanced FEM users.

In the following a fifth order Lagrangian 1D FEM shape function will
be defined, allowing for very accurate polynomial solutions with few
grid cells. The first step is to create the function file, in this
case we will call it
*sf_line_P5.m*. The
first header line of a shape function file should look like the
following

```
function [ vBase, nLDof, xLDof, sfun ] = sf_line_P5( i_eval, n_sdim, n_vert, i_dof, xi, aInvJac, vBase )
```

The input arguments are; **i_eval** the evaluation type flag,
**n_sdim** number of space dimensions, **n_vert** the number of
vertices per grid cell, **i_dof** the shape or basis function to
evaluate, **xi** local cell evaluation coordinates, **aInvJac**
inverse of the transformation Jacobian, and **vBase** the output
vector.

The output **nLDof** specifies the number of shape functions (degrees
of freedom/unknowns) associated with grid cell vertices, edges (in 2D
and 3D), faces (in 3D), and cell interiors. In this case we have two
shape functions for each of the line end points and four equally
spaced internal ones, thus **nLDof** here is

```
nLDof = [ 2 0 0 4 ];
```

**xLDof** specifies the location or position of each shape
function/dof in local cell coordinates. Barycentric (area/volume)
coordinates are used for lines and simplex cell shapes (triangles and
tetrahedra), and *-1 … 1* local reference coordinates for
quadrilateral and hexahedral cells. It is important to number and
order the shape functions and corresponding coordinates associated
with vertices first, followed by edges, faces, and cell lastly
interiors (in counterclockwise direction in 2D and 3D). For the one
dimensional $P_5$ element we thus have the following coordinates
`[ 1 0 4/5 3/5 2/5 1/5 ]`

and due to using barycentric
coordinates the second coordinate row is simply *1* minus the first
one ($\xi_2 = 1 - \xi_1$), we therefore have

```
xLDof = [ 1 0 4/5 3/5 2/5 1/5 ;
0 1 1/5 2/5 3/5 4/5 ];
```

**sfun** simply returns the function name (for quick and easy function
access during nested calls)

```
sfun = 'sf_line_P5';
```

The **i_eval** input flag is used to determine if either function
value evaluations should be computed (*i_eval = 1*), or x, y, or z
derivatives (*i_eval = 2, 3, or 4*, respectively). Furthermore,
**i_dof** specifies for which of the basis functions evaluation should
be performed, and **vBase** is the final output vector. We can now
program the rest of the shape function evaluation definitions for the
$P_5$ line element. (MATLAB source code for other shape function
definitions can be found in the featool/ellib folder.)

```
switch i_eval
case 1 % Function values.
switch i_dof % Shape/basis function.
case 1
vBase = (625*xi(1)^5)/24 - (625*xi(1)^4)/12 + (875*xi(1)^3)/24 - (125*xi(1)^2)/12 + xi(1);
case 2
vBase = - (625*xi(1)^5)/24 + (625*xi(1)^4)/8 - (2125*xi(1)^3)/24 + (375*xi(1)^2)/8 - (137*xi(1))/12 + 1;
case 3
vBase = - (3125*xi(1)^5)/24 + (6875*xi(1)^4)/24 - (5125*xi(1)^3)/24 + (1525*xi(1)^2)/24 - (25*xi(1))/4;
case 4
vBase = (3125*xi(1)^5)/12 - 625*xi(1)^4 + (6125*xi(1)^3)/12 - (325*xi(1)^2)/2 + (50*xi(1))/3;
case 5
vBase = - (3125*xi(1)^5)/12 + (8125*xi(1)^4)/12 - (7375*xi(1)^3)/12 + (2675*xi(1)^2)/12 - 25*xi(1);
case 6
vBase = (3125*xi(1)^5)/24 - (4375*xi(1)^4)/12 + (8875*xi(1)^3)/24 - (1925*xi(1)^2)/12 + 25*xi(1);
end
case 2 % x-derivative values.
switch i_dof % Shape/basis function.
case 1
dNdxi1 = (3125*xi(1)^4)/24 - (625*xi(1)^3)/3 + (875*xi(1)^2)/8 - (125*xi(1))/6 + 1;
case 2
dNdxi1 = - (3125*xi(1)^4)/24 + (625*xi(1)^3)/2 - (2125*xi(1)^2)/8 + (375*xi(1))/4 - 137/12;
case 3
dNdxi1 = - (15625*xi(1)^4)/24 + (6875*xi(1)^3)/6 - (5125*xi(1)^2)/8 + (1525*xi(1))/12 - 25/4;
case 4
dNdxi1 = (15625*xi(1)^4)/12 - 2500*xi(1)^3 + (6125*xi(1)^2)/4 - 325*xi(1) + 50/3;
case 5
dNdxi1 = - (15625*xi(1)^4)/12 + (8125*xi(1)^3)/3 - (7375*xi(1)^2)/4 + (2675*xi(1))/6 - 25;
case 6
dNdxi1 = (15625*xi(1)^4)/24 - (4375*xi(1)^3)/3 + (8875*xi(1)^2)/8 - (1925*xi(1))/6 + 25;
end
vBase = aInvJac(:,1) * dNdxi1;
otherwise
vBase = 0;
end
% end function sf_line_P5
```

Note that working in one dimension here it is enough to use one of
the barycentric coordinates in the basis function definitions, here
*xi(1)* (since *xi(2) = 1 - xi(1)*). Also, the derivative in local
coordinates *dNdxi1* is converted to global coordinates *dNdx* by
multiplication with the inverse of the grid cell transformation
Jacobian *aInvJac*, which in 1D is equal to the inverted grid cell
width.

In contrast to typical FEM codes where matrix assembly is performed
first in a loop over grid cells and then a loop over quadrature
points, the FEATool FEM assembly operations are vectorized and
performed in a single loop over quadrature points for groups of
several cells at the same time. Thus the output *vBase* is a row
vector with the same number of rows as *aInvJac*. However, as is done
here for *i_eval = 1* it is usually sufficient to compute and return a
scalar value when the shape function evaluation for a quadrature point
*xi* is independent of the real grid cell shapes. The assembly caller
function then assigns this scalar value to the corresponding output
vector entries.

To plot and visualize the defined shape functions we can use the script below

```
xi = linspace( 0, 1, 100 );
grid on, hold on
colours = {' r' 'g' 'b' 'c' 'm' 'y' };
for i_dof=1:6
for j=1:numel(xi)
[phi(i_dof,j),~,xLDof] = sf_line_P5( 1, 1, 2, i_dof, xi(j) );
end
plot( xi, phi(i_dof,:), colours{i_dof} )
text( xLDof(1,i_dof), 1, ['\phi_',num2str(i_dof),'(\xi)'] )
xlabel( '\xi' )
end
title( 'P_5 Lagrange FEM shape functions' )
set( gca, 'xdir', 'reverse' )
```

We can also confirm that the sum of the basis functions is equal to unity over the defined span 0…1

```
phi_sum = all( abs( sum(phi,1) - 1 ) < sqrt(eps) )
```

To test the new shape function we set up and solve a simple Poisson equation with a polynomial source term and exact analytical solution

\[
\left\{\begin{aligned}
\Delta u &= 140x^3 - 108x^2 + 18x - 4 \\\

u(x) &= -7x^5 + 9x^4 - 3x^3 + 2x^2 - 2x
\end{aligned}\right.
\]

```
fea.sdim = { 'x' };
fea.grid = linegrid( 1 );
fea.dvar = { 'u' };
fea.sfun = { 'sf_line_P5' };
% System matrix.
fea.eqn.a.form = { [2;2] };
fea.eqn.a.coef = { 1 };
% Source term.
fea.eqn.f.form = { 1 };
fea.eqn.f.coef = { '140*x^3 - 108*x^2 + 18*x - 4' };
% Boundary conditions.
fea.bdr.d = { 0 -1 };
fea.bdr.n = { [] [] };
% Solution.
fea.sol.u = solvestat( fea, 'icub', 5 );
```

In this case we have solved the equation on a single grid cell. We can
plot the solution and compare with the exact analytical one, and also
with linear and quadratic shape functions. (Note that the solver
parameter *icub* has been increased to *5* from the default *2* so
that the numerical integration can accurately integrate the 5th order
polynomials).

```
% Postprocessing.
x = linspace( 0, 1, 100 );
uref = evalexpr( '-7*x^5 + 9*x^4 - 3*x^3 + 2*x^2 - 2*x', x, parseprob(fea) );
u5 = evalexpr( 'u', x, parseprob(fea) );
fea.sfun = { 'sflag2' };
fea.sol.u = solvestat( fea );
u2 = evalexpr( 'u', x, parseprob(fea) );
fea.sfun = { 'sflag1' };
fea.sol.u = solvestat( fea );
u1 = evalexpr( 'u', x, parseprob(fea) );
grid on, hold on
plot( x, uref, 'k--', 'linewidth', 2 )
plot( x, u5, 'g-' )
plot( x, u2, 'b-' )
plot( x, u1, 'r-' )
hl = legend( 'exact solution', 'sf_line_P5', 'sflag2', 'sflag1' );
title( 'Solution to \Deltau = 140*x^3 - 108*x^2 + 18*x - 4' )
```

We can see that the $P_5$ shape function gives the exact solution with just a single grid cell while the $P_1$ and $P_2$ shape functions require more cells to achieve a better approximation.

It is also possible to use custom shape functions in the FEATool GUI
by entering the corresponding shape function name in the **FEM
Discretization** field of the **Equation Settings** dialog box (and
also increasing the **Numerical integration rule** order in the
**Solver Settings** dialog box).