function Fe= FEMtestK(what,Fe,dx)

switch what,
  case 'init',
   close all;  fireup; Konstructor('loadmacro','OneCube.mac');
   initialdeform=1;
   fibequal=1;
   filename = 'OneCube.FEM';
   
   Fe = readFEMfile(filename);
   
   if fibequal==1,
     for k=1:length(Fe.fiberdensity),
       Fe.fiberdensity{k} = ones(size(Fe.fiberdensity{k}));
    end;
   end;
   
   Fe.constr = zeros(size(Fe.constr));
   Fe = FEMMatrices(Fe);
   
   Fe.U = zeros(Fe.nnodes,3);
   
   con = Fe.constr;
   con = (con==0) | (con>2);
   con = ones(size(con));
   
   fu = con';
   fu = fu(:);
   nn = find(fu);
   Fe.ncond = nn;
   
   [Fe.cindxbar,Fe.dirbar,Fe.multibar]=fourthindicies;
   [Fe.cindxlag,Fe.dirlag,Fe.multilag]=lagrangeindicies;
   
   
 case 'continue',
  fibequal=1;
  con = Fe.constr;
  con = (con==0) | (con>2);
  con = ones(size(con));
  
   fu = con';
   fu = fu(:);
   nn = find(fu);
  % deformation in x dir: 
 % dz=1./sqrt(1+dx);
  Fd = [1  0   0.0;
	0  1    0;
	0 0.0   1+dx];
  uu = Fe.xyz;
  uu = (Fd*uu')'-Fe.xyz;
  
%  uu = [ Fe.xyz(:,1) * (1+dx) * 0.1 .* Fe.xyz(:,3),   Fe.xyz(:,2) * (1+dx) * 0.1 ...
%      .* Fe.xyz(:,3), Fe.xyz(:,3)];
%  uu = uu - Fe.xyz;
  Fe.U = uu.*con;
  
  El = 10;
  pois = 0.49;
  Fe.lamda = pois*El/((1+pois)*(1-2*pois));
  Fe.mu = El/(2*(1+pois));
  Fe.kappa = 0;
  
  
  Fe.mu=3.55;
  Fe.lamda = 164;
  
  activation=0;
  
  Fe = FEMFields(Fe);
  Fe = NeoHookean(Fe,activation); 
  Fe = FEMDbar(Fe);
  Fe = Stiffnessmatrix(Fe);
  
  u = uu(:);
  u = u(Fe.ncond);
  
  f = Fe.Stiff * u;
  ff = zeros(27*3,1);
  ff(Fe.ncond) = f;
  Fe.Force = reshape(ff,27,3);
  
  showforces(Fe,0.1,'nodes');
end;