function [o1,o2,o3,o4,o5]=fe_mat(varargin)

%FE_MAT general material/section property handling routine
%
%       Syntax:   out = fe_mat('unit')
%                 out = fe_mat('unitlabel',UnitSystemCode)
%                 typ=fe_mat('type','m_function',UnitCode,SubType)
%                 [m_function',UnitCode,SubType]=fe_mat('type',typ)
%                 [o1,o2,o3]=fe_mat(ElemP,ID,pl,il)
%                 out = fe_mat('convert si ba',pl);
%
%	This function is used by a number of element functions to support 
%	material property declarations. The input arguments are the
%	element function name or parent name (ElemP), ID =  [MatID ProID]
%	contains identification numbers for the current material/property,
%	PL and IL are the material/element property matrices.
%
%	Outputs depend on the element type and are detailed in the manual.
%
%	The general form for material property rows is
%	   [MatID Type Prop]
%	  with
%	   MatID a positive integer identifying a particular material property
%	   Type  identifying the type of material. Currently supported types
%	         1 standard isotropic
%	           no other type currently supported
%	   Prop  as many properties as needed
%
%	See also help m_elastic, fe_mk, beam1, bar1, tria3, ...
%		 doc  pl, m_elastic, elem0, eltfun

%	Etienne Balmes   16/01/96, 23/11/01
%       Copyright (c) 2001-2002 by INRIA and SDTools
%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       All Rights Reserved.

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

if nargin==0 help fe_mat;return; end
% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
if ischar(varargin{1})

 [CAM,Cam]=comstr(varargin{1},1); carg=2;

if comstr(Cam,'p_') Cam='type'; carg=1; end
if comstr(Cam,'m_') Cam='type'; carg=1; end

if comstr(Cam,'steel')    o1 =m_elastic('database','steel');
elseif comstr(Cam,'alum') o1 =m_elastic('database','aluminum');
elseif comstr(Cam,'types') % search available material types

 o1=matgui('search m_*');

% GetPl - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Create the pl field no matter initial format
% pl=fe_mat('getpl',model)
elseif comstr(Cam,'getpl')  [CAM,Cam]=comstr(CAM,6);

 model=varargin{carg};carg=carg+1;
 val=stack_get(model,'mat');
 if isfield(model,'pl') pl=model.pl;else pl=[];end
 for j1=1:size(val,1)
   r1=val{j1,3}; if isfield(r1,'pl')
    r2=r1.pl(:)'; if isempty(pl) i2=[];else i2=find(pl(:,1)==r2(1));end 
    if isempty(i2) pl(end+1,1:length(r2))=r2;
    elseif length(r2)==1
    else pl(i2,1:length(r2))=r2;
    end
   end
 end
 o1=pl;
 if comstr(Cam,'of_mk') error('You should use GetPlIlOfMk');end
 % for j1=1:size(o1,1)  r1=fe_mat('of_mk',o1(j1,:)); 
 %  o1(j1,:)=0;o1(j1,1:length(r1))=r1;
 % end
 %end

% GetIl - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Create the il field no matter initial format
% il=fe_mat('getil',model)
elseif comstr(Cam,'getil')  [CAM,Cam]=comstr(CAM,6);

 model=varargin{carg};carg=carg+1;
 val=stack_get(model,'pro');
 if isfield(model,'il') il=model.il;else il=[];end
 for j1=1:size(val,1)
   r1=val{j1,3}; if isfield(r1,'il')
    r2=r1.il(:)'; if isempty(il) i2=[];else i2=find(il(:,1)==r2(1));end 
    if isempty(i2) il(end+1,1:length(r2))=r2;
    elseif length(r2)==1
    else il(i2,1:length(r2))=r2;
    end
   end
 end
 o1=il;

% convert - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% convert from a specified unit system to another
% pl1=fe_mat('convert',pl,PropertyUnitTypeVector);
elseif comstr(Cam,'convert')  [CAM,Cam]=comstr(CAM,8);

  if isempty(Cam)&nargin==1
   r1=fe_mat('unitsystem');
   st=cellstr(num2str([1:length(r1)]'));
   st(:,2)={r1.name};st=st';
   fprintf('%s  %s\n',st{:})
   return;
  end
  [r1,lab]=fe_mat('unitsystems');
  if carg>nargin pl=[]; else pl=varargin{carg};carg=carg+1; end

  if carg<=nargin
    ind=varargin{carg};carg=carg+1; %PropertyUnitType
    [m_function,UnitCode,SubType]=fe_mat('type',pl(2)); 
  elseif isempty(pl)  % no argument
    ind=1:size(lab,1);ind(10)=0;ind=find(ind);
    UnitCode=0;
  elseif size(pl,1)>1  % Multiple conversion
    for j1=1:size(pl,1)
      r1=fe_mat(['convert' CAM],pl(j1,:));
      pl(j1,1:length(r1))=r1;
    end
    o1=pl;
    return; 
  else 
   [m_function,UnitCode,SubType]=fe_mat('type',pl(2)); 
   if ~exist(m_function,'file') m_function(1)='p';end
   eval(sprintf('ind= %s(''PropertyUnitType'',%i);',m_function,SubType));
  end

  if ~isempty(pl) % get unit types
    [st,i1,isub]=fe_mat('type',pl(2)); 
    if length(pl)<length(ind) pl(length(ind))=0; end
  end

  % get initial and final unit systems
  if comstr(upper(Cam(1:2)),'US')
   [CAM,Cam]=comstr(CAM,3);
   if comstr(upper(Cam(1:2)),'US')
    o1={1,1,1,1,1,1,1,1,1,1,1,1,1}';
    o1(:,end+1)=lab(:,end-1);
    return;
   end
   error('bad unit system. Can''t convert');
  end
  if length(CAM)==2 i1=UnitCode;
  else i1=strmatch(upper(Cam(1:2)),{r1.name});[CAM,Cam]=comstr(CAM,3);
  end
  i3=strmatch(upper(Cam(1:2)),{r1.name});
  if isempty(i3) | isempty(i1) error('bad unit system. Can''t convert');end
  if exist('UnitCode','var') & UnitCode & i1~=UnitCode 
   %warning('There is a mismatch in the initial unit system');
  end


  % length, force, temp, temp-offset, time
  r2=r1(i3).data(3:7)./r1(i1).data(3:7);           % basic unit conversion
  r3=reshape([lab{:,end}],length(r2),size(lab,1))';% rows in basic units

  if isempty(ind)
  else
    ind(find(ind==0))=10; % code 10 for no unit 
    in2=round(rem(ind,1)*1000);in2=find(in2);
    if ~isempty(in2); r4=r3(in2,:);end % denominator units
    r3=r3(fix(ind),:); % numerator units 
    r3=prod((r2(ones(length(ind),1),:)).^r3,2)';
    if ~isempty(in2)
     r4=prod((r2(ones(length(in2),1),:)).^r4,2)';
     r3(in2)=r3(in2)./r4;
    end

    if ~isempty(pl) % update values & type Identifier
     o1(1:length(ind))=pl(1:length(ind)).*r3;
     o1(2)=fe_mat('type',m_function,i3,isub);
    else
     pl=ones(size(ind));pl(2,:)=r3;
     o1=num2cell(pl');o1(:,3:4)=lab(ind,[i1 i3]);o1=o1(:,[1 3 2 4]);
   end
  end

% Units handling - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'unit')  [CAM,Cam]=comstr(CAM,5);

if comstr(Cam,'label') % 'UnitLabel' - - - - - - - - - - - - - - - - - -

o1={'Pa','lbf/ft^2','kgf/m^2','pdl/ft^2','milli-N/mm^2',...
                         'centi-N/cm^2','lbf/in^2(psi)','kgf/m^2','N/mm^2'};
o1(2,:)={'N','lbf','kgf','pdl','milli-N','centi-N','lbf','kgf','N'};
o1(3,:)={'kg/m^3','lbf-s^2/ft^4','kgf-s^2/m^4','lbm/ft^3',...
          'kg/mm^3','kg/cm^3','lbf-s^2/in^4','kgf-s^2/m^4','t/mm^3'};
o1(4,:)={'m','ft','m','ft','mm','cm','in','mm','mm'};
o1(5,:)={'m/s','ft/s','m/s','ft/s','mm/s','cm/s','in/s','mm/s','mm/s'};
o1(6,:)={'m/s^2','ft/s^2','m/s^2','ft/s^2','mm/s^2','cm/s^2','in/s^2','mm/s^2','mm/s^2'};
o1(7,:)={'^oC','^oF','^oC','^oF','^oC','^oC','^oF','^oC','^oC'};
o1(8,:)={'/^oC','/oF','/oC','/oF','/oC','/oC','/oF','/oC','/oC'};
o1(9,:)={'kg','lbf-sec^2/ft','kgf-sec^2/m','lbm','kg','kg',...
                 'lbf-sec^2/in','kgf-sec^2/mm','t'};
o1(10,:)={' ',' ',' ',' ',' ',' ',' ',' ',' '};
o1(11,:)={'m^4','ft^4','m^4','ft^4','mm^4','cm^4','in^4','mm^4','mm^4'};
o1(12,:)={'m^2','ft^2','m^2','ft^2','mm^2','cm^2','in^2','mm^2','mm^2'};
o1(13,:)={'N/m','lbf/ft','kgf/m','pdl/ft','milli-N/mm',...
          'centi-N/cm','lbf/in','kgf/mm','N/mm'};

o1(:,end+1)={'pressure','force','density','length','speed','acceleration',...
         'temperature','thermal coeff','mass','no unit','inertia','area',...
         'stiffness'}';

% length, force, temp, temp-offset, time
o1(:,end+1)={[-2 1 0 0 0],[0 1 0 0 0],[-4 1 0 0 2],[1 0 0 0 0],[1 0 0 0 -1],...
          [1 0 0 0 -2],[0 0 1 0 0],[0 0 -1 0 0],[-1 1 0 0 2],[0 0 0 0 0],...
          [4 0 0 0  0],[2 0 0 0 0],[-1 1 0 0 0]}';

if nargout==0; 
 o1(:,end+1)=num2cell([1:size(o1,1)]');disp(o1(:,[end end-2 1 2]));
 clear o1
end
% System of Units        Length  Time      Mass          Force       Temp(R:A)
%----------------------------------------------------------------------------- 
% SI-Meter (newton)      meter   sec   kilogram(kg)       newton(N)  deg C:K
% BG-Foot (pound f)      foot    sec   lbf-sec**2/ft     pound(lbf)  deg F:R
% MG-Meter (kilogram f)  meter   sec   kgf-sec**2/m   kilogram(kgf)  deg C:K
% BA-Foot (poundal)      foot    sec   pound(lbm)      poundal(pdl)  deg F:R
% MM-mm (milli-newton)   mm      sec   kilogram(kg)    milli-newton  deg C:K
% CM-cm (centi-newton)   cm      sec   kilogram(kg)    centi-newton  deg C:K
% IN-Inch (pound f)      inch    sec   lbf-sec**2/in     pound(lbf)  deg F:R
% GM-mm (kilogram f)     mm      sec   kgf-sec**2/mm  kilogram(kgf)  deg C:K
% TM-mm (ton- mm)        mm      sec   ton               N           deg C:K
% US-User defined        ____    sec   _____________  _____________  _____

  carg=2;if carg<=nargin
   st=varargin{carg};carg=carg+1;
   if isfield(st,'unit') st=st.unit;end
   if ~ischar(st) error('unit must be a string');end
   i1=strmatch(st,{'SI','BG','MG','BA','MM','CM','IN','GM','TM'});
   if ~isempty(i1)
    o1=o1(:,i1);
   end
  end

elseif comstr(Cam,'system')  % 'UnitSystem' - - - - - - - - - - - - - - - - - -
  [CAM,Cam]=comstr(CAM,7);
  % Code TempCode (1 absolu, 2 relative) length, force, temp, temp-offset, time
  % to convert to SI units (divide by factor)   
  o1=struct('name','SI Meter (newton)', ...
           'data',[1 2 1 1 1 1 1]);

  o1(2)=struct('name','BG Foot (pound f)', ...
           'data',[2 2 3.28083989501312350e+00  2.24808943099710480e-01 ...
                   1.80000000000000000          4.59670000000000002e+02 1]);

  o1(3)=struct('name','MG Meter (kilogram f)', ...
           'data',[3 2 1 1.01971621297792830e-01 1 2.73149999999999980e+02 1]);

  o1(4)=struct('name','BA Foot (poundal)', ...
           'data',[4 2 3.28083989501312350  7.23301385120989430  ...
                   1.80000000000000000   4.59670000000000020e+02 1]);

  o1(5)=struct('name','MM mm (milli-newton)', ...
           'data',[5 2 1.00000000000000000e+03  1.00000000000000000e+03  ...
                   1.00000000000000000  2.73149999999999980e+02 1]);

  o1(6)=struct('name','CM cm (centi-newton)', ...
           'data',[6 2 1.00000000000000000e+02  1.00000000000000000e+02  ...
                   1.00000000000000000    2.73149999999999980e+02 1]);

  o1(7)=struct('name','IN Inch (pound f)', ...
           'data',[7 2 3.93700787401574810e+01  2.24808943099710470e-01  ...
                   1.80000000000000000     4.59670000000000020e+02 1]);

  o1(8)=struct('name','GM mm (kilogram f)', ...
           'data',[8 2 1.00000000000000000e+03  1.01971621297792830e-01  ...
                   1.00000000000000000       2.73149999999999980e+02 1]);

  o1(9)=struct('name','TM mm (Newton)', ...
           'data',[9 2 1.00000000000000000e+03  1.00000000000000000e+00  ...
                   1.00000000000000000  2.73149999999999980e+02 1]);

  if nargout==2; o2=fe_mat('unitlabel'); end
  if comstr(Cam,'name'); o1={o1.name}; end

else sdtw('''Unit%s'' not know',Cam);
end

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

  [CAM,Cam]=comstr(CAM,5); if isempty(Cam) Cam='m'; end

  st1={'SI','BG','MG','BA','MM','CM','IN','GM','US'};
  r1=varargin{carg};carg=carg+1;
  if isnumeric(r1)
   % type to result
   if Cam=='p'
    if     r1==1 o1='p_beam';o2=1;o3=1;
    elseif r1==2 o1='p_shell';o2=1;o3=1;
    elseif r1<0|r1> 4.503599627370496e+001
     o1='p_null';o2=0;o3=0;
    else
     try
      st=dec2base(round(1e14*r1),36);
      if length(st)>2
       o2=comstr(st(end-1),[-1 1]); o3=comstr(st(end),[-1 1]);
       o1=['p_' comstr(st(1:end-2),-27)];
      else o2=0; o3=0; o1='';
      end
     catch
      o1='p_null';o2=0;o3=0;
     end
    end
   else
    if r1==1 o1='m_elastic';o2=1;o3=1;
    elseif r1<0|r1> 4.503599627370496e+001
     o1='m_null';o2=0;o3=0;
    else
     try
      st=dec2base(round(1e14*r1),36);
      if length(st)>2
       o2=comstr(st(end-1),[-1 1]); o3=comstr(st(end),[-1 1]);
       o1=['m_' comstr(st(1:end-2),-27)];
      else o2=0; o3=0; o1='';
      end
     catch
      o1='m_null';o2=0;o3=0;
     end
    end
   end % m or p
   if nargout==4 
    if any(o2==[1:9]) o2=st1{o2}; else o2='US';end;
    o4=o3; try; o3=feval(o1,'subtypestring',o3); catch; o3=o1; end
   end
  else % code properties into type number
   if ~comstr(r1,'m_') & ~comstr(r1,'p_')
    error('Material functions should start with ''m_''');
   end
   r1=r1(3:end);
   if carg<=nargin st=varargin{carg};carg=carg+1;else st='1'; end
   if isnumeric(st) st=num2str(st);
   elseif any(strmatch(upper(st),st1)) st=num2str(strmatch(upper(st),st1));
   end
   if length(st)>1 error('Unit code is between 0 and 9');end
   r1(end+1)=st;
   if carg<=nargin st=varargin{carg};carg=carg+1;else st='1'; end
   if isnumeric(st) st=num2str(st);end
   if length(st)>1 error('Sub type code is between 0 and 9');end
   r1(end+1)=st;
   o1=base2dec(upper(r1),36)*1e-14;
  end

% GetPlIl - - - - - - - - - - - - -
elseif comstr(Cam,'plil') [CAM,Cam]=comstr(CAM,5);% - - - - - - - - - - - - -

ID=varargin{2};pl=varargin{3};il=varargin{4};
st='';
if ~isempty(pl)
 i1=find(pl(:,1)==ID(1));
 if isempty(i1)&ID(1)~=0 st=sprintf('MatID %i not found PL',ID(1));pl=[];
 elseif ID(1)==0 pl=[]; else pl=pl(i1,:); end
end

if ~isempty(il)
 i2=find(il(:,1)==ID(2));
 if isempty(i2)&ID(2)~=0 
   st=sprintf('%s, ProID %i not found IL',st,ID(2));il=[ID(2) 0 1];
 elseif ID(2)==0 il=[0]; else il=il(i2,:); end
else
 if ID(2)==0 ID(2)=ID(1); end
 il=[ID(2) 0 1];
end

o3=[];
if ~isempty(pl)&comstr(Cam,'of_mk')
  [o1,o2]=fe_mat(Cam,pl,il); % [Constit,iopt]
  if size(o2,1)==1&length(o2)>1; o2=o2(:);end
  o2(end+1,1)=-9999; % This is a flag to use Modulef Elements
else
  if ~isempty(st); fprintf('%s\n',st);end
  o1=pl; o2=il;  if ~isempty(pl) pl=pl(3:end); end  
end

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

% this command reformats any SDT Material entry to the row input needed by
% of_mk calls to element functions
%
% type 1 : plane stress isotropic
%      2 : plane strain isotropic or 3-D isotropic
%      3 : 3-D anisotropic
%      4 : 2-D anisotropic
%      12 : reserved for acoustic fluid
%
% These are not coherent with SDT reference defined in m_elastic

pl=varargin{carg};carg=carg+1;
il=varargin{carg};carg=carg+1;
[CAM,Cam]=comstr(Cam,6);

o2=[pl(1) il(1) zeros(1,6)];
if comstr(Cam,'2d') % - - - - - - - - - - - - - - - - - - - - - - - - - -

  % INTEG Formatting for 2D elements by Amine
  %[MatId ProId NDOF NNode Type 0 Formulation IoptStiff IoptMass N]
  % of_mk_subs.f compute_matrix IOPT indices
  %              1     2    3   4  5           6        7        8 
  %
  % Values of Formulation  :
  %           IoptStiff    : always 3 (anisotropic call)
  %           IoptMass     : currently unused
  %                  N     : only used by q9a eraq2c.f element   

  if ~any(il(3)==[0 1 2]) 
    error(sprintf('%i: unknown 2-D element formulation',il(3)));
  end
  o2(7)=il(3)+1; % Formulation (plane or axisymetric)
  if length(il)>3 o2(8)=il(4);end % NFourier

elseif comstr(Cam,'3d') % - - - - - - - - - - - - - - - - - - - - - - - - - -

  % INTEG formatting for 3D elements by Marina Vidrascu
  % [MatId ProId NDOF NNode Type 0 Formulation Iopt]
  % of_mk_subs.f compute_matrix IOPT indices
  %              1     2    3   4  5           6        
  % IoptStiff  2: isotropic, 1: orthotropic

else % other % - - - - - - - - - - - - - - - - - - - - - - - - - -

end

[st,UnitCode,SubType]=fe_mat('type',pl(2));

switch st
case 'm_elastic'
 
 switch SubType
 case 1 % standard isotropic of_mk 

  if size(pl,2)<7 pl(7)=0;end

  if comstr(Cam,'2d') % For 2D only : constit=[rho eta E11 ... E33 a1 a2 a3 T0]

     %if o2(8)==0 o2(8)=o2(7);end 
     o1 = pl([5 7]); % [rho eta]
     E =pl(3); nu=pl(4); 

     switch il(3) % Formulations
     case 1 % - - - - - - - - - -  - - - - - - - plane stress
       o2(8)=3; % anisotropic element is always called
       C=E/(1.-nu^2);r1=[C C*nu C 0. 0. C*(1-nu)/2];
       o1(3:8)=r1;
       o1(9:12)=0; % a1 a2 a3 T0
     case 0 % - - - - - - - - - -  - - - - - - - plane strain
       o2(8)=3; % anisotropic element is always called
       C=E*(1-nu)/(1+nu)/(1-2*nu);
       r1=[C C*nu/(1-nu) C 0. 0. C*(1-2*nu)/2/(1-nu)];
       o1(3:8)=r1;  o1(9:12)=0; % a1 a2 a3 T0
     case 2 % - - - - - - - - - -  - - - - - - - axisymetric
 
       o1 = pl([5 7 3 4]);% xxx only isotropic case => o2(8)=0;
       o1(9:12)=0; % a1 a2 a3 T0

     end % Formulations 


  elseif  pl(6)==0 | abs(pl(6)-pl(3)/2/(1+pl(4)))<sqrt(eps)*pl(3) 

     o1 = pl([5 7 3 4]); o2(8)=2; % 3D isotropic : constit=[rho eta E nu] 

  else % G is not E/2/(1+nu) use anisotropic
     r1=pl(6);r2=pl(3)*pl(4)/(1+pl(4))/(1-2*pl(4));
     o1=[pl([5 7]) r2+[2*r1 0 2*r1 0 0 2*r1] r1 r1 r1];
     o2(8)=1; % 3D anisotropic
  end

 case 4 % 2-D Anisotropic : car = [rho eta E11 ... E33 a1 a2 a3 T0]

    if length(pl)<9 pl(14)=0.; end 
    o1=pl([9 14 3:8 10:13]);
    o2(8)=3; % anisotropic

 case 3 
 % 3-D Anisotropic : car = [rho eta G11 G12 G22 G13 G23 G33 G14 ...  G66]
 %  xxx this is not consistent with m_elastic documentation xxx
   if any(pl(2+[7 8 9 11:14 16:20]))
    % il(3) calls the old .m file assembly with fully anisotropic material
    % Constit has [ rho eta dd(:)'] 
    i1=find(triu(ones(6))); dd=zeros(6);dd(i1)=pl(3:23);
    dd=dd+tril(dd',-1); if length(pl)<25 pl(25)=0; end
    o1=[pl(24:25) dd(:)'];
    o2(8)=4; % anisotropic
   else % actually orthotropic [... eta]
    pl=[pl(1) 3 pl([24 3:8 12 17 23 40])];
    o2(8)=3;        % orthotropic
   end

 case 2 % fluid (no of_mk)

 end

otherwise error('Material function not supported by of_mk');
end

else error('Not a valid fe_mat command'); end


% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
% --------------------------------------------------------------------------
% THIS SECTION IS OBSOLETE AND WILL EVENTUALLY DISAPPEAR IN P_* functions
else % ElemF is not a string

if nargin<4 error('Not the right number of inputs');end
ElemF=varargin{1};ID=varargin{2};pl=varargin{3};il=varargin{4};
if nargin>4 t=varargin{5}; end

% find the good properties

if ~isempty(pl)
 i1=find(pl(:,1)==ID(1));
 if isempty(i1)&ID(1)~=0 warning(sprintf('MatID %i not found PL',ID(1)));pl=[];
 elseif ID(1)==0 pl=[]; else pl=pl(i1,:); end
end

if ~isempty(il)
 i2=find(il(:,1)==ID(2));
 if isempty(i2)&ID(2)~=0 warning(sprintf('ProID %i not found IL',ID(1)));il=[];
 elseif ID(2)==0 il=[]; else il=il(i2,:); end
end

if size(pl,1)>1 
  pl=pl(1,:);fprintf(1,'taking first property with MatID %i\n',pl(1));
end
if size(il,1)>1 il=il(1,:);
  fprintf(1,'taking first property with ProID %i\n',il(1));
end

if ~isempty(pl) [st,i1,i2]=fe_mat('typem',pl(2));end
if isempty(pl) 
elseif pl(2)==1|i2==1 % standard isotropic material

   if length(pl)==5 	pl=[pl 0]; 		 end
   if pl(6)==0 		pl(6)=pl(3)/2/(1+pl(4)); end

else 
 error(sprintf('MatId %i : %s,%i,%i not supported for beams',pl(1),st,i1,i2));
end

% element specific treatment of element properties

if ElemF==1 %comstr(ElemF,'beam1')

   %if il(2)~=1 error('Not a supported beam1 property'); end
   if length(il)<9 	il(9)=0; end
   o1=pl(3:length(pl)); o2=il(3:9);

elseif ElemF==2 % comstr(ElemF,'shell1')

   %  1    2   3      4   5   6      7     8     9      10   11 12 13  14
   % ProID 2 Formul Drill 0 Thick k(TS/T) MID2 12l/T^3 MID3 NSM Z1 Z2 MID4

   if isempty(pl) error('No declared material property');end
   il(15)=0;
   if il(2)==2
   else [st,i1,i2]=fe_mat('typep',il(2)); 
    if ~strcmp(st,'p_shell')
     error(sprintf('ProId %i is not a shell property (%s)',il(1),st));
    end
    if i2~=1 sdtw('p_shell subtype %i is supported by q4cs elements',i2);end 
   end
   if nargin>4 & ~isempty(t) il(6)=t; end
   if il(4)==0 il(4)==1; end % drilling DOF stiffness
   o1 = pl; o2 = il;
   o3 = pl(3)/(1-pl(4)^2); % membrane
   o3 = [o3 pl(4)*o3 0;      pl(4)*o3 o3 0;      0 0 (1-pl(4))/2*o3]*il(6);
   o4 = o3*il(6)^2/12;     % bending
   if il(9)~=0 o4=o4*il(9); end  % ratio for bending moment of inertia
   if any(il([8 10 12 13 14]))
     %sdtw('_nb','Unsupported properties in shell declaration');
   end
   if length(il)<7 il(7)=5/6; elseif il(7)==0 il(7)=5/6;end
   o5 = eye(2)*pl(6)*il(7)*il(6); % shear

%a=0; % a=0 --> Deformations planes
%     % a=1 --> Contraintes planes

%H1 = pl(1,1)*(1-a*pl(1,2))/((1+pl(1,2))*(1-pl(1,2)-a*pl(1,2)));
%Hm = H1*[1 pl(1,2)/(1-a*pl(1,2)) 0;
%         pl(1,2)/(1-a*pl(1,2)) 1 0;
%         0 0 (1-pl(1,2)-a*pl(1,2))/(2*(1-a*pl(1,2)))];

elseif ElemF==3 %comstr(ElemF,'solid1')

    o1 = pl; o2 = il;

    if pl(6)==0 pl(6)=(pl(3)/2/(1+pl(4))); end
    r1=pl(6);r2=pl(3)*pl(4)/(1+pl(4))/(1-2*pl(4));
    dd=zeros(6);
    dd([1 2 3 7 8 9 13 14 15])=r2+[2*r1 0 0  0 2*r1 0  0 0 2*r1];
    dd([22 29 36])=pl(6);
    o3=dd;

else error('Not a supported material/property call');
end

% THIS SECTION IS OBSOLETE AND WILL EVENTUALLY DISAPPEAR IN P_* functions
end % ElemF is not a string



%============================================================
%Material               : 1 - GENERIC_ISOTROPIC_STEEL
%============================================================
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : PASCAL
%  Scalar Value :  2.068000E+11 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : KILOGRAM/METER^3
%  Scalar Value :  7.820000E+03 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : KELVIN
%  Scalar Value :  2.950000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : POUND FORCE/FOOT^2
%  Scalar Value :  4.319109E+09 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : LBF SEC^2/FT^4
%  Scalar Value :  1.517332E+01 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : RANKINE
%  Scalar Value :  5.310000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : KILOGRAM FORCE/M^2
%  Scalar Value :  2.108773E+10 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : KGF SEC^2/M^4
%  Scalar Value :  7.974181E+02 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : KELVIN
%  Scalar Value :  2.950000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : POUNDAL/FOOT^2
%  Scalar Value :  1.389632E+11 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : POUND MASS/FOOT^3
%  Scalar Value :  4.881870E+02 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : RANKINE
%  Scalar Value :  5.310000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : MILLINEWTON/MILLIMETER^2
%  Scalar Value :  2.068000E+08 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : KILOGRAM/MILLIMETER^3
%  Scalar Value :  7.820000E-06 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : KELVIN
%  Scalar Value :  2.950000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : CENTINEWTON/CENTIMETER^2
%  Scalar Value :  2.068000E+09 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : KILOGRAM/CENTIMETER^3
%  Scalar Value :  7.820000E-03 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : KELVIN
%  Scalar Value :  2.950000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : POUND FORCE/INCH^2
%  Scalar Value :  2.999380E+07 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : LBF SEC^2/IN^4
%  Scalar Value :  7.317372E-04 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : RANKINE
%  Scalar Value :  5.310000E+02 
%------------------------------------------------------------------------------
%Material Property      : 1 - MODULUS OF ELASTICITY
%  Dim : PRESSURE    Unit : KILOGRAM FORCE/MM^2
%  Scalar Value :  2.108773E+04 
%Material Property      : 3 - MASS DENSITY
%  Dim : MASS DENSITY    Unit : KGF SEC^2/MM^4
%  Scalar Value :  7.974181E-10 
%Material Property      : 7 - THERMAL EXPANSION REFERENCE TEMPERATURE
%  Dim : TEMPERATURE    Unit : KELVIN
%  Scalar Value :  2.950000E+02 
%------------------------------------------------------------------------------

