stokes_infsup_constant.m

function beta = stokes_infsup_constant(model, model_data)
%function beta = stokes_infsup_constant(model, model_data)
%

old_decomp_mode = model.decomp_mode;
model.decomp_mode = 0;

% assembly required matrix
S = Fem.Assembly.stokes_fd_nl_operator(model, model_data);

model.decomp_mode = old_decomp_mode;

% get inner prodcut
K = getMatrix(model.get_inner_product_matrix(model_data));
dir_gids = model_data.bc_info.dirichlet_gids;

if ~isempty(dir_gids)
    K(dir_gids, :) = [];
    K(:, dir_gids) = [];
    S(dir_gids, :) = [];
    S(:, dir_gids) = [];
end

K = 0.5*(K + K');

% compute constant using eigs
fun = @(X) S \ (K * (S' \ X));

opts = [];
opts.tol = 1e-9;
opts.issym = 1;

beta = sqrt(eigs(fun, size(S, 1), K, 1, 'SM', opts));

end