% FUNCTION [x] = gf_compute(meshfem MF, vec U, operation [, args])
% 
%   Various computations involving the solution U of the finite element problem.
% 
%   * gf_compute(MF, U,'L2 norm', meshim mim [,CVLST])
%   Compute the L2 norm of the (real or complex) field U.  If CVLST is indicated,
%   the norm will be computed only on the listed convexes.
% 
%   * N=gf_compute(MF, U,'H1 semi norm', meshim mim [,CVLST])
%   Compute the L2 norm of grad(U).  If CVLST is given, the norm will be computed
%   only on the listed convexes.
% 
%   * N = gf_compute(MF, U,'H1 norm', meshim mim [,CVLST])
%   Compute the H1 norm of U.  If CVLST is given, the norm will be computed only
%   on the listed convexes.
% 
%   * N = gf_compute(MF, U,'H2 semi norm', meshim mim [,CVLST])
%   Compute the L2 norm of D^2(U).  If CVLST is given, the norm will be computed
%   only on the listed convexes.
% 
%   * N = gf_compute(MF, U,'H2 norm', meshim mim [,CVLST])
%   Compute the H2 norm of U.  If CVLST is given, the norm will be computed only
%   on the listed convexes.
% 
%   * DU = gf_compute(MF, U,'gradient', meshfem MFGRAD)
%   Compute the gradient of the field U defined on meshfem MFGRAD.  The gradient
%   is interpolated on the mesh_fem MFGRAD, and returned in DU.  For example, if U
%   is defined on a P2 mesh_fem, DU should be evaluated on a P1-discontinuous
%   mesh_fem. MF and MFGRAD should share the same mesh.  U may have any number of
%   dimensions (i.e. this function is not restricted to the gradient of scalar
%   fields, but may also be used for tensor fields). However the last dimension of
%   U has to be equal to the number of dof of MF. For example, if U is a
%   3x3xnbdof(MF) array, DU will be a Nx3x3[xQ]xnbdof(MFGRAD) array, where N is
%   the dimension of the mesh, and the optional Q dimension is inserted if
%   qdim(MF) != qdim(MFGRAD).
% 
%   * DU = gf_compute(MF, U,'hessian', meshfem MFHESS)
%   Compute the hessian of the field U defined on meshfem MFHESS.
% 
%   * U2 = gf_compute(MF, U,'interpolate_on', { meshfem MF2 | slice SL })
%   Interpolate a field on another mesh_fem or a slice.  * interpolation on
%   another mesh_fem MF2: MF2 has to be Lagrangian.  If MF and MF2 share the same
%   mesh object, the  interpolation will be much faster.
%   
%     * interpolation on a slice object: this is similar to interpolation on a
%   refined P1-discontinuous mesh, but it is much faster.  This can also be used
%   with gf_slice('points') to obtain field values at a given set of points.
%   
%   See also gf_asm('interpolation_matrix')
% 
%   * U2 = gf_compute(MF, U,'extrapolate_on', meshfem MF2)
%   Extrapolate a field on another mesh_fem.  If the mesh of MF2 is stricly
%   included in the mesh of MF, this function does stricly the same job as
%   gf_compute(MF, U,'interpolate_on'). However, if the mesh of MF2 is not exactly
%   included in MF (imagine interpolation between a curved refined mesh and a
%   coarse mesh), then values which are slightly outside MF will be extrapolated.
%   
%   See also gf_asm('extrapolation_matrix')
% 
%   * E=gf_compute(MF, U,'error_estimate', meshim MIM)
%   Compute an a posteriori error estimation.  Currently there is only one which
%   is available: for each convex, the jump of the normal derivative is integrated
%   on its faces.
% 
% 
%   * [U2[,MF2,[,X[,Y[,Z]]]]] = gf_compute(MF,U,'interpolate on Q1 grid', 
%                                {'regular h', hxyz | 'regular N',Nxyz |
%            			   X[,Y[,Z]]}
% 
%   Creates a cartesian Q1 mesh fem and interpolates U on it. The
%   returned field U2 is organized in a matrix such that in can be drawn
%   via the MATLAB command 'pcolor'. The first dimension is the Qdim of
%   MF (i.e.  1 if U is a scalar field)
% 
%   example (mf_u is a 2D mesh_fem): 
%    >> Uq=gf_compute(mf_u, U, 'interpolate on Q1 grid', 'regular h', [.05, .05]);
%    >> pcolor(squeeze(Uq(1,:,:)));
% 
%   * E = gf_compute(MF, U, 'mesh edges deformation', N [,vec or 
%                    mat CVLIST])
%   [OBSOLETE FUNCTION! will be removed in a future release]
%   Evaluates the deformation of the mesh caused by the field U (for a
%   2D mesh, U must be a [2 x nb_dof] matrix). N is the refinment level
%   (N>=2) of the edges.  CVLIST can be used to restrict the computation
%   to the edges of the listed convexes ( if CVLIST is a row vector ),
%   or to restrict the computations to certain faces of certain convexes
%   when CVLIST is a two-rows matrix, the first row containing convex
%   numbers and the second face numbers.
% 
%   * UP = gf_compute(MF, U, 'eval on triangulated surface', int Nrefine,
%                     [vec CVLIST])
%   [OBSOLETE FUNCTION! will be removed in a future release]
%   Utility function designed for 2D triangular meshes : returns a list
%   of triangles coordinates with interpolated U values. This can be
%   used for the accurate visualization of data defined on a
%   discontinous high order element. On output, the six first rows of UP
%   contains the triangle coordinates, and the others rows contain the
%   interpolated values of U (one for each triangle vertex) CVLIST may
%   indicate the list of convex number that should be consider, if not
%   used then all the mesh convexes will be used. U should be a row
%   vector.
% 
%   $Id: gf_compute.cc,v 1.7 2006/03/28 10:06:35 pommier Exp $
function [varargout]=gf_compute(varargin)
  if (nargin>=3 & strcmpi(varargin{3},'interpolate on Q1 grid')),
    [varargout{1:nargout}]=gf_compute_Q1grid_interp(varargin{[1 2 4:nargin]}); return;
  end;
  if (nargout),
    [varargout{1:nargout}]=gf_matlab('compute',varargin{:});
  else
    gf_matlab('compute',varargin{:});
    if (exist('ans','var') == 1), varargout{1}=ans; end;
  end;
