fv_two_phase_imp.m

function [INC, J] = fv_two_phase_imp(model, model_data, S, CU, NU_ind)
%function INC = fv_two_phase_imp(model, model_data, S, CU, NU_ind)

  grid = model_data.grid;

  if nargin <= 4 || isempty(NU_ind)
    NU_ind = cell(1,3);
    NU_ind{1} = (1:grid.nelements)';
    NU_ind{2} = (1:model_data.gn_edges)';
    NU_ind{3} = (1:grid.nelements)';
  end

  Slength = grid.nelements;
  Ulength = model_data.gn_edges;
  U = CU(1:Ulength);
  P = CU(Ulength+1:end);

  %  model.diffusivity_ptr            = @(glob, U, model) struct('K', 1, 'epsilon', 0);%model.two_phase.phi(glob, U, model);
  %  model.diffusivity_derivative_ptr = @(glob, U, model) struct('K', 0, 'epsilon', 0);%model.two_phase.phi_derivative(glob, U, model);
  %  model.laplacian_ptr              = @(glob, U, model) -getfield(model.two_phase.Phi(glob, U, model), 'K')';
  %  model.laplacian_derivative_ptr   = @(glob, U, model) -getfield(model.two_phase.phi(glob, U, model), 'K');

  num_diff_flux_u = fv_num_diff_flux_pressure_gradient(model, model_data, P, S, NU_ind{3}, NU_ind{1});

  tmppar.local_mass_matrix  = @fv_local_mass_matrix_rect;
  tmppar.element_quadrature = @cubequadrature;
  tmppar.evaluate_basis     = @fv_evaluate_basis;
  tmppar.pdeg = 0;
  S_upper = grid.A(NU_ind{1}).*l2project(@(elemin, loc, grid) model.two_phase.upper_source(elemin, loc, grid), 2, model_data.grid, tmppar);
  S_lower = grid.A(NU_ind{1}).*l2project(@(elemin, loc, grid) model.two_phase.lower_source(elemin, loc, grid), 2, model_data.grid, tmppar);

%  S_upper = grid.A(NU_ind{1}) .* model.two_phase.upper_source([grid.CX(NU_ind{1}), grid.CY(NU_ind{1})], model);
%  S_lower = grid.A(NU_ind{1}) .* model.two_phase.lower_source([grid.CX(NU_ind{1}), grid.CY(NU_ind{1})], model);

  INC = zeros(Slength+Ulength+1, 1);

  UKL = -num_diff_flux_u.G(NU_ind{2});

  INC(1:Ulength) = UKL - U;

  % eventuell + S_upper - S_lower
  INC(Ulength+1:end-1) = sum(U(model_data.gEI).*(grid.NX+grid.NY), 2) ...
    - S_upper + S_lower;

  INC(end) = P(1);
%  INC(end) = sum(grid.A(NU_ind{3}, :) .* P);

  if nargout > 1

    u_diff_p_s = fv_frechet_operators_diff_flux_pressure_gradient(model, model_data, P, S, NU_ind{3}, NU_ind{1});

    n = Slength+Ulength+1;
    m = Slength+Ulength;

    sp2 = -sparse(u_diff_p_s.pi, u_diff_p_s.pj+Ulength, u_diff_p_s.pvals, n, m);

    sp3i = repmat(1:Slength, 1, 4)+Ulength;
    sp3j = model_data.gEI(:);
    sp3vals = grid.NX+grid.NY;
    sp3 = sparse(sp3i,sp3j,sp3vals, n, m);

    sp4 = sparse(Slength+Ulength+1, Ulength+1, 1, n, m);
%    sp4 = sparse(Slength+Ulength+1*ones(1,Slength), Ulength+1:Slength+Ulength, grid.A(:), n, m);

    sp5 = sparse((1:Ulength), (1:Ulength), -1, n,m);

    J = sp2 + sp3 + sp4 + sp5;
    J(isnan(J)) = 0;
  end

end