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

%P_SHELL Shell element property function
%
%       Syntax : mat= p_shell('default') 
%                mat= p_shell('database ProId name') 
%                il = p_shell('dbval ProId name');
%                il = p_shell('dbval -unit MM ProId name');
%
%       Supported subtypes are
%       1 : standard isotropic
%	    [Id# Type f   d 0 h   k MID2 12I/T^3 MID3 NSM Z1 Z2 MID4]
%          with
%	     Type : fe_mat('p_shell','SI',1)
%            f    : formulation (see the element : quad4, ...)
%            d    : -1 no drilling stiffness, If d==0 d is set to 1.
%		    d>0 drilling DOF stiffness coefficient
%            h    : thickness
%            k    : shear factor (default value is 5/6)
%            12I/T^3 : Ratio of bending moment of inertia to nominal T^3/12
%                      (default 1)
%            NSM  : non structural mass per unit area
%            MID2, MID3, Z1, Z2, MID4 currently unused
%
%       2 : composite (follows NASTRAN's PCOMP format)
%           [ProID type   Z0 NSM SB d TREF GE f MatId1 T1 Theta1 SOUT1 ...]
%           Type fe_mat('p_shell','SI',2)
%           z0   distance from the reference plane to the bottom surface 
%                default -(total thickness)/2
%           SB   allowable shear stress (unused)
%           d    drilling coefficient (see above)
%           TREF reference temperature (unused)
%           GE   loss factor
%           f    formulation
%           MatId1,T1,Theta1,SOUT1 : material idenfier for the layer,
%                layer thickness, orientation, layer options (unused)
%
%       For automated generation of properties see, the P_SHELL
%       dbval and database commands.
%
%       See also help fe_mat
%                doc  p_shell, fem (handling materials section), pl, fe_mat

%	Etienne Balmes, 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.49 $  $Date: 2006/05/12 17:25:09 $'; return;
end
if nargin<1 help p_solid;return; end
if ischar(varargin{1}); [CAM,Cam]=comstr(varargin{1},1); il=[]; carg=2;
else il=varargin{1};[CAM,Cam]=comstr(varargin{2},1); carg=3;
end

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

 i1=varargin{carg}; 
 out1={};
 switch i1 % PropertySubType
 case 1 % standard fe_mat('unitlabel') 
     st={...
         'ProId'             0;
         'Type'              0;
         'Formulation'       0;
         'Drilling DOF'      0; 
         'unused'            0;
         't (thickness)'     4;
         'k (shear factor)'  0;
         'MID2'              0;
         '12I/T3'            0;
         'MID3'              0;
         'NSM'               9;
         'Z1'                4;
         'Z2'                4;
         'MID4'              0
          };
 case 2 % composite
     st={...
         'ProId'             0;
         'Type'              0;
         'Z0'                4;
         'NSM'               9;
         'SB'                0;
         'FT'                0;
         'TREF'              7;
         'GE'                0;
         'LAM'               0;          
         'MatId1'            0;
         'T1'                4;
         'Theta1'            0;
         'SOUT1'             0};
     out1={         'MatId%i'            0;
         'T%i'                4;
         'Theta%i'            0;
         'SOUT%i'             0};
 otherwise; st={'ProId' 0; 'Type', 0};
 end
 if ~isempty(strfind(Cam,'cell')); out=st; else; out=[st{:,2}]; end


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

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

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

 i1=strfind(comstr(Cam,-27),'-unit'); out1={};
 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

 while 1==1
  [CAM,Cam]=comstr(CAM,6);
  [i1,CAM,Cam]=comstr(CAM,'','%i');
  if isempty(CAM)&carg<=nargin st=varargin{carg};carg=carg+1; 
  else st=CAM;end
  [mat,st1,i2]=p_shell('database',st,varargin{carg:end});carg=carg+i2-3;
  if ~isempty(Unit)
   mat.il=fe_mat(sprintf('convert %s %s',mat.unit,Unit),mat.il);
  end
  r1=mat.il; if length(i1)==1 r1(1)=i1;end
  if ~isempty(il) i2=find(il(:,1)==r1(1)); else i2=[];end
  if isempty(i2) i2=size(il,1)+1;end
  il(i2,1:length(r1))=r1;
  if carg>nargin break;end
  [CAM,Cam]=comstr(varargin{carg},1);carg=carg+1;
 end
 out=il;

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


  st=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.il=[MatId fe_mat('p_shell',1,1) 1 1 0 .01]; 
  out.name='Kirch';
  out.type='p_shell';
  out.unit='SI';

  out(2).il=[MatId fe_mat('p_shell',1,1) 2 1 0 .01];
  out(2).name='Mind';
  out(2).type='p_shell';
  out(2).unit='SI';

  i1=strmatch(comstr(st,-27),comstr({out.name},-27));
  if ~isempty(i1) out=out(i1);
  elseif ~isempty(st)
   [out,carg]=create_section_shell(st,out,varargin,carg);  out.il(1)=MatId;
  end
  out1='Shell properties';out2=carg;

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

 if carg>nargin; opt=integrules('quadb');
 else; opt=varargin{carg};carg=carg+1;
 end
 if ~ischar(opt);
 elseif comstr(opt,'tria3'); opt=integrules(opt,-1); % tria3 at center
 else; opt=integrules(opt);
 end
 if carg>nargin; integ=[];else; integ=varargin{carg};carg=carg+1;end
 opt.Nw=size(opt.N,1);
 if ~isfield(opt,'Nnode'); opt.Nnode=size(opt.N,2);end
 if length(integ)<6;rule=[1 opt.Nw];
 else; rule=integ(5:6);rule=rule(:)';
 end

 opt.DofLabels={'u','v','w','\theta_u','\theta_v','\theta_w'};
 % Strain energy
 % define the deformation vector: row, NDN, DDL, NwStart, NwRule
   % define the deformation vector: row, NDN, DDL
   if size(opt.w,1)==13; rule1=[5 9];rule=[1 4];else; rule1=rule;end
   r2=[1 2 1 rule1; 2 3 2 rule1;3 3 1 rule1;3 2 2 rule1;% e_xx =N,x u,e_yy,e_xy
     4 -2 5 rule1;5 3 4 rule1;6 2 4 rule1;6 -3 5 rule1; % curvatures
     7 2 3 rule;7 1 5 rule;8 3 3 rule;8 -1 4 rule];  % shear constraint
   if size(opt.w,1)==13;
    r2=[r2; ...
    9 1 6 rule;9 1 6 rule;9 -2 2 rule;9 3 1 rule]; % dril 2\theta_z-(v,x-u,y)
   end
   % Topology
   % [rho eta Gii] in constit -> shift values of col 5 by two
    dd=[ones(6) zeros(6,2);zeros(2,6) ones(2)];dd(9,9)=1;
    dd(find(dd))=[1:length(find(dd))]+8; % see integinfo call
 opt.StrainDefinition{1}=r2;
 opt.StrainLabels{1}={'e_x','e_y','e_xy','k_x','k_y','k_xy','s_x','s_y','d_z'};
 opt.ConstitTopology{1}=dd;

 % Kinetic energy
 opt.StrainDefinition{2}= [1 1 1 rule;2 1 2 rule; 3 1 3 rule];
 opt.StrainLabels{2}={'u','v','w'};
 opt.ConstitTopology{2}=eye(3);

 opt=integrules('matrixrule',opt);
 opt.VectMap=reshape(1:6*opt.Nnode,6,opt.Nnode)'; 

 %display with :integrules('texstrain',opt);

 if carg<=nargin
  constit=varargin{carg};carg=carg+1;
  % Only keep terms needed for the current law.
  i1=find(any(constit,2))-1; 
  for j0=1:length(opt.MatrixIntegrationRule);
    r1=opt.MatrixIntegrationRule{j0};if ~isempty(r1)
     opt.MatrixIntegrationRule{j0}=r1(find(ismember(double(r1(:,5)),i1)),:);
    end
  end
 end
 if nargout==0
  integrules('texstrain',opt);
  try; opt=integrules('stressrule',opt);integrules('texstress',opt);end
 else; out=opt;
 end
 out1=23; % 3d surface


% -------------------------------------------------------------------------
%[constit,integ]=p_shell('buildconstit',ID,pl,il);
elseif comstr(Cam,'buildconstit')

ID=varargin{carg};carg=carg+1;
pl=varargin{carg};carg=carg+1;
il=varargin{carg};carg=carg+1;il0=il;
il=il(find(il(:,1)==ID(2)),:);

if size(il,1)>1; error('BuildConstit assumes one p_shell');end

dd=zeros(6); ds=zeros(2); % dd: [A B;B D]  DS shear 
rhoh=0; etah=0; out3=[];

if isempty(il); sdtw('_err','No element property row matching %i',ID(2));
elseif il(2)==2; st='p_shell';typ=1; unit='US'; % OBSOLETE
else; [st,unit,typ]=fe_mat('typep',il(2));
end

if strcmp(st,'p_shell')

switch typ
case 1 % Standard isotropic plate - - - - - - - - - - - - - - - - - -

 mat=pl(find(pl(:,1)==ID(1)),:);mat=mat(1,:);
 [st,unit,typ]=fe_mat('typem',mat(2));
 if ~strcmp(st,'m_elastic')|typ~=1; 
  sdtw('MatId %i assuming m_elastic subtype 1',ID(1));
 end
  % 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

 if isempty(mat) error(sprintf('MatId %i not found',ID(1)));end
 il(15)=0;
 if il(4)==0 il(4)==1; end % drilling DOF stiffness
 o1 = mat; o2 = il;
 o3 = mat(3)/(1-mat(4)^2); % membrane
 o3 = [o3 mat(4)*o3 0;      mat(4)*o3 o3 0;      0 0 (1-mat(4))/2*o3]*il(6);
 o4 = o3*il(6)^2/12;     % bending
 dd=[o3 o3*0;o3*0 o4];

 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 %i shell declaration',ID(2));
 end
 if length(il)<7 il(7)=5/6; elseif il(7)==0 il(7)=5/6;end
 ds = eye(2)*mat(6)*il(7)*il(6); % shear

 ddril=mean(diag(dd(4:end,4:end)));
 if il(4)==-1; ddril=0; elseif il(4)==0;ddril=ddril*1e-4;
 else; ddril=ddril*il(4);
 end

 %          rho*h       eta       f d h k 12I/t3 nsm   dd ds
 constit = [mat(5)*il(6) mat(5) il([3 4 6 7 9  11 ]) dd(:)' ds(:)' ddril];
 out3=[mat(1) -il(6)/2 il(6)/2 0]; % one layer
 % xxx deal with 5 DOF elmap for mitc4 element xxx

case 2 % Composite specified as PCOMP entry - - - - - - - - - - - - - - - - -

r1=il(10:end);i1=ceil(length(r1)/4)*4;if length(r1)<i1;r1(i1)=0;end
r1=reshape(r1,4,length(r1)/4)';  % each row defines a ply
r1=r1(find(any(r1,2)),:);
z=cumsum([0;r1(:,2)])+il(3);
other=[];

for j1=1:size(r1,1) % number of plies

   mat=pl(find(pl(:,1)==r1(j1,1)),:);
   [st,unit,typ]=fe_mat('type',mat(2));
   cz=z(j1:j1+1); 
   if any(exist(st)==[2 3 6])
     [S,rhoh,r2]=feval(st,'buildply',mat,pl,il0,cz,rhoh);
     if ~isempty(r2); other=[other;r2(:)];end
   else; error('Ply not implemented for this type of material');
   end 
  out3(j1,1:4)=[mat(1) cz(:)' r1(j1,3)]; % one layer
    %Stiffness in the principal coordinates
    C=pinv(S);
    % Stiffness matrix obtained in a arbitrary coordinate by rotation of
    % angle theta around the z axis of the main coordinates 
    % Q=(Tsigma)^-1*C*Tepsilon

    %coordinate transformation matrix
    teta=-r1(j1,3)*pi/180; ct=cos(teta); st=sin(teta);
    % ITsigma=inv(Tsigma) 
    iTsigma=[ct^2 st^2 2*ct*st 0 0;st^2 ct^2 -2*ct*st 0 0; ...
     -ct*st ct*st ct^2-st^2 0 0; 0 0 0 ct -st;0 0 0 st ct];
    Tepsilon=[ct^2 st^2 -ct*st 0 0;st^2 ct^2 ct*st 0 0; ...
     2*ct*st -2*ct*st ct^2-st^2 0 0;0 0 0 ct st;0 0 0 -st ct];
    Q=iTsigma*C*Tepsilon; 
    Q1=Q([1 6 11;2 7 12;3 8 13]);Q2=Q([19 24;20 25]);


    dd=dd+[diff(cz)*Q1 diff(cz.^2)/2*Q1;diff(cz.^2)/2*Q1 diff(cz.^3)/3*Q1];
    ds=ds+5/6*diff(cz)*Q2;

end % loop on plies
% detect zero terms
%[i1,i2]=find(dd); r2=diag(dd); 
%dd(find(abs(dd(:)./sqrt(r2(i1).*r2(i2)))<eps));
 ddril=mean(diag(dd(4:end,4:end)))*1e-4;
%          rho*h       eta       f d h k 12I/t3 nsm   dd ds dril(49)
constit = [rhoh etah il(9) 0 z(end)-z(1) 5/6 0 0  dd(:)' ds(:)' ddril];
if ~isempty(other);constit=[constit(:);other(:)];end
if sdtdef('diag')>1
 fprintf('\nMatId %i ProId %i, formulation %i\n',ID(1:2),constit(3))
 fprintf('Shell constitutive matrix membrane and bending\n')
 disp(dd)
 fprintf('Shear constitutive matrix\n')
 disp(ds)
 fprintf('Rho*h = %.3g, h= %.3g, dril = %.3g\n',constit([1 5 49]))
end

otherwise; error('Not a supported p_shell subtype');
end

else % not a p_shell property
 [st,unit,typ]=fe_mat('typep',il(2));
 [constit,ID]=feval(st,'buildconstit',ID,pl,il0);
end

out=constit(:);     % Return CONSTIT
out1=int32(ID(:));  % Return integ
out2=[]; % elmap

% -------------------------------------------------------------------------
%RunOpt=p_shell('BuildDofOpt',RunOpt,pl,il);
elseif comstr(Cam,'builddofopt')

RunOpt=varargin{carg};carg=carg+1;
pl=varargin{carg};carg=carg+1;
il=varargin{carg};carg=carg+1;
RunOpt.FieldDofs=[1:6]; 
out=RunOpt;

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

 i1=varargin{carg}; carg=carg+1;
 switch i1
 case 1;  out='Std_shell';
 case 2;  out='Comp_shell';
 otherwise out='p_shell';
 end

end

% ---------------------------------- internal functions
% --------------------------------------------------------------------
% Creation of standard section properties
function [out,carg] = create_section_shell(CAM,out,varg,carg)

[CAM,Cam]=comstr(CAM,1);
RunOpt=struct('Formulation',0,'Drilling',0);
i1=strfind(Cam,'-f'); if ~isempty(i1) % Formulation
 [RunOpt.Formulation,i2,i3,i4]=sscanf(Cam(i1+2:end),'%i',1);
 CAM=CAM([1:i1-1 i1+2+i4:end]);[CAM,Cam]=comstr(CAM,1);
end

RunOpt.unit='SI';typ=fe_mat('p_shell',RunOpt.unit,1);

if comstr(Cam,'kirchoff'); 
  if RunOpt.Formulation==0; RunOpt.Formulation=1;end
  r1=[RunOpt.Formulation RunOpt.Drilling 0 comstr(Cam,'kirchoff','%g')];
elseif comstr(Cam,'kirchhoff');  
  if RunOpt.Formulation==0; RunOpt.Formulation=1;end
  r1=[RunOpt.Formulation RunOpt.Drilling 0 comstr(Cam,'kirchhoff','%g')];
elseif  comstr(Cam,'mindlin')
  if RunOpt.Formulation==0; RunOpt.Formulation=2;end
  r1=[RunOpt.Formulation RunOpt.Drilling 0 comstr(Cam,'mindlin','%g')]; 
elseif  comstr(Cam,'laminate')
 r1=comstr(Cam,'laminate','%g');
 r1=reshape(r1,3,length(r1)/3);r1(4,:)=0;
 z=cumsum([0 r1(2,:)]);
 r1=[-z(end)/2 0 0 0 0 0 0 r1(:)'];
 typ=fe_mat('p_shell','SI',2);

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

[CAM,Cam]=comstr(CAM,10); opt=comstr(CAM,-1);
mat=varg{carg}; carg=carg+1;
if length(opt)<9 opt(9)=0; end % use G1z_min;

%opt=[hpeau1,hpeau2,hcellule,a,b,teta,t,t1,alpha_minmax]
%	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)

[st,RunOpt.unit,typ]=fe_mat('typem',mat.pl(1,2));
if ~strcmp(st,'m_elastic')|(typ~=1)
 error('Honeycomb supposes an elastic material');
end

o6=opt(6)/180*pi;

r1=[ ...
 mat.pl(3)*(opt(7)/opt(5))^3*((opt(4)/opt(5)+sin(o6))/(cos(o6))^3); %E1
 mat.pl(3)*(opt(7)/opt(5))^3*cos(o6)/((opt(4)/opt(5)+sin(o6))*sin(o6)^2); %E2
 ((opt(4)/opt(5)+sin(o6)))*sin(o6)/cos(o6)^2; % nu12
 mat.pl(3)*(opt(4)/opt(5)+sin(o6))/(2*(opt(4))^2*cos(o6))*1/(opt(5)/(2*(opt(8))^3)+opt(4)/opt(7)^3); % G12
mat.pl(6)*(opt(4)*opt(8)+opt(5)*opt(7)*sin(o6))^2/(opt(5)*cos(o6)*(opt(4)+opt(5)*sin(o6))*(2*opt(4)*opt(8)+opt(5)*opt(7))) ; %G1z
 mat.pl(6)*cos(opt(6))/(opt(4)/opt(5)+sin(opt(6)))*opt(7)/opt(5)]'; % G2z

r1(end-1)=r1(end-1)*(1-opt(9))+opt(9)* ...
 mat.pl(6)*(opt(4)*opt(8)+2*opt(5)*opt(7)*(sin(o6))^2)/ ...
  (2*opt(5)*cos(o6)*(opt(4)+opt(5)*sin(o6)));
% voir comment choisir G1z en fonction de G1z_inf et G1z_sup.

 typ=fe_mat('p_shell',mat.unit,5);

else 
  warning('Using default shell property'); out=out(1);
  return
end

 out=struct('il',[1 typ r1],'name',Cam,'type','p_shell', ...
   'unit',RunOpt.unit);

% EOF