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

%BEAM1   element function for the 2-node 12 DOF beam element
%
%	As all element functions (see ELEM0), BEAM1 is called by FE_MK for
%	model assembly, FEPLOT for structural deformation visualization, ...
%
%	In model description matrices, a group of BEAM1 elements starts with a
%	header row [Inf  abs('beam1') 0 ...] followed by element property rows
%       ELT following the format
%	   [n1 n2 MatID ProID  (vx vy vz EltID p1 p2 x1 y1 z1 x2 y2 z2)]
%	   [n1 n2 MatID ProID  (nr 0  0  EltID p1 p2 x1 y1 z1 x2 y2 z2)]
%         with required values
%	   n1,n2  node numbers of the tips of the beam element
%	   MatID  material property identification number
%	   ProID  section property identification number
%	   NR     number of node not in the beam direction defining the first
%		  bending plane (default node is 1.5 1.5 1.5) or the alternate
%		  format
%	   vx,vy,vz components of a vector in plane 1 (not collinear to the
%		  beam axis)
%	  and optional values
%	   EltID  element identifier (optional, default 0)
%	   p1,p2  pin flags (list of released DOFs, e.g. 456 is pinned)
%	   x1, ... components of offset vector at node 1
%	   x2, ... components of offset vector at node 2
%
%       PL material property matrix. Isotropic materials [MatID 1  E  nu rho G]
%	   (See FE_MAT) are the only supported.
%
%       IL element property matrix. Rows associated to BEAM1 elements follow
%	   the format
%	    [ProID Type  J  I1  I2  A  (k1 k2)]
%          with
%	     Type : 1 for standard beam definition (no other type supported)
%	     J    : torsional stiffness parameter (often different from polar
%		    moment of inertia I1+I2)
%	     I1   : moment of inertia for bending plane 1
%		    which contains the reference node NR or given by direction
%		    vx,vy,vz.
%		    For a case with a beam along x and plane 1 the xy plane
%		    I1 is equal to Iz (integral over section of rho*y^2).
%	     I2   : moment of inertia for bending in the second plane which
%		    contains the beam and is orthogonal to plane 1
%            A    : section area
%	     k1   : (optional) shear factor for direction 1 (when not 0 a
%		     Timoshenko beam element is used)
%	     k2   : (optional) shear factor for direction 2
%
%	See also help bar1, p_beam, elem0
%	         doc  bar1, eltfun, elem0

%	Etienne Balmes
%       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.36 $  $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')

  [pe,ie]=fe_mat(1,varargin{1},varargin{2},varargin{3});
  pe=pe(:)';pe(6)=0;
  constit = [pe(1:6) ie]; %[E nu Rho G eta alpha T0  J I1 I2 A k1 k2 lump]
  iopt=varargin{1};%matid proid
  out1=iopt(:); out=constit(:); 
  out2=[];  % element map

 elseif comstr(Cam,'matcall');  out=beam1('call'); out1=0; % CallSymFlag

 elseif comstr(Cam,'call')  % call for matrix assembly
   out='[k1,m1]=beam1(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,elmap,node);';
 elseif comstr(Cam,'rhscall') % call for load assembly
   out='be=beam1(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,elmap,node,defe);';

 elseif  comstr(Cam,'state');   out = ''; % no state evaluation yet
 elseif  comstr(Cam,'groupinit');   out = '';
 elseif  comstr(Cam,'node');   out = [1 2];
 elseif  comstr(Cam,'prop');   out = [3 4 8];
 elseif  comstr(Cam,'dof');    out=[1+[1:6]/100 2+[1:6]/100];
 elseif  comstr(Cam,'edge');   out =[1 2];
 elseif  comstr(Cam,'face');   out =[];
 elseif  comstr(Cam,'sci_face'); out = [1 2 2];  
 elseif  comstr(Cam,'line');   out =[1 2];
 elseif  comstr(Cam,'patch');  out= [1 2];
 elseif  comstr(Cam,'parent'); out = 'beam1';

 % Basic matrix test of the element - - - - - - - - - - - - - - - - - -
 elseif  comstr(Cam,'testmat'); [CAM,Cam] = comstr(CAM,5);

   model=femesh('testbeam1');
   [constit,integ,elmap]=beam1('integinfo',[100;110],model.pl,model.il);
   [k,m]=beam1(model.Node,model.Elt(2,:),[0 0 0 0 0 0 0 0 0],int32(integ),constit,elmap,model.Node);
   model.Elt(2,5:7)=model.Node(3,5:7)-model.Node(1,5:7);
  [k1,m1]=beam1(model.Node,model.Elt(2,:),[0 0 0 0 0 0 0 0 0],int32(integ),constit,elmap,model.Node);
   if norm(k-k1)|norm(m-m1); error('Problem with orientation');end

   out=stack_cell(m,k);
   if nargout==0;disp('TestMat passed');end

 elseif  comstr(Cam,'test');   [CAM,Cam] = comstr(Cam,5);

   [out,out1]=femesh(strcat(['teststruct beam1' Cam]));

 % Basic tests of the element - - - - - - - - - - - - - - - - - - - - - -
 else sdtw('''%s'' unknown',CAM);  end

return
end % of standard calls with one input argument

% element matrix assembly - - - - - - - - - - - - - - - - - - - - - - - - - -
nodeE=CAM; 
elt=varargin{1}; 
point=double(varargin{2});
integ=varargin{3};
%constit=varargin{4};
constit=varargin{4}; 
elmap=varargin{5};
node=varargin{6};
typ=point(5);

if isa(elmap,'int32'); elmap = double(elmap); end

%if nargin<6 NNode(node(:,1))=1:size(node,1); end
%idof = [(1:6)/100+elt(1,1) (1:6)/100+elt(1,2)];
if size(nodeE,2)>4 
 node=nodeE; NNode(node(:,1))=1:size(node,1);
 nodeE = nodeE(NNode(elt(1,1:2)),[5:7 1]);
end
x=nodeE(:,1:3);

if size(elt,2)>6
  if (any(elt(6:7))|rem(elt(5),1))  % normal given
   x0 = x(1,:)+elt(5:7);
  else
   x0 = find(node(:,1)==elt(1,5));
   if isempty(x0); x0=[1.5 1.5 1.5];
   else x0 = node(x0,5:7); 
   end
  end

elseif size(elt,2)<5 | elt(1,5)==0; x0=[1.5 1.5 1.5];  % default if nothing
else x0 = find(node(:,1)==elt(1,5));
  if isempty(x0); x0=[1.5 1.5 1.5];
  else x0 = node(x0,5:7); end
end

% off-set if any - - - - - - - - - - - - - - - - - - - - - - - - -
tr = eye(12,12);off=0;

if size(elt,2)>12
if any(elt(11:13))
  tr([51 61 38])= elt(11:13);%[3 5;1 6;2 4];ans(:,1)+12*(ans(:,2)-1)
  tr([62 39 49])=-elt(11:13);%[2 6;3 4;1 5];ans(:,1)+12*(ans(:,2)-1);ans'
  off=1;
end
end
if size(elt,2)>15 
if any(elt(14:16))
  tr([129 139 116])= elt(14:16);%[3 5;1 6;2 4]+6;ans(:,1)+12*(ans(:,2)-1)
  tr([140 117 127])=-elt(14:16);%[2 6;3 4;1 5]+6;ans(:,1)+12*(ans(:,2)-1);ans'
  off=1;
end
end

% basis function - - - - - - - - - - - - - - - - - - - - - - - - - - - -

l = norm(x(2,:)-x(1,:)); l2 = l^2;
if l==0; error('beam with zero length'); end

x=sp_util('basis',x(2,:)-x(1,:),x0-x(2,:))';


%x(1,:) = x(2,:)-x(1,:);x(2,:) = x0-x(2,:);
%x(1,:) = x(1,:)/norm(x(1,:));
%x(2,:) = x(2,:) - (x(1,:)*x(2,:)')*x(1,:);
%x(2,:) = x(2,:)/norm(x(2,:));
%x(3,:) = [ -x(2,2)*x(1,3)+x(1,2)*x(2,3) ...
%            -x(2,3)*x(1,1)+x(1,3)*x(2,1) ...
%            -x(2,1)*x(1,2)+x(1,1)*x(2,2)];

%[pe,ie]=fe_mat(1,elt(3:4),pl,il);
pe=constit(point(7)+[1:6]);pe=pe(:)';
ie=constit(point(7)+[7:12]);ie=ie(:)';

if ~any([0 1 2 3 100]==typ)
 warning(sprintf('Matrix type %i not supported by beam1',typ))
end
Tpin=[];

% stiffness - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

if any([0 1]==typ) | (size(elt,2)>9&any(elt(9:10)))

   k = zeros(12,12);
   ind = [1 7 73 79];     k(ind) = [1 -1 -1 1]*pe(1,1)*ie(1,4)/l; % 1 7
   % 4 10, torsionnal stiffness GJ/l
   ind = [40 46 112 118]; k(ind) = [1 -1 -1 1]*pe(1,4)*ie(1,1)/l; % 4 10
   ind = [14 18 20 24 62 66 68 72 86 90 92 96 134 138 140 144]; %[2 6 8 12];
   if ie(5)==0 % bernoulli
      k(ind) = (pe(1,1)*ie(1,2)/l)* ...
	[12/l2    6/l   -12/l2   6/l     6/l   4  -6/l  2 ...
	 -12/l2  -6/l    12/l2  -6/l     6/l   2  -6/l  4];
   else % timoshenko with reduced integration
      k(ind) = ...
        (pe(1,1)*ie(1,2)/l)*[0 0 0 0 0 1 0 -1 0 0 0 0 0 -1 0 1] + ... %flexion
        (pe(4)*ie(4)*ie(5)/l)*[1 l/2 -1 l/2   l/2 l2/4 -l/2 l2/4 ...
                           -1 -l/2 1 -l/2   l/2 l2/4 -l/2 l2/4];
   end      

   ind = [3 5 9 11];
   if ie(6)==0
     k(ind,ind) = (pe(1,1)*ie(1,3)/l)* ...
	[12/l2   -6/l   -12/l2  -6/l;-6/l   4   6/l  2;
	 -12/l2   6/l    12/l2   6/l;-6/l   2   6/l  4];
   else
      k(ind,ind) = ...
        pe(1,1)*ie(1,3)/l*[0 0 0 0;0 1 0 -1;0 0 0 0;0 -1 0 1] + ... % flexion
        pe(4)*ie(4)*ie(6)/l*[1 -l/2 -1 -l/2;-l/2 l2/4 l/2 l2/4;         % shear
                           -1 l/2 1 l/2;-l/2 l2/4 l/2 l2/4];
   end

   % pin flag handling (condense DOF a put zero stiffness)
   if size(elt,2)>9
    i1=abs(sprintf('%i',elt(1,9)))-48; i2=abs(sprintf('%i',elt(1,10)))-48;
    i2=[i1(find(i1)) i2(find(i2))+6];
    if ~isempty(i2)
      i1=i2(find(i2)); i2=1:12;i2(i1)=0;i2=find(i2);
      Tpin([i2(:);i1(:)],1:length(i2))=[eye(length(i2),length(i2));-k(i1,i1)\k(i1,i2)];
      k1=zeros(12,12); k1(i2,i2)=Tpin'*k*Tpin; k=k1;
    end
   end

   k=of_mk('xkx_trans',x,k);  if off; k = tr'*k*tr;end % coordinate transformation
   out=k; out1=[];
end

if any([0 2]==typ) % mass - - - - - - - - - - - - - - - - - - - - -

   m = zeros(12,12);
   ind = [1 7];    m(ind,ind) = [1 .5;.5 1]/3 * ie(4);
   ind = [4 10];   m(ind,ind) = [1 .5;.5 1]/3*(ie(2)+ie(3)); % correct 23/4/98
   if ie(2); ind = [2 6 8 12]; m(ind,ind) =  ie(4) * ...
	[78    11*l   27   -6.5*l; 11*l    2*l2    6.5*l -1.5*l2
	 27    6.5*l  78   -11*l ; -6.5*l -1.5*l2 -11*l   2*l2]/210;
   end
   if ie(3); ind = [3 5 9 11]; m(ind,ind) =  ie(4) * ...
	[78   -11*l   27    6.5*l;-11*l    2*l2   -6.5*l -1.5*l2
	 27   -6.5*l  78    11*l ; 6.5*l  -1.5*l2  11*l   2*l2]/210;
   end
   m = (pe(3) * l) * m;
   if ~isempty(Tpin);  m1=zeros(12,12); m1(i2,i2)=Tpin'*m*Tpin; m=m1; end

   % Lumped mass (sum over elements)
   if length(constit)>=point(7)+13 &constit(point(7)+13)==1

    m=of_mk('xkx_trans',x,m); if off; m = tr'*m*tr;end % coordinate transform

    r1=sum(m);
    %r2=[0 0 0 1 0 0   0 0 0 1 0 0;
    %    0 0 0 0 1 0  0 0 -L 0 1 0;0 0 0 0 0 1   0 L 0 0 0 1];
    r2=[zeros(3);x';x'*[0 0 0;0 0 l;0 -l 0];x'];
    r3=sum(r2'*m*r2)/2;if any(r3)<0 r3=diag(r2'*m*r2)'/2;end
    r1(4:6)=r3; r1(10:12)=r3; 
    m=diag(r1);
   else
     m=of_mk('xkx_trans',x,m); % coordinate transformation
     if off; m = tr'*m*tr;end % coordinate transform   
   end
   if typ==2; out=m; out1=[];  
   else out1=m; end
end

if typ==3; out=zeros(12,12); out1=[]; end 

if any(100==typ) % volumic load - - - - - - - - - - - - - - - - - -

  dire=varargin{7};x=x'; % element coordinates as colums

  if length(ie)<8; ie(8)=0; end

  FEw = [-sqrt(3/5) 0 0 .555555555555556;
                   0 0 0 .888888888888889;
           sqrt(3/5) 0 0 .555555555555556];
  %FEw=linspace(-1,1,50)';
  xi = [-1 0 0;1 0 0];
  FEna = [(1+FEw(:,1)*xi(:,1)')/2 ...         % linear shape functions
         .25*[2-3*FEw(:,1)+FEw(:,1).^3 ...  % cubic shape functions
               (1-FEw(:,1)-FEw(:,1).^2+FEw(:,1).^3) ...
               2+3*FEw(:,1)-FEw(:,1).^3 ...
               (-1-FEw(:,1)+FEw(:,1).^2+FEw(:,1).^3)]];
  dxds = ones(size(FEw,1),1)*xi(:,1)'/2*[0;l];
  FEna(:,[4 6])=FEna(:,[4 6])./dxds(:,[1 1]); % rotation dofs

  %plot(FEw,FEna)
  r1=x(:,1)*x(:,1)';  r2=x(:,2)*x(:,2)';  r3=x(:,3)*x(:,3)';

  b1=zeros(12,1);
  for j1=1:size(FEw,1) % loop on integration points
    d=zeros(12,3);
    d(1:3,:)=r1*FEna(j1,1);% node 1 translation
    d(7:9,:)=r1*FEna(j1,2);% node 2 translation
    if ie(5)==0 
       d(1:3,:)=d(1:3,:)+r2*FEna(j1,3); d(4:6,:)  =r2*FEna(j1,4);
       d(7:9,:)=d(7:9,:)+r2*FEna(j1,5); d(10:12,:)=r2*FEna(j1,6);
    else 
       d(1:3,:)=d(1:3,:)+r2*FEna(j1,1); 
       d(7:9,:)=d(7:9,:)+r2*FEna(j1,2);
    end
    if ie(6)==0 
       d(1:3,:)=d(1:3,:)+r3*FEna(j1,3); d(4:6,:)  =r3*FEna(j1,4);
       d(7:9,:)=d(7:9,:)+r3*FEna(j1,5); d(10:12,:)=r3*FEna(j1,6);
    else 
       d(1:3,:)=d(1:3,:)+r3*FEna(j1,1); 
       d(7:9,:)=d(7:9,:)+r3*FEna(j1,2);
    end
    F = dire*FEna(j1,1:2)';
    b1=b1+d*(F*dxds(j1)*ie(4)*FEw(j1,4));
  end % loop on integration points

  out=b1;

end
