gen_model_data.m

function model_data = gen_model_data(dmodel)
% function model_data = gen_model_data(model)
% method which produces all `H`-dependent data, that is required for a detailed
% simulations and is independent of the parameter `\mu`.
%
% Generated fields of model_data:
%  W      : inner product weighting matrix

ddescr = dmodel.descr;

grid = construct_grid(ddescr);
model_data.grid = grid;

model_data.W    = ddescr.inner_product_matrix_algorithm(ddescr, model_data);

[model_data.gEI, model_data.gn_inner_edges, model_data.gn_boundary_edges, ...
    i_ints_bigger, b_ints] ...
  = compute_edge_indices(grid);
model_data.gn_edges = model_data.gn_inner_edges + model_data.gn_boundary_edges;

dWdiag = zeros(1, model_data.gn_edges);
relevant_edges = [i_ints_bigger; find(b_ints)];
%dWdiag(model_data.gEI(relevant_edges)) = 1/2.*grid.DS(relevant_edges).*grid.EL(relevant_edges);
dWdiag(model_data.gEI(relevant_edges)) = grid.EL(relevant_edges);
model_data.diamondW = spdiags(dWdiag(:), 0, model_data.gn_edges, model_data.gn_edges);
model_data.diamondWinv = spdiags(1./dWdiag(:), 0, model_data.gn_edges, model_data.gn_edges);

tmppar.local_mass_matrix  = @fv_local_mass_matrix_rect;
tmppar.element_quadrature = @cubequadrature;
tmppar.evaluate_basis     = @fv_evaluate_basis;
tmppar.pdeg = 0;
model_data.opdata = [];
model_data.opdata.S_upper = grid.A .* ...
  l2project(@(elemin, loc, grid) ddescr.two_phase.upper_source(elemin, loc, grid), ...
            2, model_data.grid, tmppar);
model_data.opdata.S_lower = grid.A .* ...
   l2project(@(elemin, loc, grid) ddescr.two_phase.lower_source(elemin, loc, grid), ...
            2, model_data.grid, tmppar);