CheckMesh - Check that array sizes match for mesh structure.
   
   USAGE:

   CheckMesh(mesh)
   CheckMesh(mesh, vector)

   INPUT:

   mesh   is a MeshStructure
   vector is a vector, (optional) 
          a reduced vector of values on the mesh nodal points

   OUTPUT: calls `error' if inconsistency is found

   NOTES:

   This routine checks that the coordinates and equivalence
   array of a mesh are consistent.  If the vector argument
   is given, it checks that it is of the correct length.
   If any array sizes are inconsistent, it prints an error
   message and returns.

 DiscreteDelta - Evaluate sum of discrete delta functions.
   
   USAGE:

   fun = DiscreteDelta(mesh, l2ip, elem, ecrd, wts)

   INPUT:

   mesh is a MeshStructure
   l2ip is n x n, 
        the inner product matrix for the mesh, where
        n is the number of independent degrees of freedom
   elem is an n-vector of integers, 
        the list of elements containing the delta-function points 
   ecrd is m x n, 
        the list of barycentric coordinates of the delta-function points
   wts  is an n-vector,(optional)
        it gives the weights of each point; if not present,
        the weights are set to one

   OUTPUT:

   fun is an n-vector, 
       the nodal point values of the sum of the discrete
       delta functions associated with the points and weights

   NOTES:  

   *  The resulting function is not normalized in any way.

 EqvReduce - Reduce matrix columnwise by equivalences.
   
   USAGE:

   rmat = EqvReduce(mat, eqv)

   INPUT:

   mat is m x n, (usually sparse) a matrix, the columns of which
                 have underlying nodal point equivalences
   eqv is 2 x k, the equivalence array

   OUTPUT:

   rmat is m x (n-k), the new matrix for the reduced (condensed) 
                 set of nodes formed by adding equivalent columns 
                 to the master column

   NOTES:

   *  The columns of the unreduced nodes are added to the 
      columns of the master node.

   *  This routine only reduces along the column dimension.  To
      reduce along the rows, apply to the transpose.  To reduce
      along both dimensions, call it twice.

   *  The equivalence array can be empty, in which case nothing
      is done, and the original matrix is returned.

 EvalMeshFunc - Evaluate mesh function.
   
   USAGE:

   values = EvalMeshFunc(mesh, fun, els, ecrds)

   INPUT:

   mesh  is a MeshStructure
   fun   is an n-vector, 
         function values on the set of independent nodal points
   els   is an m-vector of integers, 
         it gives the list of elements containing the points at
         which to evaluate the function
   ecrds is k x m, 
         it is the list of barycentric coordinates of the 
         points at which to evaluate the function

   OUTPUT:

   values is an n-vector,
          the values of the function at the specified points

 H1SIPMatrix -- H^1 semi-inner product [Documentation out of date!]

   USAGE:

   h1sip = H1SIPMatrix(mesh, qrule, diff)

   INPUT:

   diff is d1 x d2 x n
        the list of tangent vectors at each of n points
   npqp is m x n,
        the matrix which takes nodal point values to
        (global) quadrature points 

   OUTPUT:

   h1sip is n x n, sparse
       it is the matrix of the H^1 semi-inner product,
	where n is the number of independent degrees of
	freedom associated with the mesh
        

 Jacobian - Compute Jacobian of linear mesh mappings.

   USAGE:

   jac = Jacobian(mesh)

   INPUT:

   mesh is a MeshStructure,
        with simplicial element type

   OUTPUT:

   jac is 1 x m, 
       the Jacobian of each element

   NOTES:

   *  The mesh may be embedded in a space of higher 
      dimension than the reference element.  In that
      case, the Jacobian is computed as (sqrt(det(J'*J))
      and is always positive.  When the target space is
      of the same dimension as the reference element,
      the Jacobian is computed as usual and can be
      positive or negative.

   *  Only simplicial (linear) element types are allowed.

 L2IPMatrix - Form matrix for L2 inner product.
   
   USAGE:

   l2ip = L2IPMatrix(gqrule, npqp)

   INPUT:

   qrule is a QRuleStructure,
         it is the global quadrature rule associated with
         the underlying mesh
   npqp  is m x n, (sparse) 
         it is the matrix which takes nodal point values
         to quadrature point values

   OUTPUT:

   l2ip is n x n, (sparse) 
        the inner product matrix associated with the mesh

 LoadMesh - Load mesh from a ascii file.
  
   USAGE:

   mesh = LoadMesh(name)

   INPUT:

   name is a string, 
       the basename of the mesh files

   OUTPUT:

   mesh is a MeshStructure, 
        for the mesh being loaded

   NOTES:

   * Expected file suffixes are:

        .crd for nodal point coordinates,
        .con for connectivity
        .eqv for equivalences (optional)
   

 LoadQuadrature - Load quadrature data rule.
   
   USAGE:

   qrule = LoadQuadrature(qname)

   INPUT:

   qname is a string, 
         the basename of the quadrature data files

   OUTPUT:

   qrule is a QRuleStructure, 
         it consists of the quadrature point locations and weights

   NOTES:

   *  It is expected that the quadrature rules are for simplices,
      and the last barycentric coordinate is appended to the file
      read from the data.

   *  Expected suffixes are .qps for the location and .qpw for
      the weights.

 MeanValue - Find integral mean value of function.
   
   USAGE:

   mv = MeanValue(f, l2ip)

   INPUT:

   f    is an n-vector, 
        an array of nodal point function values 
   l2ip is n x n, 
        the (L2) inner product matrix for the underlying mesh

   OUTPUT:

   mv is 1 x 1, 
      the mean value of `f'---the integral of f divided by
      the measure of the region

 MeshFaces - Find lower dimensional faces for a mesh.
   
   USAGE:

   [faces, (opt) mult] = MeshFaces(con)

   INPUT:

   con is d x n, 
       the mesh connectivity (tetrahedra or triangles)

   OUTPUT:

   faces is (d-1) x m, (integer)
         the connectivity of the lower dimensional
                       faces of con
   mult is 1 x m,  (integer, optional)
        the multiplicity of each face; this should
        be either 1 (for a surface face) or 2 (for
        an interior face)

 MeshInfo - Verify and print properties of mesh.
   
   USAGE:

   outinfo = MeshInfo(mesh)
   outinfo = MeshInfo(mesh, sym)
   outinfo = MeshInfo(mesh, sym, tol)

   INPUT:

   mesh is a MeshStructure
   sym  is 4 x n,  (optional)
        the symmetry group in quaternions; this argument can 
        be omitted if the equivalence array is empty or if you 
        do not desire to verify it
   tol  is a scalar, (optional, default = 1.0e-7)
        the tolerance used in verifying and generating 
        equivalences 

   OUTPUT:  none

   This function prints various information to the screen.

 MeshStructure - Create mesh structure from mesh data.
   
   USAGE:

   mesh = MeshStructure
   mesh = MeshStructure(crd, con)
   mesh = MeshStructure(crd, con, eqv)

   INPUT:

   crd is e x n, 
       the array of nodal point locations
   con is d x m, (integer)
       the mesh connectivity
   eqv is 2 x k, (integer, optional) 
       the equivalence array

   OUTPUT:

   mesh is a MeshStructure, 
        the basic MeshStructure consists of three fields,
        the nodal point coordinates (.crd), the connectivity
        (.con) and the nodal point equivalence array (.eqv).
        
   NOTES:
 
   *  With no arguments, this function returns and empty mesh
      structure.  With only two arguments, it sets the equivalence 
      array to be empty.

 NpQpMatrix - Nodal point values to quadrature points.
   
   USAGE:

   npqp = NpQpMatrix(mesh, qrule)

   INPUT:

   mesh  is a MeshStructure
   qrule is a QRuleStructure, 
         a quadrature rule on the reference element

   OUTPUT:

   npqp is m x n, (sparse) 
        it is the matrix which takes the values
        at the independent nodes of the mesh to the values 
        at the quadrature points;  here, n is the number of 
        independent nodal values, and m = q*e, where e is the 
        number of elements in the mesh and q is the number of 
        qudrature points per element
                  

 QRuleGlobal - Construct a global quadrature rule for a mesh.
   
   USAGE:

   gqrule = QRuleGlobal(mesh, qrule)
   gqrule = QRuleGlobal(mesh, qrule, @MetricFun)

   INPUT:

   mesh  is a MeshStructure
   qrule is a QRuleStructure, 
         the quadrature rule for the reference element 
  
   MetricFun is a function handle, (optional)
             It returns the volumetric integration factors
             for the applied mapping.  It is for the case
             in which the mesh is the parameter space for
             some manifold and has a subsequent mapping
             associated with each point.
            
             met = Metricfun(pts)

             pts is m x n array of points
             met is 1 x n vector of integration factors

   OUTPUT:

   NOTES:

   *  For spheres, this routine is superceded by `SphGQRule',
      which handles the metric dependency on the mesh.  This
      still is to be used for meshes on Rodrigues parameters,
      for which the metric function is independent of the mesh,
      or for meshes with no further mapping involved.
             

 QRuleStructure - Quadrature rule structure.
   
   USAGE:

   qrule = QRuleStructure
   qrule = QRuleStructure(pts, wts)

   INPUT:

   pts is m x n, 
       a list of n m-dimensional points 
   wts is an n-vector, 
       the associated weights

   OUTPUT:

   qrule is a QRuleStructure, 
         it consists of two fields, points (.pts) and
         weights (.wts)

   NOTES:

   *  With no arguments, this returns an empty structure.

 ReduceMesh - Find equivalence array for a Rodrigues mesh.
   
   USAGE:

   newmesh = ReduceMesh(mesh, sym)
   newmesh = ReduceMesh(mesh, sym, tol)

   INPUT:

   mesh is a MeshStructure on the fundamental region
        with (possibly) an empty equivalence array
   sym  is 4 x m, 
        the symmetry group in quaternions
   tol  is 1 x 1, (optional, default: 10^-14) 
        the tolerance used for comparing equivalent rotations; 

   OUTPUT:

   newmesh is a MeshStructure on the fundamental region;
           it has the same  and connectivity, but
           reordered; the equivalence array is created
           and the nodes are numbered so that all 
           the independent nodes precede the dependent nodes.

 RefDerivatives - Reference shape function derivatives.
   
   USAGE:

   der = RefDerivatives(n)

   INPUT:

   n is a positive integer,
     the dimension of the simplex

   OUTPUT:

   der is n x (n+1),
	column j contains the gradient of the j'th
       barycentric coordinate with respect to the 
       coordinate directions

 RefineMesh - Refine a simplex-based mesh.
   
   USAGE:

   rmesh = RefineMesh(mesh, n)
   rmesh = RefineMesh(mesh, n, sym)

   INPUT:

   mesh is a MeshStructure,
        using 1D, 2D, or 3D simplicial elements
   n    is a scalar, (positive integer)
        the subdivision parameter for each simplex
   sym  is 4 x s, (optional) 
        the array of quaternions giving the symmetry
        group for meshes on orientation space

   OUTPUT:

   rmesh is a MeshStructure,
         it is derived from the input mesh by subdividing each 
         simplicial element into a number of subelements 
         based on the subdivision parameter n
  
   NOTES:

   *  This routine only subdivides the simplex.  If the 
      mesh points are mapped, as in sphere meshes, then
      that mapping needs to be applied after this routine.
   

 SliceMesh - Planar slice through a 3D mesh
   
   USAGE:

   sm = SliceMesh(m, p, n)

   INPUT:

   m is a MeshStructure,
     for a general mesh
   p is a 3-vector,
     a point on the slice plane
   n is a 3-vector,
     the normal to the slice plane; it does not
     need to be unit length

   OUTPUT:

   sm is a MeshStructure,
      on the 2D planar section
   smels is an integer array,
      the list of elements containing the nodes in `sm',
      relative to the original mesh, `m'
   smecrd is 4 x n,
      the barycentric (elemental) coordinates of the nodes
      in `sm' relative to the original mesh, `m'

 SmoothFunction - Smooth by convolution.
   
   USAGE:

   fsm = SmoothFunction(f, mesh, gqrule, @SmoothFun)
   fsm = SmoothFunction(f, mesh, gqrule, @SmoothFun, arg1, ...)

   INPUT:

   f      is an n-vector, 
          the values of a function on the mesh at the
          set of independent nodal points
   mesh   is a MeshStructure,
          for any region
   qpts   is d x n,
          the barycentric coordinates of the quadrature points
   gqrule is a QRuleStructure,
          the global quadrature rule for the mesh

   SmoothFun is a function handle,
             which defines the function to use for
             smoothing; it has the form:

             SmoothFun(center, points, ...)

             center is m x 1, 
                    the center of the distribution
             points is m x n, 
                    the list of points to evaluate

   Other arguments are passed to `SmoothFun'

   OUTPUT:

   fsm is 1 x n, 
      the nodal point values of the smoothed function at
      the set of independent nodal points

   NOTES:
 
   *  The function is smoothed by convolving a fixed distribution
      (e.g. Gaussian) with the given function.  The result at 
      each nodal point value is the integral of the input function 
      times the smoothing function centered at that nodal point.

   *  The result is not normalized.  
      

 SpreadRefPts - Spread reference coordinates to all elements.
   
   USAGE:

   pts = SpreadRefPts(mesh, ref)

   INPUT:

   mesh is a MeshStructure 
        with simplicial element type
   ref is e x k, 
       a list of barycentric coordinates for points in
       the reference element; 

   OUTPUT:

   pts is d x k x m, 
       the spatial coordinates of the reference points under
       the isoparametric mapping defined by the mesh;
       `d' is the dimension of the mapped space, 'k' is the
       number of reference points, and 'm' is the number
       of elements in `mesh'

   NOTES:

   *  Typically this would be used for numerical quadrature.
   *  e = d + 1 is not necessary; cross-dimension mappings are ok.

 SubdivideSimplex - Regular subdivision of a simplex.

   USAGE:

   ref = SubdivideSimplex(dim, n)
   [ref, con] = SubdivideSimplex(dim, n)

   INPUT:

   dim is a positive integer, 
       the dimension of the simplex
   n   is a positive integer, 
       the number of subdivisions in each direction

   OUTPUT:

   ref is d x m, 
       the list of barycentric coordinates of the points 
       in the subdivision
   con is d x e, (integer)
       the connectivity of the subdivision

 ToAllNodes - Spread array at independent nodes to all nodes.
   
   USAGE:

   all = ToAllNodes(red, eqv)

   INPUT:

   red is m x n, 
       a list of n m-vectors at the independent nodes of a mesh

   eqv is 2 x l, (integer)
       the list of node equivalences (new #, old #)

   OUTPUT:

   all is m x k, 
       the array at all nodes of the mesh; k = n + l