function Fe = FEMMatrices(Fe);
%   FE = FEMMatrices(FE);
% requires a structure FE that is generated by readFEMfile().
% This program adds the various projection matrices to the structure
% Most of the work is done by FEMgradoperator(elem,X), but this
% subroutine takes of course all the credit ...
% ( ... even though all it does is making some transposition operations :-).
% Author ReWT (Aug 15, 2001)
  
  
  [nelem,nlocalnodes]=size(Fe.elem);
  [nnodes,ndim]=size(Fe.xyz);

  fprintf('FEMgradoperator ...\n');

  [ng,nshape,Fe.Nk,Fe.Ax,Fe.Ay,Fe.Az,Fe.Ass,J0,W]  = FEMgradoperator(Fe.elem,Fe.xyz);

  Fe.nelem=nelem;
  Fe.ngauss=ng*nelem;
  Fe.nnodes=nnodes;  
  
  fprintf('defining Jacobians and related things...\n');
  Fe.W = spdiags(W(:).*J0(:),[0],Fe.ngauss,Fe.ngauss);
  Fe.Wv = W(:).*J0(:);
  Fe.J0 = spdiags(J0(:),[0],Fe.ngauss,Fe.ngauss);
  
  Fe.Nk1 = Fe.Ax*Fe.Ass';
  Fe.Nk2 = Fe.Ay*Fe.Ass';
  Fe.Nk3 = Fe.Az*Fe.Ass';
  
  fprintf('define Mass matrix \n');
  Fe.Mass = Fe.Ass * Fe.Nk' * Fe.W * Fe.Nk * Fe.Ass';

  fprintf('calculate inverse Mass matrix ... skipped \n');
  % Fe.IMass = inv(Fe.Mass);
  
return;