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

%P_SOLID volume element property
%
%       Syntax : mat= p_solid('default') 
%                mat= p_solid('database matid (command)') 
%                il = p_solid('dbval ProId (command)');
%                il = p_solid('dbval -unit MM ProId (command)');
%
%       Accepted commands for database and dbval calls are
%         d3 2 : 2x2x2 integration rule for linear volumes (hexa8 ... )
%         d3 3 : 3x3x3 integration rule for quadratic volumes (hexa20 ... )
%         d2 2 : 2x2x2 integration rule for linear volumes (q4p ... )
%         d2 3 : 3x3x3 integration rule for quadratic volumes (q8p ... )
%
%       Properties supported by p_solid are
%       Subtype 1 : 3D solid volume element: hexa8, hexa20, penta6, tetra4, ...
%            [ProId Type CordM Integ Stress Isop] 
%          CordM : material coordinate system (unused)
%          Integ : if non zero, selects the integration rule in 
%            integrules(ElemP,Integ) for example integrules(hexa8('parent'),3)
%            if Integ=0 or property undefined old of_mk_subs elements are used 
%          Stress : stress integration rule (unused)
%          ISOP : geometry interpolation flag (unused)
%
%       Subtype 2 : 2D volume element : q4p, q8p, t3p, t6p, ... 
%         [ProId Type Form N Integ]
%          with 
%           Type = fe_mat('p_solid','SI',2)
%           Form : formulation (0 plane strain, 1 plane stress, 2 axisymetric)
%           N    : Fourier harmonic for axisymetric elements that support it
%           Integ: see above
%
%       Subtype 3 : ND-1 volume, surface load or pressure
%         [ProId fe_mat('p_solid','SI',3) Integ Type Ndof1 ...]
%           Integ see above
%           Type: 1 volume force, 2 density of volume force
%                 3 pressure,     4 fluid/structure coupling
%
%       See also : help fe_mat
%                  doc  fem (handling materials section), pl, fe_mat, p_shell


%	Jean-Michel Leclere, Etienne Balmes
%       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.82 $  $Date: 2006/05/29 11:14:00 $'; 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,'default')

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

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

 r1=p_solid('database');fprintf('%s\n',r1.name);

% -------------------------------------------------------------------------
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
  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=p_solid('database',st); 
  elseif isnumeric(st)
    [typ,st1,i4]=fe_mat('typep',st(2));
    mat=struct('il',st,'name',sprintf('%i',st(1)),'type',typ,'unit',st1);
  end
  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=varargin{carg};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);

  % MatId Typ                     CordM Integ Stres Isop Field
  out.il=[MatId fe_mat('p_solid','SI',1) 0 2     0     1    0]; 
  out.name='Full 2x2x2';
  out.type='p_solid';
  out.unit='SI';

  mat.il=[MatId fe_mat('p_solid','SI',1) 0 2 0 1 0]; 
  mat.name='Reduced shear';
  mat.type='p_solid';
  mat.unit='SI';
  out(end+1)=mat;

  mat.il=[MatId fe_mat('p_solid','SI',2) 1 0 0]; 
  mat.name='Plane stress';
  mat.type='p_solid';
  mat.unit='SI';
  out(end+1)=mat;

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

  if ~isempty(i1) out=out(i1);
  elseif comstr(comstr(st,-27),'d3')
   r1=comstr(st(3:end),[-1 3]);
   if length(r1)==1; r1=[0 r1 0 0];end
   out=struct('il',[MatId fe_mat('p_solid','SI',1) r1(:)'], ...
         'name',horzcat('Solid',st),'type','p_solid','unit','SI');
  elseif comstr(comstr(st,-27),'d2')
   r1=comstr(st(3:end),[-1 3]);
   if length(r1)==1; r1=[0 0 r1 0];end
   out=struct('il',[MatId fe_mat('p_solid','SI',2) r1(:)'], ...
         'name',horzcat('Solid',st),'type','p_solid','unit','SI');
  elseif ~isempty(st)
   r1=str2num(st); 
   if length(r1)<1 error('Not a consistent database call');end
   out=[];
   out.il=[MatId fe_mat('p_solid','SI',1) r1(:)']; 
   out.name=['Solid' st];
   out.type='p_solid';
   out.unit='SI';
  end

  out1='Solid';

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

 i1=varargin{carg};
 out1={};
 switch i1 % PropertySubType
 case 1 % [ProID type Coordm In Stress Isop Fctn  ]
   st={ ...
   'ProId'   0;
   'Type'    0;
   'COORDM'  0;
   'IN'      0;
   'STRESS'  0;
   'ISOP'    0;
   'FCTN'    0 };
 case 2 % 2D   [ProId Type Form N]
   st={ ...
   'ProId'   0;
   'Type'    0;
   'Form'    0;
   'N'       0 };
 otherwise; st={'ProId' 0; 'Type', 0};
 end
 if ~isempty(strfind(Cam,'cell')); out=st; else; out=[st{:,2}]; end

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

 i1=varargin{carg}; carg=carg+1;
 switch i1
 case 1;  out='p_solid';
 otherwise out='p_solid';
 end

% -------------------------------------------------------------------------
% [EltConst,NDNDim]=p_solid('Const',ElemF,integ,constit);
elseif comstr(Cam,'const')

 opt=varargin{carg};carg=carg+1;
 integ=varargin{carg};carg=carg+1;
 if carg<=nargin; constit=varargin{carg};carg=carg+1;
 else constit=[];
 end
 if ~isempty(constit)&constit(1)==-1; 
  st=fe_mat('typep',constit(2,1));
  [out,out1]=feval(st,'const',varargin{2:end});
  return;
 end
 if ischar(opt) % Allow for integrule switching here
  if size(integ,1)>4; 
    opt=integrules(opt,double(integ(5,1)));
    if any(integ(5,:)~=integ(5,1));
     warning('Cannot deal with multiple integration strategies in the same group');
    end
  else; opt=integrules(opt);
  end
 end

 opt.Nw=size(opt.N,1);out.Nnode=size(opt.N,2);
 if size(integ,1)<6|any(integ(5:6,1)==0)|integ(6)>size(opt.w,1);
       rule=[1 opt.Nw]; % By default start at 1 use all points
 else; rule=integ(5:6,1);rule=rule(:)';
 end
 if ~any(rule); 
   rule=[1 opt.Nw];
   sdtw('integ is assumed to give standard integration rule in integ(5:6)');
 elseif rule(1)==0|rule(end)>size(opt.w,1);error('Inconsistent rule'); 
 end

 % 2D SOLID - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 if integ(3)==2*double(integ(4)) 

 % Strain energy
 % define the deformation vector: row, NDN, DDL, NwStart, NwRule
 opt.StrainDefinition{1}= ...
   [1 2 1 rule; 2 3 2 rule; % e_xx = N,x u, e_yy
    3 3 1 rule;3 2 2 rule]; % shear
 opt.StrainLabels{1}={'\epsilon_x','\epsilon_y','\gamma_{xy}'};
 opt.ConstitTopology{1}=int32(reshape(3:11,3,3));
 % fprintf('D=\n');disp(constit(opt.ConstitTopology{1}));
 opt.DofLabels={'u','v'};
 out1=2; % Tell that BuildNDN rule is 2D
 % Kinetic energy
 opt.StrainDefinition{2}= [1 1 1 rule;2 1 2 rule];
 opt.StrainLabels{2}={'u','v'};
 opt.ConstitTopology{2}=eye(2);
 opt=integrules('matrixrule',opt);
 % \int f v (see elem0 for format)
 opt.RhsDefinition=int32( ...
   [101 0 2 0     0 0 -1    rule+[-1 0];
    101 1 2 0     1 0 -1    rule+[-1 0]]);
 %display with :integrules('texstrain',opt);

 % 3D SOLID - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 elseif integ(3)==3*double(integ(4)) 
 % Strain energy
 % define the deformation vector: row, NDN, DDL, NwStart, NwRule
 opt.StrainDefinition{1}= ...
   [1 2 1 rule; 2 3 2 rule;3 4 3 rule; % e_xx = N,x u, e_yy, e_zz
   4 4 2 rule;4 3 3 rule;5 4 1 rule;5 2 3 rule;6 3 1 rule;6 2 2 rule]; % shear
 opt.StrainLabels{1}={'\epsilon_x','\epsilon_y','\epsilon_z', ...
     '\gamma_{yz}','\gamma_{zx}','\gamma_{xy}'};
 opt.ConstitTopology{1}=int32(reshape(3:38,6,6));
 opt.DofLabels={'u','v','w'};
 out1=3; % Tell that BuildNDN rule is 3D

 % 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);
 % \int f v (see elem0 for format)
 opt.RhsDefinition=int32( ...
   [101 0 3 0     0 0 -1    rule+[-1 0];
    101 1 3 0     1 0 -1    rule+[-1 0];
    101 2 3 0     2 0 -1    rule+[-1 0]]);
 %display with :integrules('texstrain',opt);

 opt.VectMap=reshape(1:3*opt.Nnode,3,opt.Nnode)'; 

 % Switch on NL constitutive materials - - - - - - - - - - - - - - -
 if size(integ,1)>6; switch integ(7,1) 
 % this should match StrategyType in of_mk.c matrix assembly
 case 3; opt.material='hyperelastic';
 case 901; 
   opt.material='heart';
   cEGI=evalin('caller','varargin{end}');
   gstate=zeros(7+opt.Nnode*length(opt.DofLabels),opt.Nw*length(cEGI));
   opt.gstate = gstate;
 otherwise 
   opt.material='Elastic3DNL';
 end;end
 % 2D fluid : as many DOFs as nodes - - - - - - - - - - - - - - - - - - - -
 elseif integ(3)==double(integ(4))&~any(opt.w(:,3)) ...
   &(size(integ,1)<=4 |integ(7)==2)
 opt.DofLabels={'p'}; out1=2;% Tell that BuildNDN rule is 2D

 % Strain energy
 % define the deformation vector: row, NDN, DDL, NwStart NwRule
 opt.StrainDefinition{1}=[1 2 1 rule;2 3 1 rule];
 opt.StrainLabels{1}={'p,x','p,y'};
 opt.ConstitTopology{1}=diag([3 3 3]);

 % Kinetic energy
 opt.StrainDefinition{2}= [1 1 1 rule];
 opt.StrainLabels{2}={'p'}; opt.ConstitTopology{2}=1;

 opt=integrules('matrixrule',opt);
 % 3D fluid : as many DOFs as nodes - - - - - - - - - - - - - - - - - - - -
 elseif integ(3)==double(integ(4))&(size(integ,1)<=4|integ(7)==3)

 opt.DofLabels={'p'};
 out1=3; % Tell that BuildNDN rule is 3D
 % Strain energy
 % define the deformation vector: row, NDN, DDL, NwStart NwRule
 opt.StrainDefinition{1}=[1 2 1 rule;2 3 1 rule;3 4 1 rule];
 opt.StrainLabels{1}={'p,x','p,y','p,z'};
 opt.ConstitTopology{1}=diag([3 3 3]);

 % Kinetic energy
 opt.StrainDefinition{2}= [1 1 1 rule];
 opt.StrainLabels{2}={'p'}; opt.ConstitTopology{2}=1;

 opt=integrules('matrixrule',opt);

 elseif constit(1)==-3  % interface element
   if constit(4,1)==4 % Fluid/structure coupling
     opt.DofLabels={'u','v','w','p'};
     opt.material='fs_matrix'; opt.bas=zeros(9,size(opt.N,1));
     Ndof=4*opt.Nnode;
     opt.VectMap=int32([1:4:Ndof 2:4:Ndof 3:4:Ndof 4:4:Ndof])';
     out1=23;
   end
 else; 
  fprintf('Not a known p_solid(''const'') case'); error('Not a valid case');
 end % 2D or 3D - - - - - - - - - - - - - - - - - - - - - - - - - - -
 if size(opt.N,2)~=opt.Nnode;
  eval('opt1=p_piezo(''ReducedShear'',opt,integ,constit);opt=opt1;','')
 end
 if ~isempty(constit)&isfield(opt,'MatrixIntegrationRule')
  % 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
 else;opt.MatrixIntegrationRule={};
 end

 if nargout==0
  integrules('texstrain',opt);
  try; opt=integrules('stressrule',opt);integrules('texstress',opt);end
 else; out=opt;
 end


% -------------------------------------------------------------------------
% Implementation of elastic constitutive law building 3D
%[constit,integ]=p_solid('buildconstit',ID,pl,il);
% This is called by the element 'integinfo' command
elseif comstr(Cam,'buildconstit')

  ID=varargin{carg};carg=carg+1;
  pl=varargin{carg};carg=carg+1; mat=pl(find(pl(:,1)==ID(1)),:);
  il=varargin{carg};carg=carg+1; pro=[];try;pro=il(find(il(:,1)==ID(2)),:);end
  if size(mat,1)>1; 
   sdtw('_nb','Using first material with ID %i',ID(1));
   mat=mat(1,:);
  end
  if all(size(pro)>[0 1])
    [RunOpt.ProF,RunOpt.ProU,RunOpt.pTyp]=fe_mat('typep',pro(1,2));
  else; RunOpt.ProF='m_null';
  end

  if isempty(mat); 
   sdtw('_nb','MatId %i is not matched',ID(1));
  else;  [st,unit,typ]=fe_mat('typem',mat(2));
   % 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
  end

  % obsolete calls to of_mk_subs elements (not )
  if isequal(RunOpt.ProF,'p_solid')&RunOpt.pTyp==3|typ==2
  % new element properties should not end up using the default mechanical
  elseif ~isequal(RunOpt.ProF,'p_solid')&~isequal(RunOpt.ProF,'m_null')
  elseif length(ID)==4&((ID(3)/ID(4)==3&(length(pro)<4|pro(4)==0))| ...
    (ID(3)/ID(4)==2&(length(pro)<5|pro(5)==0)))
    if ID(3)/ID(4)==3;     [constit,out1]=fe_mat('plil of_mk 3d',ID,mat,il);
    fprintf('Fixed rule linear 3D solid p_solid Mat %i /Pro %i\n',ID(1:2));
    elseif ID(3)/ID(4)==2; [constit,out1]=fe_mat('plil of_mk 2d',ID,mat,il);
    fprintf('Fixed rule linear 2D solid p_solid Mat %i /Pro %i\n',ID(1:2));
    else; error('Not an expected case');
    end
    out1=out1(:);out1(3:4)=ID(3:4);
    i1=find(triu(ones(ID(3),ID(3)))); r1=zeros(ID(3),ID(3));r1(i1)=1:length(i1);
    out=constit(:);out2=r1+tril(r1',-1);
    % return Constit, Integ, Elmap
    return; 
  end % obsolete calls to of_mk_subs elements
  if isempty(pro); error('ProId is not matched');end

   RunOpt.Failed=0;

  % 3D solid -----------------------------------------------------------
  if isequal(RunOpt.ProF,'p_solid')&RunOpt.pTyp==1

   if strcmp(st,'m_elastic')&typ==1 % Standard isotropic
    
    r1=mat(6); % G
    r2=mat(3)*mat(4)/(1+mat(4))/(1-2*mat(4)); % En/(1+n)(1-n)
    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])=r1; % G
    %E=mat(3);n=mat(4);G=E/2/(1+n);
    %inv([1/E -n/E -n/E 0 0 0;-n/E 1/E -n/E 0 0 0;-n/E -n/E 1/E 0 0 0;
    %     0 0 0 1/G 0 0; 0 0 0 0 1/G 0;0 0 0 0 0 1/G])-dd
    if any(diag(dd)<0); sdtw('Negative diagonal term in constitutive law MatId %i',mat(1));end
    constit=[mat(5);mat(7);dd(:)]; % 1:38 fully anisotropic constitutive law
    if size(mat,2)>=9&any(mat(8:9))  % Thermal expansion and ref temp
     at=eye(3)*mat(8); constit=[constit;at(:);mat(9)];  % 39:48
    end

   elseif strcmp(st,'m_elastic')&typ==2 % 3D fluid Acoustic 

    if length(mat)<5; mat(5)=0.; end
    %  constit=[1/rho/C2 eta 1/rho]
    constit = [1/mat(1,3)/mat(1,4)^2 mat(1,5) 1/mat(1,3)]';
    ID(3)=ID(4);ID(7)=3;

   elseif strcmp(st,'m_elastic')&typ==3 % 3D anisotropic materials

    if length(mat)<25; mat(25)=0;end
    dd=triu(ones(6));dd(find(dd))=[3:23];dd=dd+triu(dd,1)';
    dd=mat(dd) 
    constit=[mat(24);mat(25);dd(:)];
    % here thermal expansion is missing 

   else; RunOpt.Failed=1;
   end
   % Integration rule saved in integ
   if size(pro,1)&size(pro,2)>3;  ID(5)=pro(1,4); end  
  % 2D solid -----------------------------------------------------------
  elseif isequal(RunOpt.ProF,'p_solid')&RunOpt.pTyp==2
   if strcmp(st,'m_elastic')&typ==1 % Standard isotropic

     %if o2(8)==0 o2(8)=o2(7);end 
     if size(mat,2)<7; mat(7)=0;end;constit = mat([5 7]); % [rho eta]
     E =mat(3); nu=mat(4); 

     switch pro(3) % Formulations
     case 1 % - - - - - - - - - -  - - - - - - - plane stress
       C=E/(1.-nu^2);r1=[constit C C*nu C 0. 0. C*(1-nu)/2];
       constit=r1([1 2 3 4 6 4 5 7 6 7 8]);
     case 0 % - - - - - - - - - -  - - - - - - - plane strain
       C=E*(1-nu)/(1+nu)/(1-2*nu);
       r1=[constit C C*nu/(1-nu) C 0. 0. C*(1-2*nu)/2/(1-nu)];
       constit=r1([1 2 3 4 6 4 5 7 6 7 8]);
     case 2 % - - - - - - - - - -  - - - - - - - axisymetric
       error('Axisymmetry is not supported by the *b element family');
     end % Formulations 

   % 2D Anisotropic : car = [rho eta E11 ... E33 a1 a2 a3 T0]
   elseif strcmp(st,'m_elastic')&typ==4 
    if length(pl)<9 pl(14)=0.; end 
    r1=mat([9 14 3:8 10:13]);constit=r1([1 2 3 4 6 4 5 7 6 7 8]);
   % 2D fluid Acoustics
   elseif strcmp(st,'m_elastic')&typ==2
    if length(mat)<5; mat(5)=0.; end
    %  constit=[1/rho/C2 eta 1/rho]
    constit = [1/mat(1,3)/mat(1,4)^2 mat(1,5) 1/mat(1,3)]';
    ID(3)=ID(4);ID(7)=2; 
   else; RunOpt.Failed=1;
   end
   % Integration rule saved in integ
   if size(pro,1)&size(pro,2)>=5;  ID(5)=pro(1,5); end  
  % Surface element -----------------------------------------------------------
  elseif isequal(RunOpt.ProF,'p_solid')&RunOpt.pTyp==3
   % Integration rule saved in integ
   ID(5)=pro(1,3);
   if pro(1,4)==4 % fluid/structure 
     ID(3)=4*ID(4); % 4 dof per node
   elseif pro(1,4)==3; ID(3)=3*ID(4); % pressure (3dof per node)
   else; error('Not a valid formulation');
   end
   constit=[-3 pro(2:end)];constit=constit(:);
  else; RunOpt.Failed=1; % 2D, 3D surface
  end

  if RunOpt.Failed % Attempts to call other p_* functions
    try; % Should be in m_* but if user wants to limit foot print p_* OK too
     eval(sprintf('[out,out1]=%s(varargin{:});',st));
    catch
     try; st(1)='p';eval(sprintf('[out,out1]=%s(varargin{:});',st));
     catch;
       if sdtdef('diag')>10; 
        eval(sprintf('[out,out1]=%s(varargin{:});',st));
       end
       error(sprintf( ...
        horzcat('Call to %s(''BuildConstit'') failed', ...
         ' for subtype (%i) with \n%s\nUse sdtdef(''diag'',12) to debug '), ...
         st,typ,lasterr));
     end
    end
    out2=elem0('elmapmat_og', ...
      [double(out1(4)) double(out1(3))/double(out1(4))]);
    return;
  end% Attempts to call other p_* functions
  out=constit(:); out1=int32(ID(:));
  if out1(4)==0;out2=[];
  else
   out2=elem0('elmapmat_og',[double(out1(4)) double(out1(3))/double(out1(4))]);
  end


% -------------------------------------------------------------------------
% Advanced multi-physic dof building
elseif comstr(Cam,'builddof')

  model=varargin{carg};carg=carg+1;
  cEGI=varargin{carg};carg=carg+1;
  nd=varargin{carg};carg=carg+1;
  ElemF=varargin{carg};carg=carg+1;
  RunOpt=struct('FieldDofs',[],'pro',fe_super('prop',ElemF), ...
   'node',fe_super('node',ElemF),'PropertyDof',[]);
  i4=[];
 %proid
  if ~isfield(model,'pl');   model.pl=zeros(0,1); end
  if ~isfield(model,'il');   model.il=zeros(0,1); end
  i1=unique(model.Elt(cEGI,RunOpt.pro(find(RunOpt.pro(1:2)))),'rows');
  for j1=1:size(i1,1) % loop on MatId ProId pairs
   if isempty(model.pl);pl=[];
   else; pl=model.pl(find(model.pl(:,1)==i1(j1,1)),:);
   end
   if isempty(pl);
      typ='m_elastic';st='SI';i2=1;
      sdtw('_nb','model.pl should containt MatId %i for BuildDof',i1(j1,1));
   else; [typ,st,i2]=fe_mat('typem',pl(1,2));
   end
   % default if undefined is a p_solid for backward compat
   il=model.il; 
   if ~isempty(il); il=il(find(il(:,1)==i1(j1,2)),:);end
   if isempty(il)&any(strcmp(ElemF,{'quad4','tria3','quadb','tria6', ...
    'beam1','beam1t'}))
     sdtw('_nb','model.il should containt ProId %i for BuildDof',i1(j1,2));
     RunOpt.pFcn='p_shell';RunOpt.pUnit='SI';RunOpt.pTyp=1;
   elseif isempty(il)
     RunOpt.pFcn='p_solid';RunOpt.pUnit='SI';RunOpt.pTyp=0;
     sdtw('_nb','model.il should containt ProId %i for BuildDof',i1(j1,2));
   else; [RunOpt.pFcn,RunOpt.pUnit,RunOpt.pTyp]=fe_mat('typep',il(1,2));
   end

   switch RunOpt.pFcn % property type
   case 'p_solid'
     RunOpt.FieldDofs=[1:3]; % this is the default
     if RunOpt.pTyp==3 % Surface elements for coupling
       if il(1,4)==4;  RunOpt.FieldDofs=[1:3 19]; % fluid structure coupling
       elseif il(1,4)==3;RunOpt.FieldDofs=[1:3]; % pressure
       else; sdtw('_nb','p_solid(3) type %i not supported',il(1,4));
       end
     elseif isequal(typ,'m_elastic')
      if i2==2; RunOpt.FieldDofs=[19];
      elseif i2==4|RunOpt.pTyp==2; RunOpt.FieldDofs=[1 2];  % Plane elements
      else % if any(i2)==[1 3]; 
        RunOpt.FieldDofs=[1:3];
      end
     elseif isequal(typ,'m_piezo');
       RunOpt.FieldDofs=[1 2 3 21];
     else; 
      try;  RunOpt.FieldDofs=feval(typ,'FieldDofs',RunOpt,i2);
      catch;
        sdtw('(p=%s and m=%s) is unknown',RunOpt.pFcn,typ);
      end
     end 
   % This calls a user p_*.m function. The call is supposed to fill in the
   % RunOpt.FieldDofs based on the contents of pl and il
   otherwise 
      try;
       RunOpt=feval(RunOpt.pFcn,'BuildDofOpt',RunOpt,model.pl,il,model);
      end
   end
  end

  % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  % Below here is a common to all element functions - - - - - - - -
  RunOpt.FieldDofs=RunOpt.FieldDofs(:);
  if isempty(nd); % build the DOF list for use in feutil('getdof')
   i2=unique(model.Elt(cEGI,RunOpt.node))*100;i2=i2(:);
   i3=RunOpt.FieldDofs(:)';
   i2=(i2(:,ones(size(i3)))+i3(ones(size(i2)),:))';i2=i2(:);
   out=[i2;round(RunOpt.PropertyDof*100)];
   out1=[];% return nodal and element DOFs
   evalin('caller',sprintf( ...
    'out2=out2+length(cEGI)*%i^2;', ...
    length(RunOpt.FieldDofs)*length(RunOpt.node)+length(RunOpt.PropertyDof)));
  else % build the DofPos for use in fe_mknl
   Case=evalin('caller','Case');
   DofPos=reshape(int32(full(nd(RunOpt.FieldDofs(:), ...
     model.Elt(cEGI,RunOpt.node)'))-1), ...
      length(RunOpt.FieldDofs)*length(RunOpt.node),length(cEGI));
   if ~isempty(RunOpt.PropertyDof); 
    i4=RunOpt.PropertyDof;i4=int32(full(nd(21,fix(i4)))-1);i4=i4(:);
    DofPos(end+[1:length(i4)],:)=i4(:,ones(size(DofPos,2),1));
   end
   out=DofPos;out1=[];
  end % Build DOF list or DofPos


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

% --------------------------------------------------------------------
