FEATool solver hook example for implementing the (summation) constraint c1 = c2 over a subdomain

featool_solver_hook_example.m

function [ fea ] = featool_solver_hook_example()
%FEATOOL_SOLVER_HOOK_EXAMPLE
%
% Solver hook example to add a constraint that sum(ci_j) = 0 for two
% dependent variables c1 and c2 (and j is the nodes/degrees of freedom).

% Copyright 2021 Precise Simulation, Ltd.

% Geometry and grid.
fea.sdim = { 'x', 'y' };
fea.geom.objects = { gobj_rectangle() };
fea.grid = rectgrid(10);

% Equation settings.
fea = addphys( fea, @convectiondiffusion, { 'c1' } );
fea.phys.cd.eqn.coef = { 'dts_cd', 'd_ts', 'Time scaling coefficient', { '1' };
                         'd_cd', 'D', 'Diffusion coefficient', { '1' };
                         'u_cd', 'u', 'Convection velocity in x-direction', { '0' };
                         'v_cd', 'v', 'Convection velocity in y-direction', { '0' };
                         'r_cd', 'R', 'Reaction rate', { '1' };
                         'c10_cd', 'c1_0', 'Initial condition for c1', { '0' } };

fea = addphys( fea, @convectiondiffusion, { 'c2' } );
fea.phys.cd2.eqn.coef = { 'dts_cd2', 'd_ts', 'Time scaling coefficient', { '1' };
                          'd_cd2', 'D', 'Diffusion coefficient', { '1' };
                          'u_cd2', 'u', 'Convection velocity in x-direction', { '0' };
                          'v_cd2', 'v', 'Convection velocity in y-direction', { '0' };
                          'r_cd2', 'R', 'Reaction rate', { '0' };
                          'c20_cd2', 'c2_0', 'Initial condition for c2', { '0' } };

% Boundary settings.
fea.phys.cd.bdr.sel = [ 1, 1, 2, 2 ];
fea.phys.cd.bdr.coef = { 'bcr_cd', 'c1 = c1_0', 'Concentration', { 'c1_0' }, { 1, 1, 1, 1 }, [], { '0', '0', '0', '0' };
                         'bcc_cd', 'n.(-Dgrad c1) = 0', 'Convective flux/outflow', [], { 0, 0, 0, 0 }, [], { 0, 0, 0, 0 };
                         'bci_cd', 'n.(-Dgrad c1 + uc1) = 0', 'Insulation/symmetry', [], { 0, 0, 0, 0 }, { '(nx*u_cd+ny*v_cd)*c1' }, { '0', '0', '0', '0' };
                         'bcf_cd', '-n.(-Dgrad c1 + uc1) = N_0', 'Flux condition', { 'N_0' }, { 0, 0, 0, 0 }, { '(nx*u_cd+ny*v_cd)*c1' }, { '0', '0', '0', '0' } };

fea.phys.cd2.bdr.sel = [ 2, 2, 1, 1 ];
fea.phys.cd2.bdr.coef = { 'bcr_cd2', 'c2 = c2_0', 'Concentration', { 'c2_0' }, { 1, 1, 1, 1 }, [], { 0, 0, '1-x', 'y' };
                          'bcc_cd2', 'n.(-Dgrad c2) = 0', 'Convective flux/outflow', [], { 0, 0, 0, 0 }, [], { 0, 0, 0, 0 };
                          'bci_cd2', 'n.(-Dgrad c2 + uc2) = 0', 'Insulation/symmetry', [], { 0, 0, 0, 0 }, { '(nx*u_cd2+ny*v_cd2)*c2' }, { '0', '0', '0', '0' };
                          'bcf_cd2', '-n.(-Dgrad c2 + uc2) = N_0', 'Flux condition', { 'N_0' }, { 0, 0, 0, 0 }, { '(nx*u_cd2+ny*v_cd2)*c2' }, { '0', '0', '0', '0' } };


% Define solver hook to call before/after solving linear system.
fea.sol.settings.hook = @my_solve_hook;


% Solver call.
fea = parsephys( fea );
fea = parseprob( fea );
fea.sol.u = solvestat( fea );


% Postprocessing.
subplot(1,3,1)
postplot( fea, 'surfexpr', 'c1', 'title', 'Concentration, c1' );
subplot(1,3,2)
postplot( fea, 'surfexpr', 'c2', 'title', 'Concentration: c2' );
subplot(1,3,3)
postplot( fea, 'surfexpr', 'c1+c2', 'title', 'c1+c2' );


% Check results.
ind_c1 = find(strcmp(fea.dvar, 'c1'));
offset = sum(fea.eqn.ndof(1:ind_c1-1));
ind_dof_c1 = unique(fea.eqn.dofm{ind_c1}(:)) + offset;

ind_c2 = find(strcmp(fea.dvar, 'c2'));
offset = sum(fea.eqn.ndof(1:ind_c2-1));
ind_dof_c2 = unique(fea.eqn.dofm{ind_c2}(:)) + offset;

c1_plus_c2 = sum(fea.sol.u(ind_dof_c1)) + ...
             sum(fea.sol.u(ind_dof_c2))
assert( abs(c1_plus_c2) < 1e-12 )


%------------------------------------------------------------------------------%
function [ A, f, prob ] = my_solve_hook( A, f, prob, solve_step )
%MY_SOLVE_HOOK Custom solver hook calling function.
%
%   [ A, F, PROB ] = SOLVE_HOOK( A, F, PROB, SOLVE_STEP )
%   Special solver hook function. Allows for special functions to be
%   called before and after the linear solver to modify the global
%   matrix A and right hand side/load vector F with the information in
%   the finite element problem struct PROB. SOLVE_STEP indicates if
%   the function call is before the linear solver call (=1), or to
%   clean-up after (=2).
%
%   A solve hook may be prescribed by entering a function name string
%   or handle in the fea.sol.settings.hook field. The corresponding
%   function will in each solve iteration then be called as
%
%       [A,f,fea] = fcn( A, f, fea, solve_step );

if( solve_step == 1 )
  disp('Calling solver hook BEFORE linear solver call.')

  ind_c1 = find(strcmp(prob.dvar, 'c1'));
  offset = sum(prob.eqn.ndof(1:ind_c1-1));
  ind_dof_c1 = unique(prob.eqn.dofm{ind_c1}(:)) + offset;
  n = length(ind_dof_c1);

  ind_c2 = find(strcmp(prob.dvar, 'c2'));
  offset = sum(prob.eqn.ndof(1:ind_c2-1));
  ind_dof_c2 = unique(prob.eqn.dofm{ind_c2}(:)) + offset;

  if( length(ind_dof_c1) ~= length(ind_dof_c2) )
    error('c1 and c2 must have the same discretization (fea.sfun)')
  end

  if( isstruct(A) )     % Sparsify matrix if necessary.
    A = sparse( A.rows, A.cols, A.vals, A.size(1), A.size(2) );
  end

  % Build constraint contributions.
  B = sparse(1, size(A,2));
  B(sub2ind(size(B), ones(1,n), ind_dof_c1(:)')) = 1;
  B(sub2ind(size(B), ones(1,n), ind_dof_c2(:)')) = 1;

  % Add row to add constraint sum_i(c1_i + c2_i) = 0.
  A = [A, B';
       B, 0 ];
  f = [f; 0];
end

if( solve_step == 2 )
  disp('Calling solver hook AFTER linear solver call.')

  % Remove system modifications.
  n0 = sum(prob.eqn.ndof);
  A = A(1:n0,1:n0);
  f = f(1:n0);

end