function [out,out1,out2,out3] =fe_time(varargin)

%FE_TIME computes time integration with Newmark
%        and Discontinuous Galerkin schemes
%
%	Syntax : [def,model]=fe_time(com,model,Case,q0,opt)
%
%       Accepted commands are : newmark (damping supported), 
%                               dg (no damping supported),
%                               nlnewmark (damping supported) : newmark + newton 
%
%       com   :  is a string or a data structure
%
%        if com is string
%        You may enter simulation informations in the command using the format 
%              'Newmark (beta) (gamma) (t0) (deltaT) (Nstep) (Nf) (Thres)'
%        (beta, alpha, t0, deltaT, Nstep, Nf and Thres are described below)
%        In this case, you ommit the last input argument opt. For example :
%              [def,model]=fe_time('newmark .25 .5 0 1e-4 50',model,Case,q0); 
% 
%        if com is data structure, fields are : 
%                  .Method : 'nlnewmark', 'dg' or 'newmark'
%                  .Opt :  [beta alpha t0 deltaT  Nstep Nf Thres]
%                  .MaxIter : maximum number of iterations
%                  .Jacobian (only for non-linear). The default jacobian is
%                   'ki=ofact(model.K{3}+2/dt*model.K{2}+4/dt^2*m);'
%                  .JacobianUpdate (only for nlnewmark) : default is 0
%                     0 modified newmark (no update in Newton iterations)
%                     1 update in Newton iteration
%                  .Residual (only for non linear). The default residual is
%                      'r = model.K{1}*a+model.K{2}*v+model.K{3}*u-fc;'
%                  .InitAcceleration : optional field to specify
%                           a particular method to init
%                           the initial guess of acceleration field
%                  .OutputFcn  (for post-processing)
%                  .SaveTimes  : optionnal time vector, saves time steps
%                  .TimeVector (if exists TimeVector is used 
%                                instead of time parameters)  
%                  .RelTol (default getpref('OpenFEM','THRESHOLD',1e-6))
%
%       Other input arguments are :
%
%       model (see doc('model')), Case (see doc('case'))
%
%       q0    : initial conditions
%               data structure containing 1 or 2 fields
%               fiels are : .DOF specifying dofs
%                           .def and .v specifying displacements and velocity 
%                                       at t=0
%               if q0 is empty, zeros initial conditions are taken
%
%       See also : help fe_mk, fe_case
%                  doc  fe_time, fe_mk, fe_load, fe_case

%	Jean-Michel Leclere, Etienne Balmes
%       Copyright (c) 2001-2006 by INRIA and SDTools,All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       Use fe_time('cvs') for revision information


%xxx opt.MaxIter, opt.SaveFcn OuputFcn => hadle with displays

% - - - - - - - - - - - - - - - - - - -  - get command
if ischar(varargin{1})

 [CAM,Cam]=comstr(varargin{1},1);carg=2;
 opt=struct('Method','','Opt',[]);
 if comstr(Cam,'newmark')    
   opt.Method='Newmark'; [CAM,Cam]=comstr(CAM,8);
 elseif comstr(Cam,'cvs');
  out='$Revision: 1.69 $  $Date: 2006/05/10 12:04:01 $';return;
 elseif comstr(Cam,'nlnewmark') 
   opt.Method='NLnewmark'; [CAM,Cam]=comstr(CAM,10);
 elseif comstr(Cam,'dg')     
   opt.Method='dg';     [CAM,Cam]=comstr(CAM,3);

 %------------------------------------------------------------------------
 %------------------------------------------------------------------------
 %------------------------------------------------------------------------
 elseif comstr(Cam,'demo') [CAM,Cam]=comstr(CAM,5);
  
  if comstr(Cam,'bar')  
   femesh('reset');
   FEnode=[1 0 0 0    0   0  0
          2 0 0 0    50  0  0];
   FEel0=[Inf abs('bar1')];
   FEel0(2,:)=[1 2 1 1 1];
   pl=[1 1 2.1e10 0 2500 0];
   il=[1 1 0 0 0 1];
   femesh(';divide500');
   model=struct('Node',FEnode,'Elt',FEel0,'pl',pl,'il',il);
   % xxx use case
   data=struct('DOF',2.01,'def',1e6,'curve','input');
   model = fe_case(model,'DOFLoad','Point load 1',data);
   model=fe_curve(model,'set','input','step 1e-4*10');
   model = fe_case(model,'KeepDof','kd',.01);
   model = fe_case(model,'FixDof','fd',1);
   out=model; out1=fe_case(model,'getcase');
   return; 
  elseif comstr(Cam,'2d')
   FEnode=[1 0 0 0    0   0  0; 2 0 0 0    10  0  0];
   FEel0=[Inf abs('bar1')]; FEel0(2,:)=[1 2 1 1 1];
   pl=[1 fe_mat('m_elastic','SI',1) 1e10 0. 2500 0]; 
   il=[1 fe_mat('p_solid','SI',2) 1];
   femesh('extrude 1 0 10 0'); femesh(';divide100 100;quad2tria');
   femesh('set groupa1 name t3p matid1 proid1')

   model=struct('Node',FEnode,'Elt',FEel0,'pl',pl,'il',il);
   [nodeID,node]=feutil('findnode x==0 & y==0',model);
   data=struct('DOF',[nodeID+.01 nodeID+.02]','def',[1e6 1e6]');
   data.curve=struct('X',linspace(1,50e-4,10),'Y',ones(10,1));

   model = fe_case(model,'DOFLoad','Point load 1',data);
   model = fe_case(model,'FixDof','fdX','x==0 & y>0 -DOF  2');
   model = fe_case(model,'FixDof','fdY','y==0 & x>0 -DOF 1');
   out=model; return;    

  else error('fe_time demo : unknown');
  end
 %------------------------------------------------------------------------
 %------------------------------------------------------------------------
 %------------------------------------------------------------------------
 %------------------------------------------------------------------------
 elseif comstr(Cam,'optimproduct'); [CAM,Cam]=comstr(CAM,13);

  % [ind,flag]=fe_time('optimproduct',k); % flag=1 => transpose

  k=varargin{carg}; carg=carg+1;
  
  if comstr(version,'7');
   st1={'','symamd','colamd','colperm','symrcm'};
  else;st1={'','symamd','colmmd','colperm','symmmd','symrcm'};
  end

  st=''; t1=Inf; i3=10; u=ones(length(k),1);ut=u';

  for j1=1:length(st1);
   if isempty(st1{1,j1}); i2=1:length(k);else; i2=feval(st1{1,j1},k); end
   k2=k(i2,i2);
   while 1
    t=cputime; for j2=1:i3; ut*k2; end;st1{3,j1}=cputime-t;
    if st1{3,j1}>1; break; else; i3=i3*10;end
   end
   t=cputime; for j2=1:i3; k2*u; end;st1{2,j1}=cputime-t;
  end
  r2=[[st1{2,:}];[st1{3,:}]]';r3=min(r2(:));
  r2=(r2-ones(size(r2,1),2)*min(r3))/r3*100;
  [r3,i2]=min(r2(:)); 
  out1=0; if i2>size(r2,1); out1=1; i2=i2-size(r2,1);end % transpose
  j1=i2;for j1=1:size(r2,1); 
   fprintf('%10s %6.1f %6.1f %% slower\n',st1{1,j1},r2(j1,:));
  end
  if isempty(st1{1,j1}); i2=1:length(k);else; i2=feval(st1{1,j1},k); end
  out=i2; out2=st1{1,j1};
  return;
 
 %------------------------------------------------------------------------
 else error('not a known type of time integration computation');
 end
 opt.Opt=comstr(Cam,[-1]);
 if isempty(opt.Opt) & nargin>=4 opt.Opt=varargin{5}; end
 if length(opt.Opt)<5 error('bad format for com.Opt'); end 

elseif isstruct(varargin{1})

 if isfield(varargin{1},'Method') opt=varargin{1};carg=2;
 elseif sp_util('issdt'); 
   model=varargin{1};opt=fe_def('defTimeOpt',model);carg=1;
 else error('First input argument must specify a method');
 end
 

else error('First input argument must be a string or a data structure');
end

% - - - - - - - - - - - - - - - - - - -  - get model, case and q0
model=varargin{carg}; carg=carg+1;
if nargin>=carg
 Case = varargin{carg}; carg=carg+1;
 if isempty(Case)      Case=fe_case(model,'getcase');
 elseif ischar(Case)    Case=fe_case(model,['getcase' Case]);
 end
else Case=fe_case(model,'getcase');
end
if nargin>=carg  q0=varargin{carg}; carg=carg+1; 
else q0=[]; end

% validate data for various methods - - - - - - - - - - - - - - - - - - - - 
if strcmpi(opt.Method,'static') |  strcmpi(opt.Method,'nlnewmark')
  t=[]; if ~isfield(opt,'MaxIter') opt.MaxIter=100; end
end
if ~strcmpi(opt.Method,'static') 
 if ~isfield(opt,'TimeVector')
  dt=opt.Opt(4); t=opt.Opt(3)+dt*[0:opt.Opt(5)]; nt=length(t);
 else t=opt.TimeVector;  nt=length(t); end
 if ~isfield(model,'DOF') | ~isfield(model,'K')  
   model=fe_mk(model); 
 end
 if length(model.K)<3&isequal(model.Opt(2,:),[2 1])
   model.K={model.K{1}  ...
     spalloc(size(model.K{1},1),size(model.K{1},1),0) model.K{2}};
   model.Opt=[1 0 0 ;2 3 1];
 end
end

% assemble the model if necessary - - - - - - - - - - - - - - - - - - - - 
if isempty(Case.DOF) Case=fe_case(model,'gett',Case);end
mdof=Case.DOF;
% f(t)
opt.nt=nt; 
[fc,ft]=fe_load('buildu time',model,Case,t,opt);
try
catch
 fc=zeros(length(Case.DOF),1); ft=ones(nt,1);
 sdtw('_nb','Building load case : Failed');
end

if isempty(fc)
 if isfield(opt,'Residual') disp('fc=residual')
 else error('You must provide fc for linear solvers');
 end
elseif isfield(fc,'DOF') fc=fe_c(mdof,fc.DOF,fc.def')'; 
end

u=zeros(length(mdof),1);v=zeros(length(mdof),1);a=zeros(length(mdof),1);
opt1=opt.Opt;

% check initial conditions - - - - - - - - - - - - - - - - - - - - 

if isempty(q0) % zeros initial conditions 
   u=zeros(length(mdof),1);v=u;
elseif isstruct(q0) & isfield(q0,'def') & isfield(q0,'DOF') 
   q1(1:length(q0.DOF),1)=q0.def;
   if  isfield(q0,'v')  q1(1:length(q0.DOF),2)=q0.v; 
   else q1(1:length(q0.DOF),2)=0.; end

   % reorder CI dofs if necessary
   if (length(mdof)~=length(q0.DOF)  | norm(q0.DOF-mdof))
    q0=fe_c(mdof,q0.DOF,q0.def')';
   else
    q0=q1;
   end
   u=q0(:,1); if size(q0,2)==2; v=q0(:,2); else; v=zeros(size(u));end
elseif size(q0,1)==length(mdof)
   u=q0(:,1); if size(q0,2)==2; v=q0(:,2); else; v=zeros(size(u));end
else error('problem on initial conditions')
end

% Give standard names to matrices
try;
 m=model.K{find(model.Opt(2,:)==2)};
 k=model.K{find(model.Opt(2,:)==1)};
 if any(model.Opt(2,:)==3); c=model.K{find(model.Opt(2,:)==3)};
 else; c=spalloc(size(k,1),size(k,2),1);
 end
catch
 [m,c,k]=deal(model.K{1:3});
end

% initialize the output
tout=[];
if ~isfield(opt,'OutputFcn') | isempty(opt.OutputFcn)
   accel=zeros(length(mdof),nt);
   out=struct('DOF',mdof,'def',accel,'v',accel,'a',accel, ...
        'data',t(:),'fun',[0 4]);
   opt.OutputFcn='out.def(:,j1+1)=u;out.v(:,j1+1)=v;out.a(:,j1+1)=a;';
   if isfield(opt,'OutInd') 
    out.OutInd=int32(opt.OutInd); 
    out.cur=[0 t(1)];
    out.def=zeros(length(opt.OutInd),nt);
    out.v=zeros(length(opt.OutInd),nt);
    out.a=zeros(length(opt.OutInd),nt);
    out.DOF=out.DOF(opt.OutInd);
    opt.OutputFcn='out.def(:,j1+1)=u(opt.OutInd);out.v(:,j1+1)=v(opt.OutInd);out.a(:,j1+1)=a(opt.OutInd);';
   end
 elseif isa(opt.OutputFcn,'double') % Vector of times is given
  tout=opt.OutputFcn(:);tout=tout(find(tout>min(t)&tout<max(t(1:end-1))));
  out=struct('DOF',Case.DOF, 'def',zeros(length(Case.DOF),length(tout)), ...
             'data',tout,    'cur',[0 t(1)],   ...
             'OutInd',int32(1:length(Case.DOF)), ...
             'fun',[0 4]);
  if isfield(opt,'OutInd') 
   out.OutInd=int32(opt.OutInd); 
   out.def=zeros(length(opt.OutInd),length(tout));
   out.v=zeros(length(opt.OutInd),length(tout));
   out.a=zeros(length(opt.OutInd),length(tout));
   out.DOF=out.DOF(opt.OutInd);
  end
 elseif isfield(opt,'OutputInit') & ~isempty(opt.OutputInit)
   eval(opt.OutputInit);
 end
 cput0=cputime;
 uva=zeros(length(u),3); 

% ----------------------------------------------------------------------------
switch comstr(opt.Method,-27)
% ----------------------------------------------------------------------------
case 'newmark'

  if ~isfield(opt,'OutputFcn') | isempty(opt.OutputFcn)
   if isfield(opt,'OutInd') 
    opt.OutputFcn=...
      'of_time (''newmarkinterp'', out, beta,gamma,uva,a,t(j1),t(j1+1) )' ;
   end
  elseif isa(opt.OutputFcn,'double') % Vector of times is given
   toutc=t(1);
   opt.OutputFcn=...
     'of_time (''newmarkinterp'', out, beta,gamma,uva,a,t(j1),t(j1+1) )' ;
  end
  [opt,SaveTimes]=init_opt_linear(opt); 
 
  % beta = 0, gamma=.5 % Explicit
  % beta =.25, gamma=.5 % Classical values for implicit 

  % initial conditions
  if isfield(out,'OutInd');
   out.def(:,1)=u(out.OutInd,1);out.v(:,1)=v(out.OutInd,1);
   out.a(:,1)=a(out.OutInd,1); 
  else
   out.def(:,1)=u;out.v(:,1)=v;out.a(:,1)=a; 
  end
  ki=[];dt0=0.; beta=opt1(1); gamma=opt1(2);

  % explicit (centered difference method shown in Geradin/Rixen
  % figure 7.4.1)
  if isequal([beta gamma],[0 0.5]) 

   v12=v+(t(2)-t(1))/2*a;

   for j1=1:1:length(t)-1 % loop on times - - - - - - - - - - - - - - - - -

    tc=t(j1+1);dt=t(j1+1)-t(j1); 

    % update iteration matrix if necessary
    if abs(dt0/dt-1)>1e-6 | isempty(ki);
      ki=newmark_update(model,ki,dt,dt0,opt1);
    end

    % increment displacement and half time step
    u=u+dt*v12;
    tc12=tc-dt/2; 

    % update acceleration
    if isempty(opt.Residual); r = (ft(j1,:)*fc'-v'*c-u'*k)';
    else;    eval(opt.Residual);
    end
    a = ki\r;
    % increment velocity
    v12 = v12 + dt*a;% xxx use dt1/2
    v = v12 + (tc-tc12)*a;

    eval(opt.OutputFcn);  dt0=dt;
    if ~isempty(SaveTimes) & (tc>=SaveTimes(1)|j1==length(t)-1);
     if ~isfinite(norm(u,inf)); error('The simulation has diverged'); end
     SaveTimes(1)=[];eval(opt.SaveFcn);  
    end

   end % loop on times - - - - - - - - - - - - - - - - -

  else % implicit case

   for j1=1:1:length(t)-1 % loop on times - - - - - - - - - - - - - - - - -

    tc=t(j1+1);dt=t(j1+1)-t(j1);
    % update iteration matrix if necessary
    if abs(dt0/dt-1)>1e-6 | isempty(ki);
      ki=newmark_update(model,ki,dt,dt0,opt1);
    end

    of_time('storelaststep',uva,u,v,a);
    % predictions
    u = u + dt*v + (.5-beta)*dt^2*a; v = v + (1-gamma)*dt*a;

    % update acceleration
    if isempty(opt.Residual); r = (ft(j1,:)*fc'-v'*c-u'*k)';
    else;    eval(opt.Residual);
    end

    a = ki\r;
    % corrections
    u = u + beta*dt^2*a;  v = v + gamma*dt*a;
    eval(opt.OutputFcn);  dt0=dt;
    if ~isempty(SaveTimes) & (tc>=SaveTimes(1)|j1==length(t)-1);
     if ~isfinite(norm(u,inf)) error('The simulation has diverged'); end
     SaveTimes(1)=[];eval(opt.SaveFcn);  
    end

   end % loop on times - - - - - - - - - - - - - - - - -

  end % of explicit/implicit choice

  ofact('clear',ki);
  out=feutil('rmfield',out,{'cur','OutInd'});
  if cputime-cput0>30;
    fprintf('Elapsed CPU time : %i s\n',fix(cputime-cput0));
  end


% ----------------------------------------------------------------------------
% - - - - - - - - - - -  DG
case 'dg'


 if ~isfield(opt,'OutputFcn') | isempty(opt.OutputFcn)
   accel=zeros(length(mdof),nt);
   out=struct('DOF',mdof,'def',accel,'v',accel,'a',accel, ...
        'data',t(:),'fun',[0 4]);
   opt.OutputFcn='out.def(:,j1+1)=u;out.v(:,j1)=v;out.a(:,j1+1)=a;';
 elseif isfield(opt,'OutputInit') & ~isempty(opt.OutputInit)
   eval(opt.OutputInit);
 end
 if ~isfield(opt,'RelTol') opt.RelTol=getpref('OpenFEM','THRESHOLD',1e-6);end
 
 opt=init_opt_linear(opt);
 u1=u;v1=v; u2=u;v2=v; a=v*0; 
 % update force
  f1=(ft(1,:)*fc')';f2=(ft(1,:)*fc')';

 % initial conditions
 out.def(:,1)=u;out.v(:,1)=v;out.a(:,1)=a;  

 ki=[];dt0=0.;

 for j1=1:1:length(t)-1
  tc=t(j1+1);dt=t(j1+1)-t(j1); 
  ki=dg_update(model,ki,dt,dt0);   
  % update forces
  f1=(ft(j1,:)*fc')';f2=(ft(j1,:)*fc')';

  %  update RHS
  f1s=dt/2*f1-dt/6*f2+5/3*m*v2-2/3*dt*k*u2;
  f2s=dt/2*f1+dt/2*f2+m*v2-dt*k*u2;
  % multicorrection
  j2=0;
  while (1)
   % update residual and correction
   v1s=ki\(f1s-2/3*m*v2);v2s=ki\(f2s-dt^2/3*model.K{3}*v1);
   r1=norm(f1s-ki*v1s-2/3*m*v2s)/norm(f1s);
   r2=norm(f2s-dt^2/3*model.K{3}*v1s-ki*v2s)/norm(f2s);
   v1=v1s; v2=v2s;
   if (r1<opt.RelTol & r2<opt.RelTol)  break; end
   if (r1>1e6 | r2>1e6)  error(1); ;end
   j2=j2+1;
  end %while
  fprintf('Step %i number of iterations =  %i\n',j1,j2);
  % update displacement
  u1=u1+dt/6*(v1-v2); u2=u2+dt/2*(v1+v2);
  u=u2; v=v2; 
  eval(opt.OutputFcn);
  dt0=dt;
 end %for

% ----------------------------------------------------------------------------
% - - - - - - - - - - -  non linear statics

case 'static'

  opt=init_newton_calls(opt);
  [Case,model.DOF]=fe_mknl('init',model,Case);
  try
   Load=fe_load(model,Case);
  catch
   Load=zeros(length(Case.DOF),1);
  end

  opt.Residual='r=Load.def-k*u;';
  ki=[]; 
  opt.Jacobian='if ~isempty(ki) ofact(''clear'',ki);end;ki=ofact(k);'; 
  if ~isfield(opt,'RelTol') opt.RelTol=getpref('OpenFEM','THRESHOLD',1e-6);end

  r=ones(size(u)); 

  ite=0; while norm(r)>opt.RelTol & ite<opt.MaxIter

   fe_mknl('assemble',model,Case,300,u); % update state
   k=fe_mknl('assemble',model,Case,1);
   eval(opt.Jacobian); 
   eval(opt.Residual); 
   du=ki\r;
   u=u+du; 
   ite=ite+1;

  end

  error(1);


% ----------------------------------------------------------------------------
% - - - - - - - - - - -  implicit + Newton
case 'nlnewmark'

  if ~isfield(opt,'OutputFcn') | isempty(opt.OutputFcn)
   if isfield(opt,'OutInd') 
    opt.OutputFcn=...
      'of_time (''newmarkinterp'', out, beta,gamma,uva,a,t(j1),t(j1+1) )' ;
   end
  elseif isa(opt.OutputFcn,'double') % Vector of times is given
   toutc=t(1);
   opt.OutputFcn=...
     'of_time (''newmarkinterp'', out, beta,gamma,uva,a,t(j1),t(j1+1) )' ;
  end
  % initial conditions
  if isfield(out,'OutInd');
   out.def(:,1)=u(out.OutInd,1);out.v(:,1)=v(out.OutInd,1);
   out.a(:,1)=a(out.OutInd,1); 
  else
   out.def(:,1)=u;out.v(:,1)=v;out.a(:,1)=a; 
  end
  ki=[];dt0=0.; beta=opt1(1); gamma=opt1(2);


  [opt,SaveTimes]=init_newton_calls(opt);

  % initial acceleration
  if ~isfield(opt,'InitAcceleration') | isempty(opt.InitAcceleration)
   sdtw('_nb','Using default acceleration init k\(f-k*u-c*v)');
   k0=ofact(m); a=k0\((ft(1,:)*fc')'- model.K{2}*v - model.K{3}*u); 
   ofact('clear',k0);
  else                   eval(opt.InitAcceleration);
  end

  % Lookup figure 7.5.1 in Geradin/Rixen
  % initial conditions
  out.def(:,1)=u;
  ki=[]; dt0=0.; fc=full(fc); 

  for j1=1:length(t)-1

   tc=t(j1+1);dt=tc-t(j1); 
   of_time('storelaststep',uva,u,v,a);

   % predictions
   u = u + dt*v + (.5-opt1(1))*dt^2*a;
   v = v + (1-opt1(2))*dt*a;
   a = zeros(length(mdof),1);

   % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   % - - - - - - - - - - - - - -  Newton iterations => correction on u, v and a 
   [u,v,a,ki,opt] = iterNewton(ki,(ft(j1,:)*fc')',u,v,a,dt,dt0,tc,model,opt,Case,j1);
   % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   eval(opt.OutputFcn);  dt0=dt;
   if ~isempty(SaveTimes) & (tc>=SaveTimes(1)|j1==length(t)-1);
    if ~isfinite(norm(u,inf)) error('The simulation has diverged'); end
    SaveTimes(1)=[];eval(opt.SaveFcn);  
    fprintf('Time = %f            Newton iterations : %i \n',tc,opt.ite);
   end

  end % for j1


% ----------------------------------------------------------------------------
otherwise error(sprintf('''%s'' not know',opt.Method)); end

out1=model; 

% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function [u,v,a,ki,opt] = iterNewton(ki,fc,u,v,a,dt,dt0,tc,model, ...
  opt,Case,j1)

 if ~isfield(opt,'RelTol'); opt.RelTol=getpref('OpenFEM','THRESHOLD',1e-6);end
 % residu
 opt.nf=0; eval(opt.IterInit);
 if opt.nf==0 % default init of residual
   eval(opt.Residual); r0=full(r);
   if opt.nf==0;
    if norm(full(fc)); opt.nf=norm(fc); else; opt.nf=norm(r0); end
   end
 end % default init of residual

 j2=0;opt1=opt.Opt;
 while(1)

  % update iteration matrix
  if j2==0 | opt.JacobianUpdate~=0  eval(opt.Jacobian); end

  % compute corrections
  dq = ki\r;

  % corrections
  u = u - dq;
  v = v - (opt1(2)/opt1(1))/dt*dq;
  a = a - 1/opt1(1)/dt^2*dq;

  % update residual
  eval(opt.Residual);

  % convergence test
  j2=j2+1;
  if norm(r)/opt.nf<opt.RelTol ;break; end
  if j2>abs(opt.MaxIter); sdtw('Max number of iterations reached');break; 
  elseif j2>abs(opt.MaxIter) & opt.MaxIter<0
   error('Max number of iterations reached');
  end 
  
 end % while
 opt.ite=j2;eval(opt.IterEnd);


% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function  [opt,SaveTimes]=init_newton_calls(opt);

  if ~isfield(opt,'Jacobian') |  isempty(opt.Jacobian)
    opt.Jacobian='ki=basic_jacobian(model,ki,dt,dt0,opt.Opt);'; 
  end
  if ~isfield(opt,'Residual') | isempty(opt.Residual)
    opt.Residual='r = model.K{1}*a+model.K{2}*v+model.K{3}*u-fc;'; 
  end
% r = (ft(j1,:)*fc'-v'*c-u'*k)';
  if ~isfield(opt,'JacobianUpdate') opt.JacobianUpdate=0; end
  if isfield(opt,'SaveTimes') 
   pw0=evalin('base','pwd');  SaveTimes=opt.SaveTimes;
   if ~isfield(opt,'SaveFcn');
    opt.SaveFcn=sprintf('save %s out', ...
       fullfile(pw0,'fe_time_save.mat'));
   end
  else; opt.SaveTimes=[]; SaveTimes=[];
  end
  if ~isfield(opt,'IterInit'); opt.IterInit='';end
  if ~isfield(opt,'IterEnd');  opt.IterEnd='';end

% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function  ki=basic_jacobian(model,ki,dt,dt0,opt);

if abs(dt0/dt-1)>1e-10 | (isnumeric(ki) &isempty(ki))
  ofact('clear',ki);
  ki=ofact(model.K{3}+(opt(2)/opt(1))/dt*model.K{2} + 1/opt(1)/dt^2*model.K{1});
end

% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function [opt,SaveTimes]=init_opt_linear(opt)

 if isfield(opt,'Jacobian')  & ~isempty(opt.Jacobian)
    error('Jacobian field not taken into account for linear problems');
 end
 if ~isfield(opt,'Residual') | isempty(opt.Residual)
    opt.Residual='';% 'r=(fc*ft(j1)-model.K{2}*v-model.K{3}*u);';
 end
 if isfield(opt,'SaveTimes') 
  pw0=evalin('base','pwd');  SaveTimes=opt.SaveTimes;
  if ~isfield(opt,'SaveFcn');
   opt.SaveFcn=sprintf('save %s out', ...
      fullfile(pw0,'fe_time_save.mat'));
  end
 else; opt.SaveTimes=[]; SaveTimes=[];
 end


% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function  ki=newmark_update(model,ki,dt,dt0,opt);   

if norm((dt-dt0)/dt)>1e-10 | (isnumeric(ki) &isempty(ki))
 ofact('clear',ki);
 % if beta is zero and C is not diagonal just keep M
 % this is to favor speed in explicit integration
 if opt(1)==0 & ~isequal(model.K{2},diag(diag(model.K{2})))
   kiter=model.K{1}; kiter=kiter+spalloc(size(kiter,1),size(kiter,2),0);
 else
                   %gamma                %beta
   kiter=model.K{1}+opt(2)*dt*model.K{2}+opt(1)*dt^2*model.K{3};
 end
 if nnz(kiter)==size(kiter,1)&isequal(kiter,diag(diag(kiter)))
  ki=ofact(kiter,'diag');
 else
  ki=ofact(kiter);
 end
end

% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function  ki=dg_update(model,ki,dt,dt0);   

if norm((dt-dt0)/dt)>1e-10 | (isnumeric(ki) &isempty(ki))
 ki=(model.K{1}+dt^2/6*model.K{3});
end

% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
function ft=get_curve_fetime(model,curve)

[st,i2]=stack_get(model,'curve',curve);
if isempty(i2) error('no curve found in model Stack')
else
 %xxx check values of t
 ft=st{3}.data;
end


% ----------------------------------------------------------------------------












