function [out,out1,out2]=basic_elt_test(varargin);

% ELEMENT TESTS
%
%% eig, load, mat for all elements => comparisons with SDT5.1
% basic_elt_test('compare')
%
%% test of basic commands (integinfo, call, etc.) for all elements
% basic_elt_test('calls')
% basic_elt_test('all')
% basic_elt_test('eig')
% basic_elt_test('mat')
% basic_elt_test('load')
% basic_elt_test('dilat') % checks stresses for uniform dilatation
% out = basic_elt_test('load energy')



%% teststress for all elements, 
% out = basic_elt_test('stress')
%
%% testenergy for all elements, 
%
%% test mat,eig,surf,... for element
%%    st=q4p,t3p,...
%%   st1=mat,eig,sur,vol
% basic_elt_test(st,st1) or basic_elt_test (all the tests)
%
%   example
%   basic_elt_test('penta15','mat') % check anisotropic material for 2D

%	Etienne Balmes, J.M. Leclere
%       Copyright (c) 2001-2005 by INRIA and SDTools,All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       $Revision: 1.25 $  $Date: 2006/03/24 14:01:50 $

carg=1; if carg<=nargin; [CAM,Cam]=comstr(varargin{carg},1);carg=carg+1;
else
 Cam='all';
end
%delete basic_elt_test.log

% List of available tests - - - - - - - - - - - - - - - - - - - - - - - 
 list={ ...
     'q4p','eig_0','load_0','mat_0','visc',
     'q8p','eig_0','load_0','mat_0','visc',
     't3p','eig_0','load_0','mat_0','visc',
     't6p','eig_0','load_0','mat_0','visc',
     'q4p','eig_1','load_1','mat_1','visc',
     'q8p','eig_1','load_1','mat_1','visc',
     't3p','eig_1','load_1','mat_1','visc', 
     't6p','eig_1','load_1','mat_1','visc',
     'q4p','eig_2','load_2','mat_2','visc',
     'q8p','eig_2','load_2','mat_2','visc',
     't3p','eig_2','load_2','mat_2','visc', 
     't6p','eig_2','load_2','mat_2','visc',
     'hexa8','eig','load','mat','visc',
     'hexa20','eig','load','mat','visc',
     'penta6','eig','load','mat','visc',
     'penta15','eig','load','mat','visc',
     'tetra4','eig','load','mat','visc',
     'tetra10','eig','load','mat','visc',
     'q4pb','eig_0','load_0','mat_0','visc',
     'q8pb','eig_0','load_0','mat_0','visc',
     't3pb','eig_0','load_0','mat_0','visc',
     't6pb','eig_0','load_0','mat_0','visc',
     'q4pb','eig_1','load_1','mat_1','visc',
     'q8pb','eig_1','load_1','mat_1','visc',
     't3pb','eig_1','load_1','mat_1','visc', 
     't6pb','eig_1','load_1','mat_1','visc',
     'hexa8b','eig','load','mat','visc',
     'hexa20b','eig','load','mat','visc',
     'penta6b','eig','load','mat','visc',
     'penta15b','eig','load','mat','visc',
     'tetra4b','eig','load','mat','visc',
     'tetra10b','eig','load','mat','visc',
     'tria3','eig','load','mat','visc',
     'tria6','','','','visc',
     'quad4','eig','load','mat','visc',
     'quadb','eig','','mat','visc',
     'quad9','','','','visc',
     'bar1','eig','load','mat','visc',
     'flui4','eig','','mat','visc',
     'flui6','eig','','mat','visc',
     'flui8','eig','','mat','visc',
     'beam1','eig','load','mat','visc',
     'beam1t','eig','','','visc',
     };
% List of available tests - - - - - - - - - - - - - - - - - - - - - - - 
if comstr(Cam,'list'); out=list; return; end
%-----------------------------------------------------------------------
% testmat for all elements

if ~isempty(strfind(Cam,'enerk')); flag='enerk';
else; flag='';
end
if ~isempty(strfind(Cam,'stress')); flag='stress';
else; flag='';
end
if ~isempty(strfind(Cam,'visc')); flag='visc';
else; flag='';
end

if comstr(Cam,'mat');     jTest=4;
elseif comstr(Cam,'eig'); jTest=2;
elseif comstr(Cam,'load energy'); jTest=5;
elseif comstr(Cam,'load'); jTest=3;
else; jTest=[];
end

if ~isempty(jTest) % basic series of tests (mat eig load)

out={};
for j1=1:size(list,1);
 if ~isempty(list{j1,jTest})
  try;
   if comstr(Cam,'mat')
     K=[];
     st=sprintf('K=%s(''test%s'');',list{j1,[1 4]}); eval(st);
     out(end+1,1:2)={st,K};
   elseif comstr(Cam,'eig')
     st=sprintf('[model,def]=%s(''test%s%s'');',list{j1,[1 2]},flag); eval(st);
     out(end+1,1:3)={st,def,def.data(1:5)};
   elseif comstr(Cam,'load')
     st=sprintf('[model,def]=%s(''test%s%s'');',list{j1,[1 3]},flag); eval(st);
     out(end+1,1:3)={st,def,norm(full(def.def),2)};
   elseif comstr(Cam,'visc')
     st=sprintf('model=(''teststruct%s -back'');',list{j1,1}); eval(st);
     [mat,mdof] = fe_mk(model,'options',[3 ] )
     [Case,model.DOF]=fe_mknl('init',model);
     mat1=fe_mknl('assemble',model,Case,3);
     out(end+1,1:3)={st,{mat,mat},{size(find(mat)),size(find(mat1))}};
   end
  catch
   out(end+1,1:2)={st,'failed'};
  end  
 end
end

%-----------------------------------------------------------------------
% Verification of dilatation for volume elements
elseif comstr(Cam,'dilat')


 % 3D  : Uniform dilatation check

 list={'hexa8b rand.1','hexa8 rand.1','hexa20 rand.1', ...
       'tetra4b rand.1'}';list=list(:);
 for j1=1:length(list)
  model=femesh(horzcat(';reset;test',list{j1,1},'struct back'));
  model.pl=m_elastic('dbval 100 steel','dbval 112 aluminum');
  model.pl(:,4)=.285;model.pl(2,6)=model.pl(2,3)/2/(1.285);
  nind=sparse(model.pl(:,1),1,1:size(model.pl)); mpid=feutil('mpid',model);
  def=feutil('geomdilat',model,model.pl(nind(mpid(2,1)),4));
  s1=fe_stres('stress sxx atcenter',model,def);
  r1=norm([s1.data./model.pl(nind(mpid(s1.EltId+1,1)),3)]-1);
  s1=fe_stres('stress syy atcenter',model,def);  r1(2)=norm(s1.data);
  list(j1,2:3)={r1(1) r1(2)};
 end
 out=list;

 % 2D  : Uniform dilatation check

 list={'q4pb rand.3','q4p rand.3'};list=list(:);
 for j1=1:length(list)
  model=femesh(horzcat(';reset;test',list{j1,1},'divide 4 4 struct back'));
  model.pl=m_elastic('dbval 100 steel','dbval 112 aluminum');
  model.pl(:,4)=.285;model.pl(2,6)=model.pl(2,3)/2/(1.285);
  nind=sparse(model.pl(:,1),1,1:size(model.pl)); mpid=feutil('mpid',model);
  def=feutil('geomdilat',model,model.pl(nind(mpid(2,1)),4));
  s1=fe_stres('stress sxx atcenter',model,def);
  r1=norm([s1.data./model.pl(nind(mpid(s1.EltId+1,1)),3)]-1);
  s1=fe_stres('stress syy atcenter',model,def);  r1(2)=norm(s1.data);
  list(j1,2:3)={r1(1) r1(2)};
 end
 out(end+[1:size(list,1)],1:size(list,2))=list;


%-----------------------------------------------------------------------
% 
elseif comstr(Cam,'compare')

 st0=list;
 ref_elt_test;
 [st0,i1]=intersect(st,st0);
 
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Matrices
 Mat  = basic_elt_test('mat'); result_Mat=Mat;
 for j1=1:size(Mat,1)
  i1=find(strcmp(Mat_ref,Mat{j1,1}));
  if ~isempty(i1) 
   if  ~ischar(Mat_ref{i1,2}) & ~ischar(Mat{j1,2}) & isa(Mat{j1,2},'cell')
    try
     result_Mat{j1,2}=[(Mat_ref{i1,2}(1)-max(svd(Mat{j1,2}{1})))/Mat_ref{i1,2}(1) ...
                  (Mat_ref{i1,2}(2)-max(svd(Mat{j1,2}{2})))/Mat_ref{i1,2}(1)  ];
    catch
      result_Mat{j1,2}='Infinite result';
    end
   else
    result_Mat{j1,2}='error';
   end
  else
    result_Mat{j1,2}='not compared';
  end
 end % j1
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - EIG
 Eig  = basic_elt_test('eig');result_Eig=Eig;
 for j1=1:size(Eig,1)
  i1=find(strcmp(Eig_ref,Eig{j1,1}));
  if ~isempty(i1) 
   if  ~ischar(Eig_ref{j1,2}) & ~ischar(Eig{j1,2}) 
    result_Eig{j1,2}=Eig_ref{i1,2}(:)'-Eig{i1,3}(:)';
   else
    result_Eig{j1,2}='error';
   end
  else
    result_Eig{j1,2}='not compared';
  end
 end % j1
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - LOAD
 Load = basic_elt_test('load');result_Load=Load;
 for j1=1:size(Load,1)
  i1=find(strcmp(Load_ref,Load{j1,1}));
  if ~isempty(i1) 
   if  ~ischar(Load_ref{j1,2}) & ~ischar(Load{j1,2})
    result_Load{j1,2}=Load_ref{i1,2}-Load{i1,3};
   else
    result_Load{j1,2}='error';
   end
  else
    result_Load{j1,2}='not compared';
  end
 end % j1
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - STRESS
 Stress = basic_elt_test('stress');result_Stress=Stress;
if 1==2 %xxx
 for j1=1:size(Stress,1)
  i1=find(strcmp(Stress_ref,Stress{j1,1}));
  if ~isempty(i1) 
   if  ~ischar(Stress_ref{j1,2}) & ~ischar(Stress{j1,2})
    result_Stress{j1,2}=Stress_ref{i1,2}-Stress{i1,3};
   else
    result_Stress{j1,2}='error';
   end
  else
    result_Stress{j1,2}='not compared';
  end
 end % j1
end
 
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

 out={result_Mat,result_Eig,result_Load};

 % - - - - - - - - - - - - - - - - - - - - - - - - - - -  Comparisons output
 for j1=1:size(out{1},1) % Mat
  if ~ischar(out{1}{j1,2}) & ~isempty(find(out{1}{j1,2}>sdtdef('epsl')))
    fprintf('%s : ',out{1}{j1,1});
    fprintf('%f ',out{1}{j1,2});
    fprintf('\n');
  elseif ~ischar(out{1}{j1,2}) & isempty(find(out{1}{j1,2}>sdtdef('epsl')))
    %fprintf('%s : comparison OK\n',out{1}{j1,1});
  else fprintf('%s : error\n',out{1}{j1,1});
  end
 end

 for j1=1:size(out{2},1) % Eig
  if ~ischar(out{2}{j1,2}) & ~isempty(find(out{2}{j1,2}>sdtdef('epsl')))
    fprintf('%s : ',out{2}{j1,1});
    fprintf('%f ',out{2}{j1,2});
    fprintf('\n');
  elseif ~ischar(out{2}{j1,2}) & isempty(find(out{2}{j1,2}>sdtdef('epsl')))
    %fprintf('%s : comparison OK\n',out{2}{j1,1});
  else fprintf('%s : error\n',out{2}{j1,1});
  end
 end

 for j1=1:size(out{3},1) % Load
  if ~ischar(out{3}{j1,2}) & ~isempty(find(out{3}{j1,2}>sdtdef('epsl')))
    fprintf('%s : ',out{3}{j1,1});
    fprintf('%f ',out{3}{j1,2});
    fprintf('\n');
  elseif ~ischar(out{3}{j1,2}) & isempty(find(out{3}{j1,2}>sdtdef('epsl')))
    %fprintf('%s : comparison OK\n',out{3}{j1,1});
  else fprintf('%s : error\n',out{3}{j1,1});
  end
 end
 return;

%-----------------------------------------------------------------------
% test each basic command of elements
elseif comstr(Cam,'call')

 st1={'call','rhscall','node','patch',...
      'dof','line','face','parent'};

 for j1=1:length(list)
  for j2=1:length(st1)
    st3=['r1=' list{j1,1} '(''' st1{j2} ''');'];
    try ;
     eval(st3)
     if isempty(r1); r1='Empty';else; r1=''; end
    catch
     r1=lasterr;
    end
  list{j1,j2+1}=r1;
  end %j2
 end %j1

 
 j2=2;while j2<=size(list,2) % empty columns
  r1=char(list(:,j2)); 
  if isempty(r1) list(:,j2)=[];st1(j2-1)=[];
  else; j2=j2+1;
  end
 end
 j2=1;while j2<=size(list,1) % empty rows
  r1=strcat(list{j2,2:end}); 
  if isempty(r1) list(j2,:)=[];
  else; j2=j2+1;
  end
 end

 st1={'' st1{:}};st1(end+[1:size(list,1)],:)=list;
 out=st1;

%-----------------------------------------------------------------------
% test each basic command of elements
elseif comstr(Cam,'all')

st='_';st=st(ones(1,60));
out=cell(1,4);

for j2=1:2 % compute then display
 %fprintf('\n%s\n\n Checking calls \n%s\n',st,st);
 %if isempty(out{1});out{1}=basic_elt_test('call');end
 %disp(out{1})

 fprintf('\n%s\n\n Checking matrix tests \n%s\n',st,st);
 if isempty(out{2});out{2}=basic_elt_test('mat');end
 ind=[];
 for j1=1:size(out{2},1)
  if ~iscell(out{2}{j1,2}); disp(out{2}{j1,1});end
 end

 fprintf('\n%s\n\n Checking eig tests \n%s\n',st,st);
 if isempty(out{3});out{3}=basic_elt_test('eig');end
 ind=[];
 for j1=1:size(out{3},1)
  if isequal(out{3}{j1,2},'failed'); disp(out{3}{j1,1});end
 end

 fprintf('\n%s\n\n Checking load tests \n%s\n',st,st);
 if isempty(out{4});out{4}=basic_elt_test('load');end
 ind=[];
 for j1=1:size(out{4},1)
  if isequal(out{4}{j1,2},'failed'); disp(out{4}{j1,1});end
 end

end % compute then display


%clear out

else error('Not a valid test')
end

% display the results more clearly - - - - - - - - - - - - - - - -
for j1=1:size(out,1);
 if isequal(out{j1,2},'failed');
  fprintf('Failed : %s\n',out{j1,1});
 end
end

return;


%-----------------------------------------------------------------------
% Actually do the test for one specified element

st2=sprintf('%s test%s;',st,st1);

try ; eval(st2)
catch
    disp(lasterr);
    fprintf(' ------------------ error in %s test\n%s\n', st2)
end

if any(strcmp(st,{'q4p','q5p','q8p','t3p','t6p','q9a'}))
 st2=sprintf('%s test%s_0;',st,st1);
 try ; eval(st2)
 catch
    disp(lasterr);
    fprintf(' ------------------ error in %s test\n%s\n',st2)
 end
 st2=sprintf('%s test%s_1;',st,st1);
 try ; eval(st2)
 catch
    disp(lasterr);
    fprintf(' ------------------ error in %s test\n%s\n',st2)
 end
 st2=sprintf('%s test%s_2;',st,st1);
 try ; eval(st2)
 catch
    disp(lasterr);
    fprintf(' ------------------ error in %s test\n%s\n',st2)
 end

 % test anisotropic material
 % plane stress
 E=2.1e11;nu=.285;C=E/(1.-nu*nu);
 e=[C C*nu C 0. 0. C*(1-nu)/2];
 pl=[100 fe_mat('m_elastic','SI',4) e 7800 0 0 0 0 .1];
 st2=sprintf('k=%s (''test%s_0'');',st,st1)
 eval(st2);
 st2=sprintf('k1=%s (''test%s_0'',pl);',st,st1)
 eval(st2);
 if isa(k1,'cell')
   if norm(k{1}-k1{1}) error(['anisotropic element ' st2]); end
 end

 % plane strain
 unmnu=1.-nu;
 C=E*unmnu/(1+nu)/(1-2*nu);
 e=[C nu*C/unmnu C 0. 0. C*(1-2*nu)/(2*unmnu)];
 pl=[100 fe_mat('m_elastic','SI',4) e 7800 0 0 0 0 .1];
 st2=sprintf('k=%s (''test%s_1'');',st,st1);
 eval(st2);
 st2=sprintf('k1=%s (''test%s_1'',pl);',st,st1);
 eval(st2);
 if isa(k1,'cell')
   if norm(k{1}-k1{1}) error(['anisotropic element ' st2]); end
 end



end

return;

%-----------------------------------------------------------------------
if 1==2
 % anisotropic
  model=femesh(';reset;testq4p');
  model.il(3)=0;

  E=2.1e11;nu=.285;C=E/(1.-nu*nu);
  e=[C C*nu C 0. 0. C*(1-nu)/2];
  model.pl=[100 fe_mat('m_elastic','SI',4) e 7800 0 0 0 0 .1];

  [constit,iopt,elmap]=q4p('integinfo',[100;110],model.pl,model.il);
  [k,m]=q4p(model.Node,model.Elt(2,:),[36 36 0 0 0 0 0 0 0],int32(iopt),constit,elmap);


     unmnu=1.-nu;
     C=E*unmnu/(1+nu)/(1-2*nu);
     e=[C nu*C/unmnu C 0. 0. C*(1-2*nu)/(2*unmnu)];
pl=[100 fe_mat('m_elastic','SI',4) e 7800 0 0 0 0 .1];
end





%--------------------------------------------------------- RHS
if 1==2
% st={'q4p','q5p','q8p','t3p','t6p',...
%      'hexa8','hexa20','penta6',...
%     'tetra4','tetra10'};

% st={'quad4','tria3','quadb'};

%st={'bar1','beam1'};
%st={'flui4','flui6','flui8'};

 st={'q4p','q5p','q8p','t3p','t6p'};
 
errors={}; 

for j2=1:4

 for j1=1:length(st)

    st1=['[model,def1]=femesh('';teststruct' st{j1} 'load'');' ];
    st2=['[model,def2]=femesh('';teststruct' st{j1} 'load'');' ];
    st3=['[model,def3]=femesh('';teststruct' st{j1} 'load'');' ];

    femesh('reset');eval(st1);
    femesh('reset');eval(st2);
    femesh('reset');eval(st3);

    if ~isequal(def1.def,def2.def) | ~isequal(def1.def,def3.def)
     errors{end+1,1}=st1;
     if ~isequal(def1.def(:,1:2),def2.def(:,1:2)) |...
                            ~isequal(def1.def(:,1:2),def3.def(:,1:2))
       errors{end,2}='gravity and/or surf';
     else
       errors{end,2}='pressure only';
     end
    end 
 end %j1
end %j2



end

%------------------------------------------ compare basic matrices
if 1==2
out=basic_elt_test('mat')
load basic_elt_test_matrices
result=out;

for j1=1:size(out,1)
  i1=find(strcmp(k,out{j1,1}));

  if ~isempty(i1) 
   if  ~ischar(out{j1,2});
    result{j1,2}=norm(out{j1,2}{1}-k{i1,2}{1});
    result{j1,3}=norm(out{j1,2}{2}-k{i1,2}{2});
   else
    result{j1,2}='error';
    result{j1,3}='error';
   end
  else
    result{j1,2}='not compared';
  end
end

end
%--------------------------------------------------------- RHS
if 1==2
out=basic_elt_test('load')
load basic_elt_test_load
result=out;

for j1=1:size(out,1)
  i1=find(strcmp(k,out{j1,1}));

  if ~isempty(i1) 
   if  ~ischar(out{j1,2});
    result{j1,2}=norm(out{j1,3}-k{i1,3});
   else
    result{j1,2}='error';
   end
  else
    result{j1,2}='not compared';
  end
end

end
%--------------------------------------------------------- EIG
if 1==2
out=basic_elt_test('eig')
load basic_elt_test_load
result=out;

for j1=1:size(out,1)
  i1=find(strcmp(k,out{j1,1}));

  if ~isempty(i1) 
   if  ~ischar(out{j1,2});
    result{j1,2}=norm(out{j1,3}-k{i1,3});
   else
    result{j1,2}='error';
   end
  else
    result{j1,2}='not compared';
  end
end

end
%---------------------------------------------------------
