function [ng,nshape,A0,Ax,Ay,Az,Ass,Jglob,W]  = FEMgradoperator(elem,X)
% [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.\
			
  [nnodes,n3]=size(X);  % total number of nodes.
  [nelem,nperelem]=size(elem); % number of elements.
  ngauss=nelem*ng;      % total number of Gauss points.
                       % assuming all the same type of element. Haha.

  locnodeInx = [1:nshape];      % index for nodes that are used. blink
                                % blink ;-)
    				
% Assembly matrix:				
   el = elem(:,locnodeInx);			
   el = el';
   ii = el(:);
   jj = [1:nelem*nshape]';
   Ass = sparse(ii,jj,ones(nelem*27,1),nnodes,nelem*27);
			
   			
% Ax = sparse(ngauss,nnodes);
% Ay = sparse(ngauss,nnodes);
% Az = sparse(ngauss,nnodes);

% for the definition of 3 sparse matrices, see below:
  iiloc = [1:ng]' * ones(1,nshape);
  jjloc = ones(ng,1) * [1:nshape];
   
  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
  
  for ne=1:nelem,
%    disp(ne);
    x = X(elem(ne,locnodeInx)',1);
    y = X(elem(ne,locnodeInx)',2);
    z = X(elem(ne,locnodeInx)',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;
    
    KK1 = (ne-1)*ng*nshape+1;
    KK2 = ne*ng*nshape;
    k1 = (ne-1)*ng+1;
    k2 = ne*ng;
    ii(KK1:KK2) = iiloc(:)+ng*(ne-1);
    jj(KK1:KK2) = jjloc(:)+nshape*(ne-1);
    
    axfld(KK1:KK2) = grx(:);    % fast index: Gauss point index
    ayfld(KK1:KK2) = gry(:);    % slow index: Shape function 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,nelem*nshape);
   Ax = sparse(ii,jj,axfld,nelem*ng,nelem*nshape);
   Ay = sparse(ii,jj,ayfld,nelem*ng,nelem*nshape);
   Az = sparse(ii,jj,azfld,nelem*ng,nelem*nshape);