function [out,out1,out2]=m_elastic(varargin)

%M_ELASTIC Elastic material function
%
%       Syntax : mat= m_elastic('default') 
%                mat= m_elastic('database name') 
%                pl = m_elastic('dbval MatId name');
%                pl = m_elastic('dbval -unit MM MatId name');
%
%       Material subtypes supported by m_elastic are
%
%	Subtype 1 : standard isotropic materials are described by a row
%	    [MatID type  E  nu rho G  eta  alpha T0]
%         with
%           type  fe_mat('m_elastic','SI',1)
%	    E     (Young's modulus)       nu    (Poisson's ratio)
%	    rho   (density)               G     (shear modulus)
%
%	Subtype 2 : acoustic fluid
%           [MatId type Rho C eta]
%           type  fe_mat('m_elastic','SI',2)
%           Rho   (Density), C (velocity), eta (loss factor)
%
%	Subtype 3 : 3-D anisotropic materials
%           [MatId type G11 G12 G22 G13 ... G14 ... G66 rho A1 ... A6 T0 eta]
%           type  fe_mat('m_elastic','SI',3)
%           Note : shear is ordered g_yz, g_zx, g_xy
%
%	Subtype 4 : 2-D anisotropic materials
%           [MatId type E11 E12 E22 E13 E23 E33 rho a1 a2 a3 T0 eta ...
%                                      ST SC SS MCSID]
%           type  fe_mat('m_elastic','SI',4)
%
%	Subtype 5 : Orthotropic material for shell
%           [MatId type E1 E2 nu12 G12 G1z G2z Rho A1 A2 TREF Xt Xc Yt Yc S ...
%             Ge F12 STRN]
%           type  fe_mat('m_elastic','SI',4)
%
%       See also help fe_mat, p_shell, p_beam
%                doc m_elastic, pl, fem


%	Etienne Balmes, Jean-Michel Leclere, Corine Florens
%       Copyright (c) 2001-2005 by INRIA and SDTools, All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license

if nargin==1 & comstr(varargin{1},'cvs')
 out='$Revision: 1.33 $  $Date: 2006/04/27 13:42:18 $'; return;
end

if nargin<1 help m_elastic;return; end
if ischar(varargin{1}) [CAM,Cam]=comstr(varargin{1},1); pl=[]; carg=2;
else pl=varargin{1};[CAM,Cam]=comstr(varargin{2},1); carg=3;
end

sdtdef('diag',['m_elastic(''' Cam ''')'])





% -------------------------------------------------------------------------
if comstr(Cam,'info')

 r1=m_elastic('database');
 disp(char({r1.name}))

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

  out=m_elastic('database'); out=out(1);

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


 while 1==1
  i1=strfind(comstr(Cam,-27),'-unit'); out1={}; % See if unit specified
  if ~isempty(i1)
   [Unit,i2,i3,i4]=sscanf(CAM(i1+5:end),'%s',1);
   i4=i1+[0:4+i4];CAM(i4)=''; [CAM,Cam]=comstr(CAM,1);
  else Unit=''; end
  if ischar(CAM); [i1,CAM,Cam]=comstr(CAM,'dbval','%i');else; i1=[];end
  if isempty(CAM)&carg<=nargin st=varargin{carg};carg=carg+1; 
  else st=CAM;end
  if isempty(st)
  elseif ischar(st); mat=m_elastic('database',st); 
  elseif isnumeric(st)
    [typ,st1,i4]=fe_mat('typem',st(2));
    mat=struct('pl',st,'name',sprintf('%i',st(1)),'type',typ,'unit',st1);
  end
  if ~isempty(Unit)
   mat.pl=fe_mat(sprintf('convert %s %s',mat.unit,Unit),mat.pl);
  end
  if (length(i1)==1); mat.pl(1)=i1;end
  r1=mat.pl; 
  if ~isempty(pl) ; i2=find(pl(:,1)==r1(1)); else i2=[];end
  if isempty(i2)  ; i2=size(pl,1)+1;end
  pl(i2,1:length(r1))=r1; 
  out1(end+1,1:3)={'mat',sprintf('%i_%s',mat.pl(1),mat.name),mat};
  if carg>nargin break;end
  CAM=varargin{carg};carg=carg+1;
 end
 out=pl;

% -------------------------------------------------------------------------
elseif comstr(Cam,'database'); [st,Cam]=comstr(CAM,9);

  if isempty(st)&carg<=nargin st=varargin{carg}; carg=carg+1;end
  [MatId,i2,i3,i4]=sscanf(st,'%i',1); if i2~=1 MatId=1;end
  st=comstr(st,i4);

  out.pl=[MatId fe_mat('type','m_elastic','SI',1) ... 
          210e9 .285 7800 210e9/2/(1.285)];
  out.name='Steel';
  out.type='m_elastic';
  out.unit='SI';

  out(2).pl=[MatId fe_mat('type','m_elastic','SI',1) 72e9 .3 2700 72e9/2/(1.3)];
  out(2).name='Aluminum';
  out(2).type='m_elastic';
  out(2).unit='SI';

  out(3).pl=[MatId fe_mat('type','m_elastic','SI',1) 2.4e9 .3  1.8e3 2.4e9/2/1.3];
  out(3).name='Nylon';
  out(3).type='m_elastic';
  out(3).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',1) ...
    120e9 .33 4700 120e9/2/(1.33)];
  out(end).name='Titanium';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',2) 1.225 330 0];
  out(end).name='Air';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',2) 1000 1500 0];
  out(end).name='Water';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',1) ...
    4e9 .4 1200 4e9/2/(1.4) 0 1e-4 20];
  out(end).name='Epoxy_Resin';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',1) ...
    3e9 .4 1200 3e9/2/(1.4) 0 8.5e-5 20];
  out(end).name='Polyamide_Resin';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',1) ...
          4.e10 .17 2200 4.e10/2/(1.17) .01 1.2e-5 20];
  out(end).name='Concrete';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',1) ...
          15e9 .3 1380 15e9/2/(0.3)]; 
  % found in "techniques de l'ingénieur", grande incertitude pour Nu
  out(end).name='Nomex';
  out(end).type='m_elastic';
  out(end).unit='SI';

  out(end+1).pl=[MatId fe_mat('type','m_elastic','SI',1) 1 0.0 1 0.5]; 
  out(end).name='Strain';
  out(end).type='m_elastic';
  out(end).unit='SI';


  i1=strmatch(comstr(st,-27),comstr({out.name},-27),'exact');
  out1='Elastic';

  if isempty(i1) & isempty(st); return; end

  % match a name 
  if ~isempty(i1) out=out(i1);
  elseif comstr(st,'lamina')

% From Jones p.91 Estimate composite lamina properties
   
 opt=comstr(st(7:end),-1);v=[opt(1) 1-opt(1)];
 if opt(5)==0; opt(5)=opt(2)/2/(1+opt(3));end
 if opt(9)==0; opt(9)=opt(6)/2/(1+opt(7));end
 out=struct('pl',[MatId fe_mat('type','m_elastic','US',5) ...
   opt(2)*v(1)+opt(6)*v(2) ... % E_1=E_f V_f + E_m V_m (3.7)
   opt(2)*opt(6)/(v(2)*opt(2)+v(1)*opt(6)) ... % E_2=E_fE_m/V_mE_f+V_fE_m 3.13
   opt(3)*v(1)+opt(7)*v(2) ... % \nu_{12}=V_m\nu_m+V_f\nu_f
   opt(5)*opt(9)/(v(2)*opt(5)+v(1)*opt(9)) ...% G_12=G_fG_m/V_mG_f+G_fG_m 3.27
   0 0 opt(4)*v(1)+opt(8)*v(2)  ... % rho
   ],'unit','US','type','m_elastic','name',sprintf('lamina_%i',MatId));

  else % assume values given
    error('Not a supported material');
  end


% -------------------------------------------------------------------------
elseif comstr(Cam,'buildconstit');[out,out1,out2]=p_solid(varargin);
% -------------------------------------------------------------------------
%[S,RhoH]=p_shell('buildply',mat,pl,il,cz,rhoh);
elseif comstr(Cam,'buildply')

mat=varargin{carg};carg=carg+1;
pl=varargin{carg};carg=carg+1;
il=varargin{carg};carg=carg+1;
cz=varargin{carg};carg=carg+1;
rhoh=varargin{carg};carg=carg+1;

[st,unit,typ]=fe_mat('type',mat(2));
S=zeros(5);
   if strcmp(st,'m_elastic')&typ==5  
     % Orthotropic material for shell. See    help m_elastic
     % ply compliance matrix
     %S=[1/E1(i) -Nu12(i)/E1(i) 0 0 0;-Nu12(i)/E1(i) 1/E2(i) 0 0 0;
     %    0 0  1/G12(i) 0 0;0 0 0 1/G23(i) 0;0 0 0 0 1/G13(i)];
     %S=eye(5);S(1:2,1:2)=1;sprintf('%i ',find(S))
     S([1 2 6 7 13])= [1/mat(3) -mat(5)/mat(3) -mat(5)/mat(3) 1/mat(4)  1/mat(6)];
     if mat(7); S(19)=1/mat(7); else; S(19)=5/6/mat(6); end; % xxx
     if mat(8); S(25)=1/mat(8); else; S(25)=5/6/mat(6); end
     rhoh=rhoh+diff(cz)*mat(9);
   elseif strcmp(st,'m_elastic')&typ==1 % isotropic material
     % defaults for elastic istropic material
     if strcmp(st,'m_elastic')&typ==1 % Standard isotropic
      if length(mat)<6|mat(6)==0 mat(6)=mat(3)/2/(1+mat(4));end % G
      if size(mat,2)<7;mat(7)=0;end
     end
     %S=[1/E -Nu/E 0 0 0;-Nu/E1 1/E 0 0 0;
     %    0 0  1/G 0 0;0 0 0 1/G 0;0 0 0 0 1/G];
     S([1 2 6 7 13])= [1/mat(3) -mat(4)/mat(3) -mat(4)/mat(3) 1/mat(3) ...
      1/mat(6)];S(19)=5/6/mat(6);S(25)=5/6/mat(6); % xxx
     rhoh=rhoh+diff(cz)*mat(5);
   elseif any(exist(st)==[2 3 6])
     [S,rhoh]=feval(st,'buildply',mat,pl,il0);
   else; error('Ply not implemented for this type of material');
   end 
   out=S; out1=rhoh;out2=[];

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

 if carg<2 error('bad number of input arguments'); end
 %out=PropertyUnitType_elastic(varargin{2});
MaterialSubType=varargin{carg};
% the indices match rows in
% fe_mat('convertSIGM');[num2cell(1:size(ans,1))' ans(:,1:2)]
out1={};

switch MaterialSubType % fe_mat('unitlabel','SI')
case 1 % isotropic   [MatID typ E nu rho G eta alpha T0]
               ind = [0     0   1 0  3   1 0   8     7]; 
 st=...
 {'MatId'                  0;
  'Type'                   0;
  'E : Youngs Modulus'     1;
  'Poisson''s ratio'       0;
  'Density'                3;
  'Shear Modulus'          1;
  'Loss factor'            0;
  'Thermal expansion coef' 8;
  'Reference temperature'  7 };

case 2 % acoustic fluid [MatId typ rho C eta]
               ind =    [0     0   3   5 0];         
 st=...
 {'MatId'          0
  'Type'           0
  'Rho (density)'  3
  'C (velocity)'   5
  'Loss factor'    0};
case 3 % 3-D anisotropic solid [MatId typ Gij           rho A1... A6 T0 eta]
               ind =           [0     0    1*ones(1,21) 3  8*ones(1,6) 7 0];    
 st=...
 {'MatId' 0;
  'Type'  0;
  'G11'   1;
  'G12'   1;
  'G22'   1;
  'G13'   1;
  'G23'   1;
  'G33'   1;
  'G14'   1;
  'G24'   1;
  'G34'   1;
  'G44'   1;
  'G15'   1;
  'G25'   1;
  'G35'   1;
  'G45'   1;
  'G55'   1;
  'G16'   1;
  'G26'   1;
  'G36'   1;
  'G46'   1;
  'G56'   1;
  'G66'   1;
  'Rho (density)' 3;
  'A1'    8;
  'A2'    8;
  'A3'    8;
  'A4'    8;
  'A5'    8;
  'A6'    8;
  'Reference temperature' 8;
  'Loss factor'           0}; 

case 4 % 2-D anisotropic solid 
       % [MatId typ E11 E12 E22 E13 E23 E33 rho eta a1 a2 a3]
   ind = [0     0   1*ones(1,6)             3   0   8*ones(1,3)]; 
 st=...
 {'MatId' 0;
  'Type'  0;
  'E11'   1;
  'E12'   1;
  'E22'   1;
  'E13'   1;
  'E23'   1;
  'E33'   1;
  'Density' 3;
  'Loss factor' 0;
  'Alpha1'  8;
  'Alpha2'  8;
  'Alpha3'  8 };

case 5 % Orthotropic material for shell
 %[MatId type E1 E2 nu12 G12 G1z G2z Rho A1 A2 TREF Xt Xc Yt Yc S Ge F12 STRN]
 %                                                  xxx not valid xxx
 ind=[0 0     1 1   0    1   1   1   3   8  8  7    0  0  0  0  0 0  0   0];
 st=...
 {'MatId' 0;
  'Type'  0;
  'E11'   1;
  'E12'   1 };
otherwise; st={'MatId' 0; 'Type', 0};
end

if ~isempty(strfind(Cam,'cell')); out=st; else; out=[st{:,2}]; end

% -------------------------------------------------------------------------
elseif comstr(Cam,'subtype');[CAM,Cam]=comstr(CAM,8);

st={'1 (Isotropic)','2 (Acoustic fluid)',...
    '3 (3-D anisotropic solid)','4 (2-D anisotropic solid)'}; 
if carg<=nargin;
 i1=varargin{carg}; carg=carg+1;
 if ischar(i1)
  out=strmatch(i1,st);
 else;
  try; out=st{i1};
  catch; out=sprintf('m_elastic %i',i1);
  end
 end
end

% -------------------------------------------------------------------------
else; sdtw('''%s'' not known',CAM);
end % commands

