comsol_evol_detailed_simulation.m

function sim_data = comsol_evol_detailed_simulation(model,model_data)
%function sim_data = comsol_evol_detailed_simulation(model,model_data)
%
% This function computes a detailed solution!
% Depending on model.use_comsol = 0 or 1 it uses COMSOL to calculate a
% solution or does it "by hand".
%
% Oliver Zeeb, 2013


sim_data.grid_info = model_data.grid_info;
if model.use_comsol
    comsol_model=comsol_set_mu_from_rbmatlab_model(model,model_data.comsol_model);
    comsol_model.sol(model.comsol_tags.sol).feature('t1').active(true);
    comsol_model.sol(model.comsol_tags.sol).runAll;
    sol_size = comsol_model.sol(model.comsol_tags.sol).getSize();
    nr_dofs = sol_size(1);
    nt = sol_size(2);
    sim_data.t =  comsol_model.sol(model.comsol_tags.sol).getPVals();
    sim_data.U = zeros(nr_dofs, nt);
    for k = 1:length(sim_data.t)
        sim_data.U(:,k) = comsol_model.sol(model.comsol_tags.sol).getU(k);
    end
    
    % calculation of the derivative --> calculate only if needed!
    % needs twice the time than storing only sim_data.U!!!
    %sim_data.Udot = zeros(nr_dofs, nt);
    %for k = 1:length(sim_data.t)
    %    sim_data.Udot(:,k) = comsol_model.sol(model.comsol_tags.sol).getUDot(k); 
    %end
    if model.compute_output_functional
        %sim_data.s = model.output_functional(model,sim_data);
        %get operators:
        v = model.operators_output(model,sim_data);
        sim_data.s = v'*sim_data.U;
    end
    
    
else %~model.use_comsol
    nt = model.nt;
    nr_DOFs = length(model_data.operators.init_values);
    U = zeros(nr_DOFs,nt+1); 
    
    % Initialize array for eliminated solution 
    Uc = zeros(size(model_data.ind_vectors.Null,2),nt+1); 

    model.decomp_mode = 1;
    [L_I_comp, L_E_comp, b_comp] = model.operators_ptr(model,model_data);

    model.decomp_mode = 2;
    [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model,model_data);

    L_I = lincomb_sequence(L_I_comp, L_I_coeff);
    L_E = lincomb_sequence(L_E_comp, L_E_coeff);    
    b = lincomb_sequence(b_comp, b_coeff);

% %%%%
% %FOR TESTING: ARE THE MATRICES OK?
% comsol_set_mu_from_rbmatlab_model(model,model_data.comsol_model);
% str = get_assembled_matrices(model_data.comsol_model);
% LI_comsol = str.Dc+model.dt*str.Kc;
% LE_comsol = str.Dc;
% b_comsol  = model.dt * str.Lc;
% %Differenz bilden
% LI_diff = L_I-LI_comsol;
% LE_diff = L_E-LE_comsol;
% b_diff = b-b_comsol;
% %maximum der  differenz ermitteln
% max_LI = max(max(abs(LI_diff)))
% max_LE = max(max(abs(LE_diff)))
% max_b  = max(abs(b_diff))
% %norm der differenz ermitteln
% norm(LI_diff,'fro')
% norm(LE_diff,'fro')
% norm(b_diff,'fro')
% keyboard
% %ENDE: ZUM TESTEN OB MATRITZEN PASSEN!
% %%%%



    % Time stepping iteration
    for t_ind = 1:nt;
        Uc(:,t_ind+1) = L_I \ (b + L_E*Uc(:,t_ind)) ;
    end
    
   % blow up the vectors to full size 
   for t_ind = 1:nt+1;
       U(:,t_ind) = model_data.ind_vectors.Null*Uc(:,t_ind);
   end

    sim_data.comsol_model = [];
    sim_data.U = U;
    
    if model.compute_output_functional
        %sim_data.s = model.output_functional(model,sim_data);
        %get operators:
        v = model.operators_output(model,model_data);
        sim_data.s = v'*sim_data.U;
    end
    
end