function [A0,Ax,Ay,Az,Axt,Ayt,Azt,Ass,Jglob,W]  = FEMgradoperator(elem,X)
  % ATTENTION: This is the old file with the above function name.
% [A0,Ax,Ay,Az,J0]  = FEMgradoperator(elem,X)
% This function returns the three differential 
% operators in for spatial x ,y and z direction in Ax,Ay,Az.
% It further returns the spatial interpolation operator in A0.
% And further the integration constants 
% input: elem is the assembly, and X are the reference
% coordinates of the nodes.
% if specified, it also returns the det(J0) in the 4th field
% The purpose of the matrices Ax,Ay,Az is to obtain the 
% material derivative of any field e defined at the nodes
% The operation Ax*e results in the differential with respect
% to the first material coordinate at all Gausspoints. 
% If x is the current configuration, Ax*x gives the 
% first column of F, the deformation gradient, Ay*x, the second
% and Az*x the third column of F. 
% Author ReWT (Aug 14, 2001) 
  
%  [xp,w]=FEM3x3Gausspts;      % Gausspoint coords and weights in master
%  element
%  [xp,w]=FEM3x3Lobatto;      % Lobatto isn't any better
 
  [xp,w]=FEM4x4Gausspts; 
  DNk1=FEMinterpol(xp,1);     % gradients in master element.
  DNk2=FEMinterpol(xp,2);
  DNk3=FEMinterpol(xp,3);
  NK = FEMinterpol(xp,0);     % shape functions.
  % as originally returned: 
  % 
  %   NK  --k-->  horizontal (second index) corresponds to node number.
  %    |
  %    v 
  %    |
  %    V
  
  [ng,nshape]=size(NK) % ng: num of Gauss pts per element.
                        % nshape num of Shape functions.
			% to change this change FEMinterpol.\
  NK = NK';             % analog to grx,gry,grz, see in loop below.
			
  [nnodes,n3]=size(X)  % total number of nodes.
  [nelem,n27]=size(elem) % number of elements.
  ngauss=nelem*ng      % total number of Gauss points.
                        % assuming all the same type of element. Haha.
  
% Ax = sparse(ngauss,nnodes);
% Ay = sparse(ngauss,nnodes);
% Az = sparse(ngauss,nnodes);

% for the definition of 3 sparse matrices, see below:
  
  ii = zeros(ng*nshape*nelem,1);
  jj = zeros(ng*nshape*nelem,1);

  nkfld = zeros(ng*nshape*nelem,1);
  axfld = zeros(ng*nshape*nelem,1);
  ayfld = zeros(ng*nshape*nelem,1);
  azfld = zeros(ng*nshape*nelem,1);
  Jglob=zeros(ng*nelem,1);      % Jacobian (det) at each Gauss point.
  W = zeros(ng*nelem,1);        % Weights at Gauss points
  onesshape = ones(nshape,1);   % used in the loop below.
  onesng = ones(1,ng);          % ditto
  
  for ne=1:nelem,
%    disp(ne);
    x = X(elem(ne,:)',1);
    y = X(elem(ne,:)',2);
    z = X(elem(ne,:)',3);
% Calculate three gradient fields on the master element:
    a11=DNk1 * x;
    a12=DNk2 * x;
    a13=DNk3 * x;
    
    a21=DNk1 * y;
    a22=DNk2 * y;
    a23=DNk3 * y;
    
    a31=DNk1 * z;
    a32=DNk2 * z; 
    a33=DNk3 * z;
    
% These guys are the lines of the Lamba_ij. 

    J0 = a12.*a23.*a31 + a13.*a21.*a32 - a11.*a23.*a32 - ...
	 a12.*a21.*a33 + a11.*a22.*a33 - a13.*a22.*a31;
    
    
    b11 =   ( a22.*a33 - a23.*a32)./J0;
    b12 =   ( a13.*a32 - a12.*a33)./J0; 
    b13 =   ( a12.*a23 - a13.*a22)./J0;
    b21 =   ( a23.*a31 - a21.*a33)./J0;
    b22 =   ( a11.*a33 - a13.*a31)./J0;
    b23 =   ( a13.*a21 - a11.*a23)./J0;
    b31 =   ( a21.*a32 - a22.*a31)./J0;
    b32 =   ( a12.*a31 - a11.*a32)./J0;
    b33 =   ( a11.*a22 - a12.*a21)./J0;
    
    % next tansform basic gradients: by multiplying with inverse
    % transposed of Lambda. 
    
    for k=1:nshape,  % k runs over the interpolation functions
      grx(:,k)  = DNk1(:,k).*b11 + DNk2(:,k).*b21 + DNk3(:,k).*b31;
      gry(:,k)  = DNk1(:,k).*b12 + DNk2(:,k).*b22 + DNk3(:,k).*b32;
      grz(:,k)  = DNk1(:,k).*b13 + DNk2(:,k).*b23 + DNk3(:,k).*b33;
    end;
    grx=grx';
    gry=gry';
    grz=grz';
    
    k1 = (ne-1)*ng+1;
    k2 = ne*ng;
    KK1 = (ne-1)*ng*nshape+1;
    KK2 = ne*ng*nshape;

    ik = onesshape * [k1:k2]; 
    jk = elem(ne,:)'* onesng; 
    iik = ik;                    % new
    jjk = iik';                   % new

    ii(KK1:KK2) = ik(:);
    jj(KK1:KK2) = jk(:);
    iii(KK1:KK2) = iik(:);
    jjj(KK1:KK2) = jjk(:);

    axfld(KK1:KK2) = grx(:);    % fast index: shape function index.
    ayfld(KK1:KK2) = gry(:);    % slow index: Gauss point index.
    azfld(KK1:KK2) = grz(:);
    nkfld(KK1:KK2) = NK(:);
    Jglob(k1:k2) = J0;
    W(k1:k2) = w(:);   % weights
  end;
   
   A0 = sparse(ii,jj,nkfld,nelem*ng,nnodes);
   Ax = sparse(ii,jj,axfld,nelem*ng,nnodes);
   Ay = sparse(ii,jj,ayfld,nelem*ng,nnodes);
   Az = sparse(ii,jj,azfld,nelem*ng,nnodes);
   
   Axt = sparse(iii,jjj,axfld,nelem*ng,nelem*ng);
   Ayt = sparse(iii,jjj,ayfld,nelem*ng,nelem*ng);
   Azt = sparse(iii,jjj,azfld,nelem*ng,nelem*ng);
   
   el = elem';
   ii = el(:);
   jj = [1:nelem*27]';
   Ass = sparse(ii,jj,ones(nelem*27,1),nnodes,nelem*27);