function Fe=FEMFields(Fe) 
% Fe=FEMFields(Fe) 
% 
% Assumes that Fe.U, the displacement is known and 
% further that all the projection matrices are defined.
% compute the following fields:
%   J      
%   Jtild     (J^(-2/3))
%   F      
%   FinvT  
%   C      
%   Cinv   
  
% Computation of F, the deformation tensor.
  
  x_t = Fe.xyz+Fe.U ;
  
  F = [Fe.Nk1*x_t, Fe.Nk2*x_t, Fe.Nk3*x_t];
  
  
  % This is understood as the following computation:
  % F=[F11,F21,F31, F12,F22,F32, F13,F23,F33]
  % F=[F1, F2, F3,  F4, F5, F6,  F7, F8, F9];
  % Or in matrix notation:
  %  [F1,F4,F7
  %   F2,F5,F8
  %   F3,F6,F9]
  
  % Written with the second method, from Mathematica, the determinate is:
  % J =     -F3 F5 F7 + F2 F6 F7 + F3 F4 F8 - F1 F6 F8 - F2 F4 F9 + F1 F5 F9
  
  
  J  =   F(:,2).*F(:,6).*F(:,7) + F(:,3).*F(:,4).*F(:,8) ...
	 - F(:,1).*F(:,6).*F(:,8) - F(:,2).*F(:,4).*F(:,9) ...
	 + F(:,1).*F(:,5).*F(:,9) - F(:,3).*F(:,5).*F(:,7);
  
% The following may not be needed: computation of F^(-T):  
  FinvT(:,1) = ( F(:,5).*F(:,9) - F(:,6).*F(:,8) )./J;
  FinvT(:,2) = ( F(:,6).*F(:,7) - F(:,4).*F(:,9) )./J;
  FinvT(:,3) = ( F(:,4).*F(:,8) - F(:,5).*F(:,7) )./J;
  
  FinvT(:,4) = ( F(:,3).*F(:,8) - F(:,2).*F(:,9) )./J;
  FinvT(:,5) = ( F(:,1).*F(:,9) - F(:,3).*F(:,7) )./J;
  FinvT(:,6) = ( F(:,2).*F(:,7) - F(:,1).*F(:,8) )./J;
  
  FinvT(:,7) = ( F(:,2).*F(:,6) - F(:,3).*F(:,5) )./J;
  FinvT(:,8) = ( F(:,3).*F(:,4) - F(:,1).*F(:,6) )./J;
  FinvT(:,9) = ( F(:,1).*F(:,5) - F(:,2).*F(:,4) )./J;
  
% we need C and Cinv, and they are symmetric.  
  Jtild = J.^(-2/3);

% Cauchy tensor:

% C = [c1 c2 c3]
%     [c2 c4 c5]
%     [c3 c5 c6]
% ditto Cinv.
  C(:,1) = F(:,1).^2 + F(:,2).^2 + F(:,3).^2;
  C(:,2) = F(:,1).*F(:,4) + F(:,2).*F(:,5) + F(:,3).*F(:,6);
  C(:,3) = F(:,1).*F(:,7) + F(:,2).*F(:,8) + F(:,3).*F(:,9);
  C(:,4) = F(:,4).^2 + F(:,5).^2 + F(:,6).^2;
  C(:,5) = F(:,4).*F(:,7) + F(:,5).*F(:,8) + F(:,6).*F(:,9);
  C(:,6) = F(:,7).^2 + F(:,8).^2 + F(:,9).^2;
  
% inverse Cauchy tensor  
  Cinv(:,1) =( -C(:,5).^2 + C(:,4).*C(:,6))./J.^2;
  Cinv(:,2) =( C(:,3).*C(:,5) - C(:,2).*C(:,6))./J.^2; 
  Cinv(:,3) =( -(C(:,3).*C(:,4)) + C(:,2).*C(:,5))./J.^2;
  Cinv(:,4) =( -C(:,3).^2 + C(:,1).*C(:,6))./J.^2;
  Cinv(:,5) =(  C(:,2).*C(:,3) - C(:,1).*C(:,5))./J.^2;
  Cinv(:,6) =(  -C(:,2).^2 + C(:,1).*C(:,4))./J.^2;
  
  
  Fe.J = J;
  Fe.Jtild=Jtild;
  Fe.F = F;
  Fe.FinvT=FinvT;
  Fe.C = C;
  Fe.Cinv=Cinv;  
  si = [1 2 3;
	2 4 5;
	3 5 6];
  for kk=1:36,
    I=Fe.dirbar(kk,1);
    J=Fe.dirbar(kk,2);
    K=Fe.dirbar(kk,3);
    L=Fe.dirbar(kk,4);
    Fe.II(:,kk) = (1/2)*(Fe.Cinv(:,si(I,K)).*Fe.Cinv(:,si(J,L)) + ...
			 Fe.Cinv(:,si(J,K)).*Fe.Cinv(:,si(I,L)));
  end;
    
  return;
  
  
%A_ij B_kl
%A_ji B_kl -
%A_ij B_lk
%A_ji B_lk -
%A_kl B_ij