% FUNCTION gf_asm(operation[, arg])
% 
%   General assembly function.
% 
%   Many of the functions below use more than one mesh_fem: the main
%   mesh_fem (mf_u) used for the main unknow, and data mesh_fem (mf_d)
%   used for the data. It is always assumed that the Qdim of mf_d is
%   equal to 1: if mf_d is used to describe vector or tensor data, you
%   just have to "stack" (in fortran ordering) as many scalar fields as
%   necessary.
%   
%   * [vec V]=gf_asm('volumic_source', meshim mim, meshfem mf_u, meshfem mf_d,
%   @dcvec Fd)
%   Assembly of a volumic source term.  Output a vector V, assembled on the
%   mesh_fem MF_U, using the data vector FD defined on the data mesh_fem MF_D. FD
%   may be real or complex-valued.
% 
%   * gf_asm('boundary_source', meshim mim, boundary_num, meshfem mf_u, meshfem
%   mf_d, vec G)
%   Assembly of a boundary source term.  G should be a [QDIM x N] matrix, where N
%   is the number of degree of freedom of mf_d, and
%   QDIM=MESHFEM::GET{mf_u}('qdim').
% 
%   * gf_asm('mass_matrix', mim, mf1[, mf2])
%   Assembly of a mass matrix.
% 
%   * gf_asm('laplacian', mim, meshfem mf_u, meshfem mf_d, vec A)
%   Assembly of the Laplacian ( div(a(x)*grad U) with a scalar).
% 
%   * gf_asm('linear_elasticity', meshim mim, meshfem mf_u, meshfem mf_d, vec
%   lambda_d, vec mu_d)
%   Assembles of linear elasticity.  Lambda and mu are the usual Lamé
%   coefficients.
% 
%   * []=gf_asm('nonlinear_elasticity',meshim mim, meshfem mf_u, vec U, string law,
%   meshfem mf_d, mat params, {'tangent matrix' | 'rhs' | 'incompressible tangent
%   matrix',mf_p,P | 'incompressible rhs', mf_p, P})
%   Assembles terms (tangent matrix and right hand side) for nonlinear elasticity.
%   The solution U is required at the current time-step. The law may be choosen
%   among:    - 'SaintVenant Kirchhoff' (linearized law, should be   avoided).
%   This law has the two usual Lamé coefficients as   parameters: params =
%   [lambda; mu]    - 'Mooney Rivlin' (only for incompressibility). This law has
%   two parameters, called C1 and C2.    - 'Ciarlet Geymonat'. This law has 3
%   parameters : params =   [lambda; mu; gamma'(1)], with gamma'(1) chosen such
%   that   gamma'(1) is in ]-lambda/2-mu, -mu[    The parameters of the material
%   law are described of the   mesh_fem mf_d. The matrix 'params' should have
%   nbdof(mf_d)   columns, each row correspounds to a parameter.    The last
%   argument selects what is to be built: either the   tangent matrix, or the
%   right hand side. If the   incompressibility is considered, it should be
%   followed by a   mesh_fem 'mf_p', for the pression.
% 
%   * [K,B]=gf_asm('stokes', meshim mim, meshfem mf_u, meshfem mf_p, meshfem
%   mf_data, vec visc)
%   Assembly of matrices for the Stokes problem { visc(x).Delta(U) - p = 0, div(U)
%   = 0 }.  On output, K is the usual linear elasticity stiffness matrix with
%   lambda=0 and 2*mu=visc. B is a sparse matrix corresponding to $\int p.div v$.
% 
%   * [A]=gf_asm('helmholtz', meshim mim, meshfem mf_u, meshfem mf_data, @cvec k2)
%   Assembly of the matrix for the Helmholtz problem { Laplacian(u) + k^2u = 0 }
% 
%   * [A]=gf_asm('bilaplacian', meshim mim, meshfem mf_u, meshfem mf_data, vec a)
%   Assembly of the matrix for the Bilaplacian problem {
%   Laplacian(a(x)*Laplacian(u)) }
% 
%   * [HH,RR]=gf_asm('dirichlet', int boundary_num, meshim mim, meshfem mf_u,
%   meshfem mf_d, mat H, vec R [, threshold])
%   Assembly of Dirichlet conditions of type 'h.u = r'.  Handle 'h.u = r' where h
%   is a square matrix (of any rank) whose size is equal to
%   MESHFEM:GET{mf_u}('qdim'). This matrix is stored in 'H', one column per dof in
%   'mf_d', each column containing the values of the matrix 'h' stored in fortran
%   order: H(:,j) = [h11(x_j) h21(x_j) h12(x_j) h22(x_j)] if u is a 2D vector
%   field.
%   
%     Of course, if the unknown is a scalar field, you just have to set H =
%   ones(1, MESFEM:GET{mf_d}('nbdof')).
%   
%     This is basically the same than calling gf_asm('boundary_qu_term') for H and
%   calling gf_asm('neumann') for R, except that this function tries to produce a
%   'better' (more diagonal) constraints matrix (when possible).  see also
%   gf_spmat_get(M,'Dirichlet_nullspace').
% 
%   * M=gf_asm('boundary_qu_term', int boundary_num, meshim mim, meshfem mf_u,
%   meshfem mf_d, mat q)
%    Q should be be a [Qdim x Qdim x N] array, where N is the number of degree of
%   freedom of mf_d, and QDIM is dimension of the unkown u (that is set when
%   creating the mesh_fem).
% 
%   * [Q,G,H,R,F]=gf_asm('pdetool_boundary_conditions',mf_u, mf_d, b, e[, f_expr])
%   Assembly of pdetool boundary conditions.   B is the boundary matrix exported
%   by pdetool, and E is the edges array. 'F_EXPR' is an optionnal expression (or
%   vector) for the volumic term. On return Q,G,H,R,F contain the assembled
%   boundary conditions (Q and H are matrices), similar to the ones returned by
%   the function ASSEMB from PDETOOL.
% 
%   * [...]=gf_asm('volumic'[,CVLST], expr [[, mesh_ims, mesh_fems, data...]])
%   Generic assembly procedure for volumic assembly.   The expression 'expr' is
%   evaluated over the mesh_fems listed in the arguments (with optional data) and
%   assigned to the output arguments. For details about the syntax of assembly
%   expressions, please refer to the getfem user manual (or look at the file
%   getfem_assembling.h in the getfem++ sources).
%   
%     For example, the L2 norm of a field can be computed with
%   
%    gf_compute(MF, U,'L2 norm') or with:
%   
%     gf_asm('volumic','u=data(#1);
%   V()+=u(i).u(j).comp(Base(#1).Base(#1))(i,j)',mim,mf,U)
%   
%     The Laplacian stiffness matrix can be evaluated with
%   
%    gf_asm('Laplacian',mim, mf, A) or equivalently with:
%   
%    gf_asm('volumic','a=data(#2);M(#1,#1)+=sym(comp(Grad(#1).Grad(#1).Base(#2))(:
%   ,i,:,i,j).a(j))', mim,mf, A);
% 
%   * [...]=gf_asm('boundary', @tint bnum, string expr [[, mesh_ims, mesh_fems,
%   data...]])
%   Generic boundary assembly.  See the help for gf_asm('volumic').
% 
%   * M=gf_asm('interpolation_matrix', MF1, MF2)
%   Build the interpolation matrix from a meshfem onto another meshfem.  Return a
%   sparse matrix M, such that V=M*U is equal to
%   ::COMPUTE{MF1,U}('interpolate_on',MF2). Useful for repeated interpolations.
%   Note that this is just interpolation, no elementary integrations are involved
%   here, and MF2 has to be lagrangian. In the more general case, you would have
%   to do a L2 projection via the mass matrix.
% 
%   * M=gf_asm('extrapolation_matrix', MF1, MF2)
%   Build the extrapolation matrix from a meshfem onto another meshfem.  Return a
%   sparse matrix M, such that V=M*U is equal to
%   ::COMPUTE{MF1,U}('extrapolate_on',MF2). Useful for repeated extrapolations.
% 
function [varargout]=gf_asm(varargin)
  if (nargin>=1 & strcmpi(varargin{1},'pdetool boundary conditions')),
    [varargout{1:nargout}]=gf_asm_pdetoolbc(varargin{[1 3:nargin]}); return;
  end;
  if (nargout),
    [varargout{1:nargout}]=gf_matlab('asm',varargin{:});
  else
    gf_matlab('asm',varargin{:});
    if (exist('ans','var') == 1), varargout{1}=ans; end;
  end;
