% ----------------------------------------------------------------------------
%        simple test, bar  traction-compression et echelon
if 1==2
 FEnode=[1 0 0 0    0   0  0
         2 0 0 0    50  0  0];
 FEel0=[Inf abs('bar1')]; %femesh('testbar1 divide 100')
 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);
 model = fe_case(model,'AddToCase 1','DOFLoad','Point load 1',data);
 model = fe_case(model,'AddToCase 1','KeepDof','kd',.01);
 model = fe_case(model,'AddToCase 1','FixDof','fd',1);

 model=fe_mk(model);
 model.K{3}=model.K{2};
 model.K{2}=0.;
 mdof=model.DOF;

 %q0=zeros(length(mdof));
 q0=[];
 com.Method='newmark';
 com.Opt=[.25 .5 3e-4 1e-4 50 10];
 com.Residual='';
 com.OutputFcn='';
 %com.InputVector=[ones(1,10) zeros(1,40)];

 [def,model]=fe_time(com,model,'Case 1',q0);

if 1==2
 [def,model]=fe_time('dg',model,'Case 1',q0,[.25 .5 0 1e-4 50 10]);
 [def,model]=fe_time('newmark',model,'Case 1',q0,[.25 .5 0 1e-4 50 10]);
end
 cf=feplot;cf.model=model;cf.def={def.v,mdof+.01};fecom(';view2;animtime');
end
% ----------------------------------------------------------------------------
% same test with outputFcn
 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);
 data=struct('DOF',2.01,'def',1e6);
 model = fe_case(model,'AddToCase 1','DOFLoad','Point load 1',data);
 model = fe_case(model,'AddToCase 1','KeepDof','kd',.01);
 model = fe_case(model,'AddToCase 1','FixDof','fd',1);

 model=fe_mk(model); model.K{3}=model.K{2}; model.K{2}=0.;
 mdof=model.DOF;

 q0=[];
 com.Method='newmark';
 com.Opt=[.25 .5 3e-4 1e-4 50 10];
 com.Residual='';


 com.OutInd=1;
 [def1,model]=fe_time(com,model,'Case 1',q0); % reference

 dt=com.Opt(4); t=com.Opt(3)+dt*[0:com.Opt(5)]; 
 com.OutputFcn=t(1:end)+dt;
 [def2,model]=fe_time(com,model,'Case 1',q0); %xxx
 figure(1); plot(def2.data,def2.def(1,:),'o-');
 error('newmark interp problem')


 % test of_time newmarkinterp
 dt=com.Opt(4); t=com.Opt(3)+dt*[0:com.Opt(5)]; 
 com.OutputFcn=t(1:end)-dt/4;
 com.OutInd=1;
 [def,model]=fe_time(com,model,'Case 1',q0);

figure(1); plot(def.data,def.def(1,:),'o-');
com.OutputFcn=t(1:end)+dt/4;
[def,model]=fe_time(com,model,'Case 1',q0);
figure(1); plot(def.data,def.def(1,:),'o-');


% ----------------------------------------------------------------------------
% 2D test
 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 1 2.1e10 0.3 2500 0];
 il=p_shell('default');il=il.il;
 femesh('extrude 1 0 10 0');
 femesh(';divide100 100');

 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]');
 model = fe_case(model,'AddToCase 1','DOFLoad','Point load 1',data);
 model = fe_case(model,'AddToCase 1','KeepDof','kd',[.01 .02]');
 [nodeID,node]=feutil('findnode x=0 & y>0',model);
 model = fe_case(model,'AddToCase 1','FixDof','fdX',nodeID+.02);
 [nodeID,node]=feutil('findnode y=0 & x>0',model);
 model = fe_case(model,'AddToCase 1','FixDof','fdY',nodeID+.01);

 model=fe_mk(model);
 model.K{3}=model.K{2};
 model.K{2}=sparse(size(model.K{1},1),size(model.K{1},1));
 mdof=model.DOF;

 %q0=zeros(length(mdof));
 q0=[];
 com.Method='newmark';
 com.Opt=[.25 .5 3e-4 1e-4 50 10];
 com.Residual='';
% com.OutputFcn='';
% com.InputVector=[ones(1,10) zeros(1,40)];

 [def,model]=fe_time(com,model,'Case 1',q0);

if 1==2
 [def,model]=fe_time('dg',model,'Case 1',q0,[.25 .5 0 1e-4 50 10]);
 [def,model]=fe_time('newmark',model,'Case 1',q0,[.25 .5 0 1e-4 50 10]);
end

cf=feplot;cf.model=model;cf.def=def;fecom(';view2;animtime');
fecom colordataa

% ----------------------------------------------------------------------------
% to test extrapolation of curves
  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','ricker_1');
   model = fe_case(model,'AddToCase 1','DOFLoad','Point load 1',data);

   curve_1=fe_curve('TestRicker 10e-4 120 3. 1e-3') % here curve to extrapolate
   model=fe_curve(model,'set','ricker_1',curve_1)

   model = fe_case(model,'AddToCase 1','KeepDof','kd',.01);
   model = fe_case(model,'AddToCase 1','FixDof','fd',1);
   model=fe_mk(model);
   Case=fe_case(model,'getcase');

   def=fe_time('newmark .25 .5 0 1e-4 100',model,Case,[]);
   def_v=def;def_v.def=def_v.v; def_v.DOF=def.DOF+.01;
   feplot(model,def_v); fecom(';view2;animtime;ch30;scd3');
% ----------------------------------------------------------------------------
