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.