This post describes how to implement finite element FEM models with custom periodic boundary conditions in FEATool. A periodic boundary condition can be defined for opposing boundaries so that their values are linked in some defined way. The typical case for two periodic boundaries is to require them to have identical values, thus representing a partly infinite domain. They can also be used to offset boundary values by a constant, such as for flow boundaries with a given pressure differential.

As periodic boundary conditions are currently not part of the pre-defined physics modes, they have to be implemented by custom modifications. This can be accomplished with the solver hook functionality described below.

## Solve Hook Functionality

The solver hook functionality is a part of the
solvestat and
solvetime functions to allow for custom
matrix and right hand side vector modifications. The functionality can
currently only be activated for Neumann boundary conditions defined
with a string coefficient named `solve_hook_myfun`

(where
*myfun* can be any custom function name). This is also applicable to
the MATLAB GUI for physics modes where pure Neumann boundary
coefficients can be prescribed. If present, the solver will discard
the *solve_hook_* prefix and try to call the function
*myfun* both before and after the linear solver calls (```
u =
A\b
```

). The solve hook function header and calls take the form

```
function [ A, f, fea ] = myfun( A, f, fea, i_dvar, j_bdr, solve_step )
```

where *A* and *f* are the assembled system matrix and right hand
side/load vector, respectively. *fea* is the finite element problem
struct, *i_dvar* the current dependent variable, *j_bdr* the
boundary. *solve_step* is an integer flag which indicates if the
call is before the linear solve step (1), or after it (2). The
(optionally) modified outputs are returned in *A*, *f*, and *fea*. In
pseudo code the solve hook calling step looks something like the
following

```
...
if( exist('myfun') )
[A,f,fea] = myfun( A, f, fea, i_dvar, j_bdr, 1 );
end
u = A\f;
if( exist('myfun') )
[A,f,fea] = myfun( A, f, fea, i_dvar, j_bdr, 2 );
end
...
```

Besides being used to set and apply periodic boundary conditions as explained below, the solve hook functionality is also automatically used to set velocity slip boundary conditions (non-aligned with the coordinate system) for Navier-Stokes and Brinkman equation fluid flow physics modes, and can also to couple dependent variables between different active subdomains and physics modes.

## Moving Pulse with Periodic Boundary Conditions

This example shows how to use the solve hook functionality to define a
periodic pulse that is moving out from one side of a one dimensional
domain and reentering the other. The first step to define the model on
the command line is to assign a name or label for the space dimension
(here *x*), and generate a 100 cell grid for the unit line

```
fea.sdim = {'x'};
fea.grid = linegrid( 100 );
```

Next, the *convection and diffusion physics* mode is added with a
diffusion coefficient of *1*, and a constant convection velocity of
*0.002* in the forward *x*-direction

```
fea = addphys( fea, @convectiondiffusion );
fea.phys.cd.eqn.coef{3,end} = { 1 };
fea.phys.cd.eqn.coef{2,end} = { 0.002 };
fea = parsephys( fea );
```

The *parsephys* call processes the physics mode and translates the
information to the global finite element *fea.eqn* and *fea.bdr*
fields.

Proceeding, any prescribed Dirichlet boundary conditions are cleared
(which if present take priority over and overrides Neumann flux
boundary conditions), and a *solve_hook* function, here named
`periodic_bc`

, is prescribed for the right node (boundary
2).

```
i_dvar = 1; % Dependent variable number.
j_bdr = 2; % Boundary selection.
fea.bdr.d{i_dvar} = { [] [] };
fea.bdr.n{i_dvar}{j_bdr} = 'solve_hook_periodic_bc';
```

Before the problem can be solved, the actual *periodic_bc*
function must first be created and defined. To do this, create a
separate file named `periodic_bc.m`

with the following
contents

```
function [ A, f, prob ] = periodic_bc( A, f, prob, i_dvar, j_bdr, solve_step )
% Only process step before the linear solver call.
if( solve_step~=1 ), return, end
% Sparsify matrix if necessary.
if( isstruct(A) ), A = sparse( A.rows, A.cols, A.vals, A.size(1), A.size(2) ); end
% Degrees of freedom to couple.
dof_i = 1; % Left side node/dof.
dof_j = prob.eqn.ndof; % Right side node/dof.
% Add i and j rows (to preserve both i and j equations).
A(dof_i,:) = A(dof_i,:) + A(dof_j,:);
f(dof_i) = f(dof_i) + f(dof_j);
% Modify j rows so that A(j,i) - A(j,j) = 0.
A(dof_j,:) = 0; % Zero out periodic bc rows.
A(dof_j,[dof_i dof_j]) = [ 1 -1 ];
f(dof_j) = 0;
```

The *periodic_bc* function is here customized for this exact
problem to couple the first (left side) node (degree of freedom) with
the last right side one. In essence this modifies the *A* matrix so
that the relation *u(1) - u(100) = 0* is enforced.

Finally, assuming the *periodic_bc* function exists and can be
found in any of the the MATLAB paths, the problem can now be checked
for error, parsed, and solved with a unit pulse defined between
*x=0.4* and *x=0.6*

```
fea = parseprob( fea );
fea.sol.u = solvetime( fea, 'init', '(x>=0.4)*(x<=0.6)', ...
'tstep', 0.005, 'tmax', 1 );
```

After the solver has finished, the solution at different times can be plotted and visualized with the following code

```
postplot( fea, 'surfexpr', 'c', 'solnum', 1, 'color', 'b' )
postplot( fea, 'surfexpr', 'c', 'solnum', floor(size(fea.sol.u,2)/4), 'color', 'c' )
postplot( fea, 'surfexpr', 'c', 'solnum', floor(size(fea.sol.u,2)/2), 'color', 'm' )
postplot( fea, 'surfexpr', 'c', 'solnum', floor(size(fea.sol.u,2)*3/4), 'color', 'g' )
postplot( fea, 'surfexpr', 'c', 'solnum', size(fea.sol.u,2) )
axis( [0 1 0 1] )
grid on
```

Alternatively, snapshots and an animated gif movie of the moving
pulse can be created with these commands (assuming that the *convert*
tool from ImageMagick is installed)

```
x = linspace(0,1,101);
for i=1:size(fea.sol.u,2)
y = fea.sol.u(:,i); clf
patch( 'vertices', [[0;x(:);1] [0;y(:);0]], 'faces', [1:103 1], ...
'facecolor', [0.5 0.5 0.8], 'linewidth', 3 );
set( gca, 'fontsize', 20, 'ytick', [], 'ycolor', [1 1 1] )
axis( [0 1 0 1.05] )
print( '-r200', '-dpng', sprintf( 'img%03i.jpg', i ) )
end
system( ['convert -loop 0 img*.png -resize 640x480 -dither none' ...
'-colors 8 -fuzz "40%" -layers OptimizeFrame output.gif'] )
```

As can be seen the solution depicts an initially square pulse being smoothed due to diffusion while moving to the right. As it reaches the right hand side edge it reenters the left side as expected from a periodic domain.

## Other Periodic Boundary Condition Examples

For more examples defining and using periodic boundary the conditions,
see the axisymmetric Taylor-Couette swirl flow model, and the two dimensional
periodic Poisson equation example which is available in the FEATool
model and *examples* directory as the ex_periodic2
MATLAB script file.