matrix_neumann.m

function res = matrix_neumann(x, model, df_info, i, j, elind, edgeind)
%function res = matrix_neumann(x, model, df_info, i ,j, elind, edgeind)
%
% auxiliary function for integral kernel for A_neumann
% integral kernel on all neumann-edges simultaneously
% here x is a scalar value on unit interval
% f(x) = (b^T n) hatphi_i hatphi_j 
% multiplication with |edge| is realized in caller after quadrature.
% function can handle cell-array valued data

% B. Haasdonk 22.2.2011
% IM 21.07.2016

grid = df_info.grid;
res = zeros(0,1);
for local_edge_id = 1:3
  % transform 1d coordinate to 2d local coordinate
  lcoord = llocal2local(grid,local_edge_id,x); % local coord for x on edge 1    
  hat_phi_i = evaluate_basis_function(df_info,lcoord,i);
  hat_phi_j = evaluate_basis_function(df_info,lcoord,j);
  elinds = elind(edgeind==local_edge_id);
  b = cache_function(@(xelind, xarg, xmodel)xmodel.velocity(grid, ...
      xelind, xarg, xmodel), elinds, lcoord, model);
  if ~iscell(b)
    res = [res;...
           ( b(:,1).*grid.NX(elinds,local_edge_id) ...
             + b(:,2).*grid.NY(elinds,local_edge_id) ...
             ) * hat_phi_i * hat_phi_j; ...
          ];
  else % iscell!!
    if ~iscell(res)
      res = cell(1,length(b));
      for q = 1:length(b)
        res{q} = zeros(0,1);
      end;
    end;
    for q = 1:length(b)      
      res{q} = [res{q};...
             ( b{q}(:,1).*grid.NX(elinds,local_edge_id) ...
             + b{q}(:,2).*grid.NY(elinds,local_edge_id) ...
             ) * hat_phi_i * hat_phi_j; ...
          ];
    end;
  
  end;
end;