refine_t_part.m

function [model, detailed_data] = refine_t_part(model, detailed_data, refine_p, sim_data)
%function refine_t_part(model, detailed_data, refine_p)
%
% function refines the t_partitions with numbers 'refine_p'.
% The function identifies the time step where half of the error growth in
% this partition is reached and refines the partition at this time step.
% refine_p is a vector of t-partition ids.
%
% parameters:
%   refine_p:       array of partition indices which are to be refined.
%   sim_data:       a maximum error curve over a test sample of mus
%
% Required fields of model:
%   t_part_range:   cell-array, giving information about the t-partitions
%   transition_model: boolean flag indicating whether a projection should take
%                     place at every overlap (false), of whether at an overlap
%                     both reduced basis shall be used (true). In the second
%                     case the t_part_ranges do not overlap.
%   
% Required fields of detailed_data:
%   t_part_detailed_data: the t-partition detailed data
%
%

% Markus Dihlmann 22.02.10

if ~isfield(model,'t_part_middle_point_bisection')
    model.t_part_middle_point_bisection =0;
end

t_part_detailed_data         = detailed_data.t_part_detailed_data;
t_parts_for_basis_generation = [];
t_part_range                 = model.t_part_range;
t_part_range_old             = t_part_range;
offset=0;
for p = 1:length(refine_p)
    
    % get detailed data for current t-partition
    tp_detailed_data = detailed_data.t_part_detailed_data{refine_p(p)};
    
    %identify for which mu the maximum error has been detected
    %mu=tp_detailed_data.RB_info.mu_sequence(:,end);
    
    %old_mu = get_mu(model);
    %t_part_old = model.t_part_for_simulation;
    %prepare a reduced simulation for this parameter
    %model = set_mu(model, mu);
    %reduced_data = gen_reduced_data(model, detailed_data);
    %model.t_part_for_simulation = refine_p(p);
    
    %sim_data = rb_simulation(model, reduced_data);
    
    %find k where the average error is reached
    if ~model.t_part_middle_point_bisection
        %Bisection at timestep k where half the error is reached
        delta_diff = sim_data.t_part_sim_data{refine_p(p)}.Delta(end)...
                     -sim_data.t_part_sim_data{refine_p(p)}.Delta(1);

        v          = find(sim_data.t_part_sim_data{refine_p(p)}.Delta>(delta_diff/2));
        k_cut      = v(1);
    else
       %bisection in the middle of the interval
       interval_length=t_part_range{refine_p(p)+offset}(2)-t_part_range{refine_p(p)+offset}(1);
       k_cut = floor(interval_length/2);
        
    end
    
    % check wether t-partition can be refined at k_cut
    if (k_cut>3) && (k_cut<(length(sim_data.t_part_sim_data{refine_p(p)}.Delta)-3))

        % new t partition interval
        new_t_partition                     = [k_cut + t_part_range{refine_p(p) + offset}(1), ...
                                               t_part_range{refine_p(p) + offset}(2)];
        % 
        t_part_range{refine_p(p)+offset}(2) = new_t_partition(1);
        % add new t partition interval to old range list
        t_part_range                        = add_to_array(new_t_partition,...
                                                           refine_p(p) + 1 + offset,...
                                                           t_part_range);
        % also copy the detailed data with RB to this new interval
        t_part_detailed_data                = add_to_array(tp_detailed_data,...
                                                           refine_p(p)+1+offset,...
                                                           t_part_detailed_data);
        % 
        t_parts_for_basis_generation        = [t_parts_for_basis_generation,...
                                               refine_p(p)+offset,...
                                               refine_p(p)+1+offset];
        offset = offset+1;
    else
        disp('no segmentation possible because of too big projection error');
    end
    
%    model = set_mu(model, old_mu);
%    model.t_part_for_simulation = t_part_old;
end

model.t_part_range                 = t_part_range;
detailed_data.t_part_detailed_data = t_part_detailed_data;
detailed_data.t_part_range         = t_part_range;
model.t_parts_for_basis_generation = t_parts_for_basis_generation;

%store partition history
if ~isfield(detailed_data,'t_part_ranges_history')
    detailed_data.t_part_ranges_history{1} = t_part_range_old;
else
    len = length(detailed_data.t_part_ranges_history);
    detailed_data.t_part_ranges_history{len+1} = t_part_range_old;
end

%store max error history
sim_data.Delta=sim_data.t_part_sim_data{1}.Delta;
if isfield(sim_data.t_part_sim_data{1},'s')
  sim_data.s       = sim_data.t_part_sim_data{1}.s;
  sim_data.Delta_s = sim_data.t_part_sim_data{1}.Delta_s;
end

for ind = 2:length(sim_data.t_part_sim_data)
  len = length(sim_data.t_part_sim_data{ind}.a(1,:));
  %sim_data.a=[sim_data.a,t_part_sim_data{ind}.a(:,2:len)];
  if ~model.transition_model
    sim_data.Delta=[sim_data.Delta,...
                    sim_data.t_part_sim_data{ind}.Delta(2:len)];
    if isfield(sim_data.t_part_sim_data{1},'s')
      sim_data.s       = [sim_data.s,...
                          sim_data.t_part_sim_data{ind}.s(2:len)];
      sim_data.Delta_s = [sim_data.Delta_s,...
                          sim_data.t_part_sim_data{ind}.Delta_s(2:len)];
    end
  else
    sim_data.Delta=[sim_data.Delta,...
                    sim_data.t_part_sim_data{ind}.Delta(1:len)];
    if isfield(sim_data.t_part_sim_data{1},'s')
      sim_data.s       = [sim_data.s,sim_data.t_part_sim_data{ind}.s(1:len)];
      sim_data.Delta_s = [sim_data.Delta_s,sim_data.t_part_sim_data{ind}.Delta_s(1:len)];
    end
  end
end

if ~isfield(detailed_data,'max_error_history')
    detailed_data.max_error_history{1} = sim_data.Delta;
else
    len                                    = length(detailed_data.max_error_history);
    detailed_data.max_error_history{len+1} = sim_data.Delta;
end