dom_dec_gen_model_data.m

function model_data = dom_dec_gen_model_data(model)
% function model_data = dom_dec_model_data(model)
%
% function computing required data for the simulation of model
% currently only supporting polynomial degree 1
%
% Generated fields of model_data:
%  grids: cell containing the grids for both subdomains
%  df_infos: cell containing the infos for the discrete function
%  spaces on the subdomains
%  gamma_dofs: cell containing indices of the dofs on the interface
%  gamma in same order
%  masks: cell containing vectors for extracting a discrete
%  function on a subdomain out of a function on the whole domain
%  base_model_data: required data for the simulation of the
%  base_model

% I.Maier, 19.07.2011

model_data = [];
model_data.df_infos = cell(1,2);
model_data.gamma_dofs = cell(1,2);
model_data.masks = cell(1,2);
model_data.gamma_inner_product_matrices = cell(1,2);
model_data.mu_ref_extension_matrices = cell(1,2);

% full grid:
grid = construct_grid(model.base_model); %structure storing much
                                         %geometrical information

model_data.base_model_data = [];
model_data.base_model_data.df_info = feminfo(model.base_model, ...
                                             grid);

inds = cell(1,2);
% find indices of elements in omega_1 
for i = 1:length(model.dd_rect_corner1)
  ind_tmp = find( (grid.CX>model.dd_rect_corner1{i}(1)) & ...
                  (grid.CY>model.dd_rect_corner1{i}(2)) & ...
                  (grid.CX<model.dd_rect_corner2{i}(1)) & ...
                  (grid.CY<model.dd_rect_corner2{i}(2)) );
  inds{1} = [inds{1};ind_tmp];
end;

% indices of elements in omega_2
inds{2} = setdiff(1:grid.nelements,inds{1});

% subgrids:
grids = cell(1,2);
for i = 1:2
  grids{i} = gridpart(grid,inds{i});
  % structures containing information about the fe-function space
  model_data.df_infos{i} = feminfo(model.base_model,grids{i});
end;

% generate index vectors of vertices on gamma in same order
masks = cell(1,2);
for i = 1:2
  masks{i} = zeros(grid.nvertices,1);
  masks{i}(grid.VI(inds{i},:)) = 1;
  masks{i} = logical(masks{i});
end;
gamma_vind = (masks{1}==masks{2});

% supports only pdeg = 1!
for i = 1:2
  model_data.gamma_dofs{i} = ...
      setdiff(find(gamma_vind(masks{i})), ...
              model_data.df_infos{i}.dirichlet_gids);
end;

% the following is needed for error computation
for i = 1:2
  model_data.masks{i} = masks{i};
end;

% compute inner product matrices on gamma (L2 and H10)
no_dir = mod(model.dirichlet_side,2)+1;
[elind,edgeind] = find(grids{no_dir}.NBI == -10);

model_data.gamma_inner_product_matrices{1} = ...
    fem_l2_boundary_inner_product_matrix( ...
        model_data.df_infos{no_dir},elind,edgeind);

model_data.gamma_inner_product_matrices{2} = ...
    fem_h10_boundary_inner_product_matrix( ...
        model_data.df_infos{no_dir},elind,edgeind);

model = set_mu(model, model.mu_ref);
[dummy1, dummy2, model_data.mu_ref_extension_matrices, dummy3] = ...
    model.operators(model, model_data);