evol_TB_2D_weak_form_operators_output.m

function [v, s_l2norm] = evol_TB_2D_weak_form_operators_output(model,model_data)
%function [v, s_l2norm] = evol_TB_2D_weak_form_operators_output(model,model_data)
%
% This function returns the output operator v for the Comsol 2D evolution
% Thermal Block model.
%
% The output is defined as the integral over the bottom of the Thermal
% block:
% s(\mu) = l(u(\mu)) = \int_0^1 u(x,0) dx
% and can be comupted s(\mu) = v' * (rb_)sim_data.U;
%
% This function checks whether linear or quadratic finfite elements are
% used and uses a Simpson-rule (pdeg=2) or an trapezoidal rule (pdeg=1)
% If other Finite Elements are used, a trapezoidal rule is used (BUT THIS
% CASE HAS NOT BEEN TESTED, SO BE CAREFUL!!!)
%
% The L2-Norm of the output operator is the square root of the integration
% area, so in this case should be equal to 1.
%
% Input:
%   - model
%   - model_data
%
% Output:
%   - v        : output_operator: s=v'*sim-data.U;
%   - s_l2norm : L2-Norm of the ouput operator (can be used for error
%                estimation)
%
% Oliver Zeeb, 2013/11/21



if model.use_comsol
        grid_info = model_data.grid_info;

        % get dof indices of the bottom
        bnd_dof_ind = find(grid_info.dofs.coords(2,:) <= 10*eps);

        %get the according coords and values
        bnd_coords = grid_info.dofs.coords(1,bnd_dof_ind);

        %resort for quadrature
        [bnd_coord_sort,bnd_coord_sort_ind] = sort(bnd_coords); 
        %bnd_coord_sort is now: 
        %for pdeg=1: x1 x2 x3 x4 ...
        %for pdeg=2: x1 (x1+x2)/2 x2 (x2+x3)/2 x3 ... --> coords in between of the mesh edges!

        if strcmp(model.element_type{1},'shlag') %check lagrange finite elements are used
            v_dof = zeros(1,size(bnd_coords,2));
            switch model.pdeg{1}
                case 2 %fuer pdeg = 2: Simpson-Rule
                    v_dof(1) = 1/6*(bnd_coord_sort(3)-bnd_coord_sort(1));
                    v_dof(2:2:end-1) = diff(bnd_coord_sort(1:2:end))*2/3;
                    v_dof(3:2:end-2) = 1/6* (bnd_coord_sort(5:2:end) - bnd_coord_sort(1:2:end-3));
                    v_dof(end) = 1/6*(bnd_coord_sort(end)-bnd_coord_sort(end-2));
                otherwise %pdeg = 1 or anything else: trapezoidal rule
                    v_dof(1) = 0.5*(bnd_coord_sort(2) - bnd_coord_sort(1));
                    v_dof(2:end-1) = 0.5 * (bnd_coord_sort(3:end) - bnd_coord_sort(1:end-2));
                    v_dof(end) = 0.5 * (bnd_coord_sort(end) - bnd_coord_sort(end-1));
            end
        else %no lagrangian elements!
            warning('using trapezodial rule for output integration! Might be inaccurate since used element type is unknown!')
            warning('THIS CASE HAS NOT BEEN TESTED EXTENSEVELY SO FAR...')
            v_dof(1) = 0.5*(bnd_coord_sort(2) - bnd_coord_sort(1));
            v_dof(2:end-1) = 0.5 * (bnd_coord_sort(3:end) - bnd_coord_sort(1:end-2));
            v_dof(end) = 0.5 * (bnd_coord_sort(end) - bnd_coord_sort(end-1));
        end

        %Blow up the Vector again to be able to calculate v' * sim_data.U
        v = zeros(grid_info.ndofs,1);
        bnd_dof_ind_sort = bnd_dof_ind(bnd_coord_sort_ind); %since v is sorted, the indices of the bnd_dofs have to be sorted accordingly
        %now blow up the vector:
        v(bnd_dof_ind_sort) = v_dof;
        s_l2norm = sqrt(max(bnd_coords)-min(bnd_coords));


else %~model.use_comsol
    v = model_data.operators.operators_output;
    s_l2norm = model_data.operators.output_operator_l2norm;
end