%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mechanics of growth %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_tube1; % vol
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_tube2; % vol
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_tube3; % vol
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_beams; % vol 
 [q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_humr1; % den 
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_humr2; % den
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_femur; % den
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_frame; % den 
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_cylin; % den
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_tubed; % den
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_bimat; % den
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_unity; % den 
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_punch; % den
%[q0,edof,emat,u_pre,F_pre,mat,ndim,nel,node,ndof,nip,nlod] = ex_block; % den
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
dt     = 1;                       % init pseudo time
time   = dt;                      % init time
tol    = 1e-10;                   % tolerance of newton iteration
e_mat  = extr_dof(edof,q0);       % init material coordinates
dof    = q0;                      % init spatial coordinates
dofn   = dof;                     % init vector of global dofs @time n
i_var  = zeros(nel,nip);          % init internal variables
nsteps = 0;                       % init no of time steps
F_inc  = F_pre;                   % init incremental neuman bc's
u_inc  = u_pre;                   % init incremental dirichlet bc's

[n,nedof] = size(edof);
nedof     = nedof-1;
plot_int(e_mat,i_var,nel,ndim,0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for im=1:1000 
 macro = input('macro:','s'); [ir,ic] = size(macro);
 if ic<4; disp('@ least 4 letters needed for input'); else
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% read data from restart file
  if (strcmp(macro(1:4),'rest') == 1)
    disp('... loading data from restart_time and restart_data file...')
    f01=fopen('restart_time');nsteps=fscanf(f01,'%i\n',1);mat=fscanf(f01,'%f\n',7)';
    f02=fopen('restart_data');dof   =fscanf(f02,'%f\n',ndof); dofn=dof; 
    f03=fopen('restart_ivar');i_var =fscanf(f03,'%f\n',nel*nip);
	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% apply load in instep increments of instep*F_pre through step,,instep 
  elseif (strcmp(macro(1:4),'step') == 1);
    [ir,ic] = size(macro);
    if ic==4; inpstep = 1; else; inpstep = str2num(macro(7:ic)); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% loop over all load steps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for is = (nsteps+1):(nsteps+inpstep);
      iter     = 0; 
      residuum = 1;
      disp(['step: ',num2str(is)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% global newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    while residuum > tol
      iter=iter+1;

      I = zeros(nel,nip);             % init internal variables
      R = zeros(ndof,1);              % init global residual
      K = sparse(ndof,ndof);          % init global stema
      e_spa = extr_dof(edof,dof);     % extract dofs of global vector
      e_par = mat(emat(:),:);         % extract material parameters
	  kkv = zeros(nedof*nedof,nel);   % matrix coefficients
	  kki = zeros(nedof*nedof,nel);   % matrix i indexes
	  kkj = zeros(nedof*nedof,nel);   % matrix j indexes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%% loop over all elements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      for ie = 1:nel
        if ndim==2 
	    [Ke,Re,Ie] = quads_2d(e_mat(ie,:),e_spa(ie,:),i_var(ie,:),e_par(ie,:));
        else
		if nedof==12  
		[Ke,Re,Ie] = tetra_3d(e_mat(ie,:),e_spa(ie,:),i_var(ie,:),e_par(ie,:));  
		else
	    [Ke,Re,Ie] = brick_3d(e_mat(ie,:),e_spa(ie,:),i_var(ie,:),e_par(ie,:));
		end
		end
%%%%%%%% assembly with sparse solver %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           
%     	[K, R, I ] = assm_sys(edof(ie,:),K,Ke,R,Re,I,Ie);
		index       = edof(ie,2:end);
		rows		= index(1:nedof)'*ones(1,nedof);
		cols		= ones(nedof,1)*index(1:nedof);
		kkv(:,ie) = Ke(:);
		kki(:,ie) = rows(:);
		kkj(:,ie) = cols(:);
		R(edof(ie,2:end)) = R(edof(ie,2:end)) + Re;
		I(ie,:)=Ie(:);
	  end
	  K = sparse(kki(:),kkj(:),kkv(:),ndof,ndof);    
%%%%% loop over all elements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
      if ( is<=nlod )
        u_inc(:,2) =   dt * u_pre(:,2);   % scale dirichlet bc's with dt
        F_inc      = time * F_pre;        % scale neumann   bc's with dt
      else
	    u_inc(:,2) = 0.0;    
	  end
      R = R - F_inc;                    % add ext load to righthand side
      dofold = dof; 
      [dof,F] = solve_nr(K,R,dof,iter,u_inc);     % solution and  update     
      residuum=res_norm((dof-dofold),u_inc);
      disp(['*** iteration: ',num2str(iter),...
            ' *** residual: ',num2str(residuum)]);  
	  if(iter>20); disp(['*** NO GLOBAL CONVERGENCE ***']); return; end;		    
    end
%%% global newton-raphson iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   i_var = I;
   time = time + dt;
   plot_int(e_spa,i_var,nel,ndim,is);

%%% store restart data and overwrite after every step %%%%%%%%%%%%%%%%%
    f01=fopen('restart_time','w');fprintf(f01,'%i\n',is);fprintf(f01,'%f\n',mat);fclose(f01);
    f02=fopen('restart_data','w');fprintf(f02,'%f\n',dof);fclose(f02);
    f03=fopen('restart_ivar','w');fprintf(f03,'%f\n',i_var);fclose(f03);
end  
%%% loop over all load steps %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
nsteps = nsteps + inpstep;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% plot material configuration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    elseif (strcmp(macro(1:4),'pmat') == 1);
      e_mat = extr_dof(edof,q0);    % init material coordinates per elem
      figure(1)
      plot_mat(e_mat,ndof)
%%% plot material configuration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% plot internal variables  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    elseif (strcmp(macro(1:4),'pint') == 1);
      e_spa = extr_dof(edof,dof);   % init spatial  coordinates per elem
      figure(1)
      plot_int(e_spa,i_var,nel,ndim,is);
%%% plot internal variables  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% quit %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    elseif (strcmp(macro(1:4),'quit') == 1);
      return
    else
      disp('step    ... apply one load step')
      disp('step,,n ... apply  n  load steps')  
      disp('pmat    ... plot material configuration')
      disp('pspa    ... plot spatial  configuration')
      disp('quit    ... quit fe analyses')
    end
  end
end

hold off; fclose(f01); fclose(f02); fclose(f03);
%%% quit %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%








