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

%P_HEAT property function of heat and other diffusion equations
%
%       Syntax : mat= p_heat('default') 
%
%       Supported ELEMENT property subtypes are
%       1 : 3D volume element for heat diffusion
%           [ProId fe_mat('p_heat','SI',1) CordM Integ DIM] 
%       2 : 3D surface element for heat exchange
%           [ProId fe_mat('p_heat','SI',2) CordM Integ DIM] 
%         Integ : is rule number in integrules
%         DIM is problem dimension 2 or 3 D
%
%       Supported MATERIAL property subtypes are
%       1 : 3D heat diffusion 
%           [MatId fe_mat('m_heat','SI',2) k rho C alpha]
%
%       See also help fe_mat
%                doc  p_shell, fem (handling materials section), pl, fe_mat

%	Etienne Balmes, Frederic Bourquin, A. Nassiopoulos
%       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.18 $  $Date: 2006/04/27 13:42:18 $'; return;
end
if nargin<1; help p_heat;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_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_heat('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_heat','SI',1) fe_mat('m_heat','SI',1) 400 0 0 1 ];
  % k rho C alpha 
  out.name='default heat';
  out.type='p_heat';
  out.unit='SI';

  out1='Heat properties';out2=carg;


  i1=strmatch(comstr(st,-27),comstr({out.name},-27));
  if ~isempty(i1) out=out(i1);
  elseif comstr(st,'shell') % piezo shell
   r1=comstr(st(6:end),-1,[1 1000 .01 0]);
   out=out(1);out.il(2+[1:length(r1)])=r1;
  end

% -------------------------------------------------------------------------
%[constit,integ]=p_heat('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;
pl=pl(find(pl(:,1)==ID(1)),:);
il=il(find(il(:,1)==ID(2)),:);
% Here one defines how a CONSTIT vector is specified - - - - -
if il(1,5)==3; K=diag(pl(1,3)*[1 1 1]);
else; K=diag(pl(1,3)*[1 1]);
end
       % -1 typ     rho*C           alpha    K  
constit=[-1;pl(1,2);pl(1,4)*pl(1,5);pl(1,6);K(:)];

[st,i1,i2]=fe_mat('typep',il(1,2));
ID=ID(:); ID(3)=ID(4); % only one DOF per node
ID(5:7)=[il(1,4);il(1,5);i2];  % Integrule, DIM, subtype
out=constit(:);     % Return CONSTIT
out1=int32(ID(:));  % Return integ

% -------------------------------------------------------------------------
% EltConst=p_solid('ConstSolid',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 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);opt.Nnode=size(opt.N,2);
 if size(integ,1)<9;rule=[1 opt.Nw]; % By default start at 1 use all points
 else; rule=integ(8:9,1);rule=rule(:)';
 end
 if ~any(rule); 
   rule=[1 opt.Nw];
   sdtw('integ is assumed to give standard integration rule in integ(5:6)');
 end
 opt.DofLabels={'T'}; 

 % 3D \int KgradT gradv - - - - - - - - - - - - - - - - - -
 if isequal(double(integ(6:7,1)),[3;1])
  % 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}={'T,x','T,y','T,z'};
  opt.ConstitTopology{1}=reshape(4+[1:9],3,3);
  % Kinetic energy
  opt.StrainDefinition{2}= [1 1 1 rule];
  opt.StrainLabels{2}={'T'}; opt.ConstitTopology{2}=3;
  % \int f v (see elem0 for format)
  opt.RhsDefinition=int32([101 0 1 0     0 0 -1    rule+[-1 0]]);
  out1=3;
 % 2D \int KgradT gradv - - - - - - - - - - - - - - - - - -
 elseif isequal(double(integ(6:7,1)),[2;1]) 
  % 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}={'T,x','T,y'};
  opt.ConstitTopology{1}=reshape(4+[1:4],2,2);
  % Kinetic energy
  opt.StrainDefinition{2}= [1 1 1 rule];
  opt.StrainLabels{2}={'T'}; opt.ConstitTopology{2}=3;
  % \int f v
  %opt.Rhs
  out1=2;
 % 3D \int_{\partial \omega} \alpha T v
 elseif isequal(double(integ(6:7,1)),[3;2]) 
  opt.StrainDefinition{1}= [1 1 1 rule];opt.StrainDefinition{2}=[];
  opt.StrainLabels{1}={'T'};
  opt.ConstitTopology{1}=eye(2)*4;
  % \int_{\partial \omega} (g+\alpha \theta_ext) v
  opt.RhsDefinition=int32([101 0 1 0     0 0 -1    rule+[-1 0]]);
  out1=23;
 % 2D \int_{\partial \omega} \alpha T v
 elseif isequal(double(integ(6:7,1)),[2;2]) 
  opt.StrainDefinition{1}= [1 1 1 rule];opt.StrainDefinition{2}=[];
  opt.StrainLabels{1}={'T'};
  opt.ConstitTopology{1}=4;
  % \int_{\partial \omega} g v
  %opt.Rhs
  out1=13;
 else; error('Not a supported case');
 end
 opt.VectMap=reshape(1:opt.Nnode*length(opt.DofLabels), ...
   length(opt.DofLabels),opt.Nnode)';
 opt=integrules('matrixrule',opt);
 opt.material='multi';

 if ~isempty(constit)
  % 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
% -------------------------------------------------------------------------
%RunOpt=p_heat('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=20; RunOpt.PropertyDof=[];
out=RunOpt;

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

 i1=varargin{carg}; 
 out1={};
 switch i1
 case {1,2} % [ProId fe_mat('p_heat','SI',1) CordM Integ DIM] 
 st={ ...
         'ProId'             0;
         'Type'              0;
         'CorM'              0;
         'Integ'              0;
         'DIM'              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='3D volume';
 case 2;  out='3D surface element';
 otherwise out='?';
 end

end

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