function [idof,k,m]=cbush(node,elt,pl,il,opt,bas)

%CBUSH element function for the 2-node spring/damper element
%
%       The element property row of a cbush element takes the form
%
%       n1 n2 MatId ProId EltId x1 x2 x3 CID S OCID S1 S2 S3
%
%           OCID : coordinate system for offset definition
%           Si   : components of offset in OCID
%
%	See also help beam1, elem0
%	         doc  p_spring

%	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 comstr(node,'cvs')
 idof='$Revision: 1.17 $  $Date: 2006/01/27 14:28:50 $'; return;
end
% standard calls with one input argument
if ischar(node)
 [CAM,Cam]=comstr(node,1);
 if comstr(Cam,'integinfo')

   idof=[];k=[];m=[];

 elseif      comstr(Cam,'call')
   idof = ['[i1,k1,m1] = cbush(node(NNode(elt(cEGI(jElt),[1 2])),:),elt(cEGI(jElt),:),pl,il,[opt(1) jGroup jElt],bas);'];
 elseif  comstr(Cam,'groupinit');  idof = '';
 elseif  comstr(Cam,'node');  idof = [1 2];
 elseif  comstr(Cam,'prop');  idof = [3 4 5];
 elseif  comstr(Cam,'dof');   idof=[1+[1:6]/100 2+[1:6]/100];
 elseif  comstr(Cam,'line');  idof = [1 2];
 elseif  comstr(Cam,'face');   idof =[];
 elseif  comstr(Cam,'edge');   idof =[1 2];
 elseif  comstr(Cam,'patch');  idof = [1 2];
 elseif  comstr(Cam,'parent');   idof = 'beam1';
 elseif  comstr(Cam,'sci_face'); idof = [1 2 2];
 elseif  comstr(Cam,'test')
  warning('cbush : tests not implemented ');return
 end
return
end % of standard calls with one input argument

% element matrix assembly - - - - - - - - - - - - - - - - - - - - - - - - -

if size(node,1)>2
 i1=node(:,1);node=node([find(i1==elt(1));find(i1==elt(2))],:);
end

il=il(find(il(:,1)==elt(4)),:);

idof=[elt(1)+[1:6]/100 elt(2)+[1:6]/100]';
if  size(elt,2)<10; elt(10:14)=[.5 0 0 0 0]; end
% element coordinate system and S used
if size(elt,2)<11 | (size(elt,2)>10&elt(11)<1)
   S=[elt(10) 1-elt(10)]*node(:,5:7);
else % offset defined
   i2=find(bas(:,1)==elt(11));
   r1=basis([1 bas(i2,1) 0 0 elt(12:14)],bas(i2,:));
   S=r1(5:7);
end
if elt(9)==0;b=eye(3);
elseif elt(9)<1  % Coordinate system defined by nodes 
  if any(elt(7:8))|rem(elt(6),1)~=0
    b=sp_util('basis',diff(node(:,5:7)),elt(6:8)-node(1,:));
  else
    error('External G0 reference not implemented')
  end
else  % Externally defined coordinate system
  i2=find(bas(:,1)==elt(9));
  b=reshape(bas(i2,7:15),3,3)';
end

x=S-node(1,5:7);
tr1=[1 0 0   0 x(3) -x(2);0 1 0 -x(3) 0 x(1);0 0 1 x(2) -x(1) 0;
     0 0 0  1 0 0;0 0 0  0 1 0;0 0 0  0 0 1];
x=S-node(2,5:7);
tr2=[1 0 0   0 x(3) -x(2);0 1 0 -x(3) 0 x(1);0 0 1 x(2) -x(1) 0;
     0 0 0  1 0 0;0 0 0  0 1 0;0 0 0  0 0 1];

tr=[tr1 -tr2];
% standard cbush element

if any([0 1 4]==opt(1))	% stiffness or hysteretic damping
  k = diag(il(3:8)); if opt(1)==4; k=k*il(15);end
  k=of_mk('xkx_trans',b,k);
  k=tr'*k*tr;
else   k = zeros(12,12); end
if opt(1)==3  % viscous damping
  k = diag(il(9:14)); k=of_mk('xkx_trans',b,k); k=tr'*k*tr;
end

m=zeros(12,12);

if any([0 1 2 3 4 100]==opt(1))
elseif opt(1)==-1 % return CBUSH without any local basis reference

  be=b(:,3);[r1,i1]=max(abs(be));if be(i1)<0; be=-be;end
  be=basis(diff(node(:,5:7)),be);

 % n1 n2  MatId ProId EltId x1 x2 x3   CID S   OCID   S1 S2 S3
 idof=[elt(1:5)              be(:,2)'  -2 0.5  -1  (S-node(1,5:7))*be'];
 return;
 
else
 error(sprintf('Matrix type %i not supported by cbush',opt(1)))
end

