function Fe=FEMForces(Fe,how)
% Fe=FEMForces(Fe,how);
% if how is 1, or missing, the first Piola tensor is used
% directly. If it is 2, the first Piola tensor is computed
% by the operation P=F*S, where F is the deformation tensor,
% S is the 2nd Piola tensor and P is the first Piola tensor. 
%
% Piola: 
%   P11 P12 P13          1  4  7
%   P21 P22 P23          2  5  8
%   P31 P32 P33          3  6  9
%    

  if nargin>2,
    how=1;
  end;
  

    Fe.Piola(:,1)= Fe.F(:,1).*Fe.S(:,1)+Fe.F(:,4).*Fe.S(:,2)+Fe.F(:,7).*Fe.S(:,3);
    Fe.Piola(:,2)= Fe.F(:,2).*Fe.S(:,1)+Fe.F(:,5).*Fe.S(:,2)+Fe.F(:,8).*Fe.S(:,3);
    Fe.Piola(:,3)= Fe.F(:,3).*Fe.S(:,1)+Fe.F(:,6).*Fe.S(:,2)+Fe.F(:,9).*Fe.S(:,3);
    
    Fe.Piola(:,4)= Fe.F(:,1).*Fe.S(:,2)+Fe.F(:,4).*Fe.S(:,4)+Fe.F(:,7).*Fe.S(:,5);
    Fe.Piola(:,5)= Fe.F(:,2).*Fe.S(:,2)+Fe.F(:,5).*Fe.S(:,4)+Fe.F(:,8).*Fe.S(:,5);
    Fe.Piola(:,6)= Fe.F(:,3).*Fe.S(:,2)+Fe.F(:,6).*Fe.S(:,4)+Fe.F(:,9).*Fe.S(:,5);
    
    Fe.Piola(:,7)= Fe.F(:,1).*Fe.S(:,3)+Fe.F(:,4).*Fe.S(:,5)+Fe.F(:,7).*Fe.S(:,6);
    Fe.Piola(:,8)= Fe.F(:,2).*Fe.S(:,3)+Fe.F(:,5).*Fe.S(:,5)+Fe.F(:,8).*Fe.S(:,6);
    Fe.Piola(:,9)= Fe.F(:,3).*Fe.S(:,3)+Fe.F(:,6).*Fe.S(:,5)+Fe.F(:,9).*Fe.S(:,6);

  % in the operation P_ij N,j,:
  % N,1 multiplies with the first column only and is stored in F1  
  % N,2 with the second ..

%  Fe.Piola = Fe.W* Fe.Piola;
  
  F1 =  Fe.Ax' * Fe.W * Fe.Piola(:,[1 2 3]);
  F2 =  Fe.Ay' * Fe.W * Fe.Piola(:,[4 5 6]);
  F3 =  Fe.Az' * Fe.W * Fe.Piola(:,[7 8 9]);
  Fe.F1 = F1;
  Fe.F2 = F2;
  Fe.F3 = F3;
  Fe.forcy = F1+F2+F3;

  % new with stiffness matrix
  forc = Fe.Ass * Fe.forcy;
  Fe.Force = forc;
  
  % force is actually the 
  
  was=2;
  if was == 1,
    forc = forc';
    forc = forc(:);
    
    Fe.fff = forc(:);
  
    fffI = zeros(Fe.nnodes*3,1);
    fffI(Fe.ncond) = Fe.Stiff \ Fe.fff(Fe.ncond);
    Fe.Force=(reshape(fffI,3,Fe.nnodes))';
  else
    forc=forc(:);
    Fe.fff = forc(:);
    fffI = zeros(Fe.nnodes*3,1);
    fffI(Fe.ncond) = Fe.Stiff\Fe.fff(Fe.ncond);
    Fe.deltaU = reshape(fffI,Fe.nnodes,3);
    disp(Fe.deltaU);
  end;
  
  return;