stokes_sobolev_embedding_constant.m

function rho_sqr = stokes_sobolev_embedding_constant(model, model_data, tol, maxiter)
%function rho_sqr = stokes_sobolev_embedding_constant(model, model_data, tol, maxiter)
%

% IM 02.12.2014

if model.verbose > 1
    disp('start sobolev embedding constant algorithm');
end
start_iteration = tic;

if nargin < 4
    maxiter = [];
end
if nargin < 3
    tol = [];
end

if isempty(tol)
    tol = 1e-9;
end
if isempty(maxiter)
    maxiter = 20;
end

vel_index = get_index(model_data.df_info, 'velocity');
df_info = get(model_data.df_info, vel_index);
v = Fem.DiscFunc([], df_info);

B = get(model.get_inner_product_matrix(model_data), vel_index) ...
    + eps*speye(df_info.ndofs);

dirichlet_gids = [];
if model.has_dirichlet_values
    dirichlet_gids = model_data.bc_info.dirichlet_gids;
    B(dirichlet_gids, :) = [];
    B(:, dirichlet_gids) = [];
end
B = 0.5*(B+B');

int_dofs = setdiff(1:df_info.ndofs, dirichlet_gids);

tmodel = [];
tmodel.reaction = @(grid, eindices, loc, params) ones(length(eindices), 1);
tmodel.qdeg = 2*model.pdeg;

rho_sqr = 0;
rho_old = -1;
niter = 0;

while abs(rho_old - rho_sqr) > tol
    
    rho_old = rho_sqr;
    
    A = Fem.Assembly.matrix_part(@Fem.IntegralKernels.matrix_reaction, ...
                                        tmodel, df_info) + eps*speye(df_info.ndofs);
    if ~isempty(dirichlet_gids)
        A(dirichlet_gids, :) = [];
        A(:, dirichlet_gids) = [];
    end
    
    [v_int_dofs, rho_sqr] = eigs(A, B, 1, 'LM');
    niter = niter+1;
    if model.verbose > 1
        disp(['deviation to previous step: ', num2str(abs(rho_old - rho_sqr))]);
    end
    
    % compute factor
    v.dofs(int_dofs) = v_int_dofs;    
    v1 = get_component(v, 1);
    v2 = get_component(v, 2);
    v1.dofs = v1.dofs.^2 + v2.dofs.^2;
    factor = 1/l2_norm(v1);
    
    tmodel.reaction = @(grid, eindices, loc, params) ...
        factor * sum(evaluate(v, eindices, loc, params).^2, 2);
    
    if niter >= maxiter
        warning('stokes:sobolev_embedding_constant', ...
            ['iteration maximum reached, norm(defect): ', num2str(abs(rho_old - rho_sqr))])
        break
    end
end
if model.verbose > 1
    disp(['terminate sobolev embedding constant algorithm after ', ...
        num2str(toc(start_iteration)), ' seconds']);
end

end