function [out,out1,out2]=dktp(CAM,varargin);

%DKTP element function for the 3-node 9-DOF 
%    triangular 3D element 
%
%	In an model description matrix a group of DKTP elements starts with a
%	header row [Inf  abs('dktp') 0 ...] followed by element property rows
%       ELT following the format
%	    [n1 n2 n3 MatID ProID EltID]
%         with
%	   n1 ... n3  identification numbers for the element nodes
%	   MatID  material property identification number
%	   ProID  element property identification number (only used by upcom)
%	   EltID  optional element identifier
%
%     PL material property rows are either isotropic materials
%           [MatId type E nu rho G eta alpha T0]
%           type=fe_mat('m_elastic',1,1)
%      or 2-D anisotropic materials
%           [MatId type E11 E12 E22 E13 E23 E33 rho eta a1 a2 a3]
%           type = fe_mat('m_elastic',1,4)
%
%     IL Element property rows for 2-D elements follow the format
%      [ProId Type  f   d 0 h]
%     with 
%       Type = fe_mat('p_solid',1,2)
%       f    : not used        	    
%       d    : not used
%       h    : thickness
%
%     Standard tests available with dktp('testeig') (mat,eig,load) http://www.mathworks.fr/store/default.do
%
%	See also help t3p, ...
% 	         doc  eltfun, elem0

%	E. Balmes, J. Leclere, H. Assime
%       Copyright (c) 2001-2003 by INRIA and SDTools,All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license

if comstr(CAM,'cvs')
 out='$Revision: 1.15 $  $Date: 2006/01/27 14:28:50 $'; return;
end

% standard calls with one input argument
if ischar(CAM)

 [CAM,Cam]=comstr(CAM,1);
 % Build Constit, Integ, and Elmap for later integration
 if comstr(Cam,'integinfo')

if 1==2
  [constit,iopt] = ...
     fe_mat('plil of_mk 2d',varargin{1},varargin{2},varargin{3}); %solid
end
  
  pl=varargin{2}; pl=pl(find(pl(:,1)==varargin{1}(1)),:);
  if length(pl)<7; pl(7)=0.; end
  il=varargin{3}; il=il(find(pl(:,1)==varargin{1}(1)),:);
  constit = [pl(5) il(6) pl([3 4]) il(6)]; % constit =[rho h E nu h]

  iopt(3:6)=[9 3 0 0];
  out1=iopt(:); out=constit(:); 

  if isempty(constit); error('Material properties are not defined');end
  % element map
  i1=find(triu(ones(9,9))); r1=zeros(9,9);r1(i1)=1:length(i1);
  out2=r1+tril(r1',-1);

 elseif comstr(Cam,'matcall');  out='mat_of';out1=1; % SymFlag
 elseif comstr(Cam,'call')
   out='[k1,m1]=dktp(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,elmap);';

 elseif comstr(Cam,'rhscall') % call for load assembly
   out='be=dktp(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,gstate,defe,EltConst,InfoAtNode);';

  elseif  comstr(Cam,'groupinit');   out = '';
  elseif  comstr(Cam,'node');   out = [1 2 3];
  elseif  comstr(Cam,'prop');   out = [4 5 6];
  elseif  comstr(Cam,'dof')    
   out =[1:3];k=[3 4 5]'/100;
   out=out(ones(3,1),:)+k(:,ones(3,1));out=out(:);
  elseif  comstr(Cam,'line');   out = [1 2 3 1];
  elseif  comstr(Cam,'patch');  out = [1 2 3];
  elseif  comstr(Cam,'edge');   out = [1 2; 2 3; 3 1];
  elseif  comstr(Cam,'face');   out = [1 2 3];
  elseif comstr(Cam,'sci_face'); out = [1 2 3];
  elseif  comstr(Cam,'parent'); out = 'tria3'; 

 % Basic matrix test of the element - - - - - - - - - - - - - - - - - -
  elseif  comstr(Cam,'testmat'); [CAM,Cam] = comstr(CAM,5);
   model=femesh('testdktp');
   [constit,iopt,elmap]=dktp('integinfo',[100;110],model.pl,model.il);
   [k,m]=dktp(model.Node,model.Elt(2,:),[45 45 0 0 0 0 0 0 0],int32(iopt),constit,elmap);

 iopt(5)=4; constit(2)=.5;
   k1=dktp(model.Node,model.Elt(2,:),[45 45 0 0 0 0 0 0 0],int32(iopt),constit,elmap);

   out=stack_cell(k,m);
   disp('TestMat passed');
  elseif  comstr(Cam,'test');   [CAM,Cam] = comstr(Cam,5);

   [out,out1]=femesh(strcat(['teststruct dktp' Cam]));
   % Basic tests of the element - - - - - - - - - - - - - - - - - - - - - -

  end
  return
end % of standard calls with one input argument

% element matrix assembly - - - - - - - - - - - - - - - - - - - - - - - - - -
node=CAM; 
elt=varargin{1}; 
point=varargin{2};
integ=varargin{3};
constit=varargin{4};
elmap=varargin{5};
if isa(elmap,'int32'); elmap = double(elmap); end


if size(node,2)~=4
 NNode(node(:,1))=1:size(node,1);
 node=node(NNode(elt(1,1:3)),[5:7 1]);
end

typ=point(5);
if (typ==0) 
  [k1,m1]=of_mk('dktp',int32(point),integ,constit,node);
  k=reshape(k1(elmap),size(elmap,1),size(elmap,2)); m=diag(m1(1:9));
  out=k; out1=m;
elseif typ==100
  point(5)=100; point(25)=0;
  out=of_mk('dktp',int32(point),varargin{3:4},node,varargin{[5 6 7 8]});

else

   k1=of_mk('dktp',int32(point),integ,constit,node);
   k=reshape(k1(elmap),size(elmap,1),size(elmap,2));m=[];
   out=k; out1=[];

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




