function [idof,k,m]=celas(node,elt,pl,il,opt,Case)

%CELAS	element function for scalar springs and penalized rigid links
%
%	As all element functions (see ELEM0), CELAS is called by FE_MK for
%	model assembly, FEPLOT for structural deformation visualization, ...
%
%	In an model description matrix a group of CELAS elements starts with a
%	header row [Inf  abs('celas') 0 ...] followed by element property rows
%       ELT following the format
%	    [n1 n2 DofID1 DofID2 ProID EltID  Kv Mv Cv Bv]
%         with
%	   n1,n2  node numbers of the connected nodes. Grounded springs are
%		  obtained by setting n1 or n2 to 0. 
%	   DofID  Identification of selected DOFs.
%		For RIGID LINKS, the first node defines the rigid body motion.
%		  DofID (>0) defines which DOFs of the slave are connected
%		  by the constraint. Thus [1 2 123 0 0 0 1e14] will only impose
%		  translations of node 2 are imposed by motion of node 1, while
%		  [1 2 123456 0 0 0 1e14] will also penalize the difference in
%		  rotations.
%		For SCALAR SPRINGS, DofID1 (negative) defines which DOFs of
%		  node 1 are connected to which of node 2. DofID2 can be used
%		  to specify different DOFs on the 2 nodes. For example
%		  [1 2 -123 231 0 0 1e14] connects DOFs 1.01 to 2.02, etc.
%	   ProID  Property identification number (see format below)
%	   Kv     Stiffness value used as a weighting  associated with the 
%		  constraint. If Kv is zero (or not given), the default value
%		  in the element property declaration is used. If this
%		  is still 0, Kv is set to sdtdef('kcelas').
%
%	Element property rows for CELAS elements take the form
%	   [ProID 4  KvDefault]
%
%	Warning : when seeking rigid rotations, be sure that the master node
%	  n1 has rotational DOFs. 
%
%	See also help rigid
%	         doc  rigid, bar1, eltfun, elem0

%	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.14 $  $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(node,'rhscall'); idof='';
 elseif comstr(node,'call')
   idof = ['[i1,k1,m1] = celas(nodeE, elt(cEGI(jElt),:),pl,il,[opt(1) jGroup jElt],Case);'];
 elseif comstr(node,'matcall')
   idof=celas('call');k=0;
 elseif  comstr(Cam,'groupinit');   idof = '';
 elseif  comstr(node,'node');    idof = [1 2];
 elseif  comstr(node,'prop');    idof = [0 5 6];
 elseif  comstr(node,'dof');    idof=[1+[1:6]/100 2+[1:6]/100];
 elseif  comstr(node,'line');    idof =[1 2];
 elseif  comstr(node,'patch');   idof= [1 2];
 elseif  comstr(node,'edge');   idof =[1 2];
 elseif  comstr(node,'face');    idof= [];
 elseif  comstr(node,'sci_face'); idof = [1 2];
 elseif  comstr(node,'parent');  idof = 'beam1';
 end
 return
elseif nargin==0; idof=[];return; end %of standard calls with 1 input argument

if nargin<6; NNode(node(:,1))=1:size(node,1);end

i2=abs(sprintf('%i',abs(elt(1,3))))-48; % DOF selected for node 1
i3=abs(sprintf('%i',abs(elt(1,4))))-48; % DOF selected for node 2

if ~isempty(il); i1 = find(il(:,1)==elt(1,5));
  if ~isempty(i1); elt(1,7)=il(i1,3); end
end
if size(elt,2)<9 elt(1,10)=0;end
if elt(1,7)==0; elt(1,7)=sdtdef('kcelas'); end
if isempty(i2); warning('No DOF selected'); end

idof=[];
if ~isempty(Case.cGL) 

  error('Local displacement coordinates not yet supported');
  in1=full(FEC.nd(round(idof*100)-100));
  T=FEC.cGL(in1,in1);
  k=T*k*T';

elseif elt(1,3)>0&size(node,1)>1  % penalized linearized rigid link
  
       x = node(:,1:3);x=x(2,:)-x(1,:);
       x = [0 x(3) -x(2) ;-x(3) 0 x(1);x(2) -x(1) 0];
       r = [eye(3,3) x;zeros(3,3) eye(3,3)];
       r = sparse([r(i2,:) -eye(length(i2),length(i2))]);
       k=zeros(12,12);m=k;i3=[1:6 i2+6];
       k(i3,i3)=(elt(1,7)*r'*r);
       m(i3,i3)=(elt(1,8)*r'*r);

else % standard spring connection

       if isequal(i3,0); i3=i2; end
       if length(i3)~=length(i2)
         error('non-matching numbers of DOFs for standard spring connection');
       end
       r = sparse([eye(length(i2),length(i2)) -eye(length(i2),length(i2))]);
       if any(elt(1,1:2)==0)
        idof =[elt(1,1)+i2(:)/100; elt(1,2)+i3(:)/100];
        i1=find(idof<1);r(:,i1)=0;%k(:,i1)=0;k(i1,:)=0;
       end
       k=zeros(12,12);m=k;i3=[i2 i3+6]; 
       k(i3,i3)=(elt(1,7)*r'*r);
       m(i3,i3)=(elt(1,8)*r'*r);

end

if opt(1)==2; k=m;
elseif opt(1)==5; % NL tangent stiffness
elseif opt(1)>2; idof=[]; k=[]; end
