for loo = 1:Nloop,
  fprintf('Step %4i  of %4i\n',loo,Nloop);
  
  Fe = FEMFields(Fe);
% given U compute deformation gradient F, Cauchy C, 
% and other related fields (F^(-1), det(F), etc.) 
  
  figure(1);

  Fe = NeoHookean(Fe,activation); 
  Fe = FEMDbar(Fe);
  Fe = Stiffnessmatrix(Fe);
  
% Overwrite
%  Fe.S = zeros(Fe.ngauss,6);
%  Fe.S(:,6) = Fe.S(:,6)-100;
  

%  Fe.S = zeros(Fe.ngauss,6);
%  Fe.S(:,1) = 100*ones(Fe.ngauss,1);
  
  
%  Fe = Deviatoric(Fe);
  
  Fe = FEMForces(Fe,2);
  scal=0.001;

  deltaU = Fe.deltaU;
%  showforces(Fe,0.1,'nodes');
  
%  fprintf('min lfib: %d  /max:  %d / max tau: %d\n',min(Fe.lfib),max(Fe.lfib), ...
%	  max(Fe.muscltau));
  
  Konstructor('deform', deltaU);
  disp(Fe.U(8,:));
  
  Fe.U = Fe.U + deltaU;
  disp(Fe.U(8,:));
end;