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 and Octave 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 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 (), 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 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

```
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 shape function gives the exact solution with just a single grid cell while the and 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).