disc_gradient.m

function [gradx,grady] = disc_gradient(glob,fdf)
%function [gradx,grady] = disc_gradient(glob,fdf)
%
% compute the gradient of the femdiscfunc fdf
% 
% input:
% glob: list of coordinates of points, in which the gradient shall be
%       calculated. Must not be points in grid vertices or edges, but
%       inside elements.
% fdf: femdiscfunc (FEM solution of a PDE)
% 
% output:
% gradx : values of the gradient in x direction
% grady : values of the gradient in y direction
%
% this computation is only exact for pdeg = 1!
%
% Immanuel Maier, 07.02.2011, Oliver Zeeb, 10.02.11

elind = zeros(size(glob,1),1);

% loop over all elements:
for i = 1:fdf.grid.nelements
  % compute local coordinates:
  loc = global2local(fdf.grid,i,glob);
  % if points lie in element i,
  % save element-index in elind:
  iselement = find((loc(:,1) >= 0-5000*eps) & (loc(:,2) >= 0-5000*eps) & ...
                   (loc(:,1)+loc(:,2) <= 1+5000*eps)); 
  elind(iselement) = i;
end;

if ~isempty(find(elind == 0))
  error('cannot find elements for all points!');
end;

% compute grad_hat u_hat
grad_hatx = fdf.dofs(fdf.grid.VI(elind,2)) - fdf.dofs(fdf.grid.VI(elind,1));
grad_haty = fdf.dofs(fdf.grid.VI(elind,3)) - fdf.dofs(fdf.grid.VI(elind,1));

% compute grad u with coordinate-transformation-derivative-rule
gradx = fdf.grid.JIT(elind,1,1) .* grad_hatx ...
        + fdf.grid.JIT(elind,1,2) .* grad_haty;
grady = fdf.grid.JIT(elind,2,1) .* grad_hatx ...
        + fdf.grid.JIT(elind,2,2) .* grad_haty;