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

%BEAM1T 2 node beam with pretension. For non-linear cable statics and dynamics
%
%       This element has an internal state where each colum of 
%       Case.GroupInfo{5} gives the local basis and element tension
%        [bas(:) T]

%	Etienne Balmes, Jean-Michel Leclere
%       Copyright (c) 2001-2004 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.18 $  $Date: 2006/05/02 19:42:39 $'; return;
end

if ischar(CAM)

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

  [out,out1,out2]=beam1('integinfo',varargin{:});

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

           % State  ElementConstants        Reindex
   out='[gstate,Case.GroupInfo{jGroup,8}]=beam1t(''BasisAndConstants'',model,node,model.Elt(EGroup(jGroup)+1:EGroup(jGroup+1)-1,:),integ);';

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

node=varargin{2};elt=varargin{3}; 
state=zeros(11,size(elt,1));
NNode=sparse(node(:,1),1,1:size(node,1));

for jElt=1:size(elt,1) 
 x=node(NNode(elt(jElt,[1 2])),5:7);
 if size(elt,2)>6
  if (any(elt(jElt,6:7))|rem(elt(jElt,5),1))  % normal given
   x0 = x(1,:)+elt(jElt,5:7);
  else
   x0 = find(node(:,1)==elt(jElt,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(jElt,5)==0; x0=[1.5 1.5 1.5];  % if nothing
 else x0 = find(node(:,1)==elt(jElt,5));
  if isempty(x0); x0=[1.5 1.5 1.5];
  else x0 = node(x0,5:7); end
 end
 L = norm(x(2,:)-x(1,:));
 x=sp_util('basis',x(2,:)-x(1,:),x0-x(1,:))';
 state(1:9,jElt)=x(:);
 state(10,jElt)=L;

end
out=state; out1=integrules('beam1'); % constants
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% returns a normal map for first bending plane
 elseif  comstr(Cam,'map');

  model=varargin{1}; if isa(model,'v_handle');model=model.GetData;end
  model.Elt=feutil('setgroup beam1 name beam1t',model);
  [eltid,model.Elt]=feutil('eltidfix',model);
  Case=fe_mknl('initnocon',model);
  [EGroup,nGroup]=getegroup(model.Elt);
  out=struct('ID',[],'normal',[],'vertex',[],'n2',[]);
  NNode=sparse(model.Node(:,1),1,1:size(model.Node,1));
  for jGroup=1:nGroup
     ElemF= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
     cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
     if strcmp(ElemF,'beam1t');
      out.ID=[out.ID;eltid(cEGI)]; ind=[1:length(cEGI)];
      out.vertex(end+[1:length(cEGI)],1:3)= ...
       (model.Node(NNode(model.Elt(cEGI,1)),5:7)+ ...
         model.Node(NNode(model.Elt(cEGI,2)),5:7))/2;
      out.normal(end+ind,1:3)=Case.GroupInfo{jGroup,5}([2 5 8],:)';
      out.n2(end+ind,1:3)=Case.GroupInfo{jGroup,5}([3 6 9],:)';
     end
  end
  if sp_util('issdt')&nargout==0;cf=feplot;cf.model=model;fecom(cf,'showmap',out);end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'state');   % call for state update

   out='beam1t(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,gstate,FieldAtEltDof,int32(jElt));';

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'call')  % call for matrix assembly

   out='[k1,m1]=beam1t(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,gstate(:,jElt),Case.GroupInfo{jGroup,8});';

 elseif comstr(Cam,'rhscall') % call for load assembly

   out='be=beam1t(nodeE,elt(cEGI(jElt),:),pointers(:,jElt),integ,constit,gstate(:,jElt),Case.GroupInfo{jGroup,8},defe);';

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

  model=femesh('testbeam1 divide 50');
  model.Elt=feutil('set group 1 name  beam1t',model);

  model=fe_case(model, 'reset', ...
   'fixdof','pinned',[1+[.01;.02;.03];2+[.01;.02;.03]], ...
   'fixdof','2D',[.01; .02;.04;.06]);
  model.il=[112 fe_mat('p_beam','SI',1) [1e-12 1e-12 1e-12] .1];
  model.pl(1,5)=1000; % 1000 kg/m3

  if comstr(Cam,'eig')
    def=fe_eig(model);
    out=model;out1=def; return;
  end


  [Case,model.DOF]=fe_mknl('init',model);
  T=1e8; Case.GroupInfo{1,5}(11,:)=T;

  k=fe_mknl('assemble',model,Case,1);
  m=fe_mknl('assemble',model,Case,2);

  g=(m*sum(fe_c(Case.DOF,.03))');
  kd=ofact(k); q=kd\g;ofact('clear',kd)

  i1=feutil('findnode groupall',model);
  x=sort(model.Node(i1,5));
  i2=fe_c(Case.DOF,.03,'ind');
  x=model.Node(fix(Case.DOF(i2)),5);
  q(i2,2)=model.il(1,6)*model.pl(1,5)*1/T*(x.*(1-x))/2;

  figure(1);plot(x,q(i2,1)./q(i2,2))
  ylabel('Computed/true deflection');xlabel('x(m)');

  model.il=p_beam('dbval 112 circle 1e-2');

  [Case,model.DOF]=fe_mknl('init',model);
  T=1e3; Case.GroupInfo{1,5}(11,:)=T;
  k=fe_mknl('assemble',model,Case,1);
  f=fe_c(Case.DOF,28.03)'*100;
  kd=ofact(k); q=kd\f;ofact('clear',kd)
  r1=q(i2,1);

  model=fe_case(model,'fixdof','2D',[.01;.03;.04;.05]);
  [Case,model.DOF]=fe_mknl('init',model);
  T=1e3; Case.GroupInfo{1,5}(11,:)=T;
  k=fe_mknl('assemble',model,Case,1);
  f=fe_c(Case.DOF,28.02)'*100;
  kd=ofact(k); q=kd\f;ofact('clear',kd)
  r1(:,2)=q(i2,1);
  figure(102);h=plot(x,r1);set(h(2),'linestyle',':','marker','x');

 if 1==2 % arch test to correct normal computations - - - - - - -

 FEnode=[1 0 0 0  -1 -1 0];
 femesh('object mass 1'); 
 femesh('rev 20 o 0 0 0 180 -1.0 1.0 0')
 femesh plotel0
 FEel0(2:end,3:4)=1; FEel0(2:end,7)=1;
 model=femesh('model0');
 MAP=beam1t('map',model) 

 end% arch test to correct normal computations - - - - - - -

 else out=beam1(CAM);
 end

return
end % of standard calls with one input argument


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

typ=point(5);

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


if any(300==typ) % state update - - - - - - - - - - - - - - - - - -

 gstate=varargin{5};
 FieldAtEltDof=varargin{6}; 
 jElt=varargin{7};

 state=gstate(:,jElt);
 L=state(10); x=reshape(state(1:9),3,3);

 r1=constit(point(7)+1)*constit(point(7)+10)/L; % EA/L

 r1=(-x(1,:)*FieldAtEltDof(1:3)+x(1,:)*FieldAtEltDof(7:9))*r1;
 sp_util('setinput',gstate,r1,int32(11),jElt);
 out=[];
 return; 

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

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

state=varargin{5};rules=varargin{6};

L=state(10); x=reshape(state(1:9),3,3); Tpin=[];

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

   k = zeros(12,12);
   r1=constit(point(7)+1)*constit(point(7)+10)/L; %EA/L
   ind = [1 7 73 79];     k(ind) = [1 -1 -1 1]*r1; % 1 7
   % 4 10, torsionnal stiffness GJ/l
   r1=constit(point(7)+4)*constit(point(7)+7)/L; %GJ/L
   ind = [40 46 112 118]; k(ind) = [1 -1 -1 1]*r1; % 4 10

   wx=rules.Nr;wx([1 2 3 7 8 9])=wx([1 2 3 7 8 9])/L;
   wxx=rules.Nrr/L;wxx([1 2 3 7 8 9])=wxx([1 2 3 7 8 9])/L;
   r1=state(11)*L;                               % T  * L

   r2=constit(point(7)+1)*constit(point(7)+8)*L; % EI * L
   dk=zeros(4);for j1=1:size(rules.w,1)
    dk=dk+rules.w(j1,4)*r1*wx(j1,:)'*wx(j1,:)+  ... % Pre-tension stiffness
       rules.w(j1,4)*r2*wxx(j1,:)'*wxx(j1,:);       % bending stiffness
   end
   k([14 18 20 24 62 66 68 72 86 90 92 96 134 138 140 144])=dk;%[2 6 8 12];

   wx([1 2 3 7 8 9])=-wx([1 2 3 7 8 9]);  % theta_y = - dz/dx
   wxx([1 2 3 7 8 9])=-wxx([1 2 3 7 8 9]);

   r2=constit(point(7)+1)*constit(point(7)+9)*L; % EI * L
   dk=zeros(4);for j1=1:size(rules.w,1)
    dk=dk+rules.w(j1,4)*r1*wx(j1,:)'*wx(j1,:)+  ... % Pre-tension stiffness
       rules.w(j1,4)*r2*wxx(j1,:)'*wxx(j1,:);       % bending stiffness
   end
   k([27 29 33 35 51 53 57 59 99 101 105 107 123 125 129 131])=dk;%[3 5 9 11];

   % 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;
    i1=[i1(find(i1)) i2(find(i2))+6];
    if ~isempty(i1)
      i2=1:12;i2(i1)=0;i2=find(i2);
      Tpin=zeros(12); Tpin([i2(:);i1(:)],i2)= ...
        [eye(length(i2),length(i2));-k(i1,i1)\k(i1,i2)];
      k0=k;k=Tpin'*k*Tpin;
    end
   end

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

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


   rho=constit(point(7)+3);A=constit(point(7)+10);
   m = zeros(12,12);

    % consistent mass
    %ind=[1 7];m=zeros(12);m(ind,ind)=ones(2);sprintf('% i',find(m))
    r1=rho*A*L/3; % Rho * A * L
    m([1 7 73 79]) = [1 .5 .5 1]*r1;

    r1=rho*sum(constit(point(7)+[8:9]))*L/3; % Rho * (i1+i2) * L
    m([40 46 112 118]) = [1 .5 .5 1]*r1;

    L2=L^2; r1=constit(point(7)+8);
    if r1; ind = [2 6 8 12];
      r1=rho*A*L/210;  % rho I L /210
       m(ind,ind) = r1 * ...
	[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];
    end
    r1=constit(point(7)+9);
    if r1; ind = [3 5 9 11];
      r1=rho*A*L/210; % rho I L /210
      m(ind,ind) =  r1 * ...
	[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];
    end
    if ~isempty(Tpin); m=Tpin'*m*Tpin; end

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

    m=of_mk('xkx_trans',x,m); 
    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
   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}; % element coordinates as colums

  dxds=L;
  Nu=[1-rules.w(:,1) rules.w(:,1)];
  Nw=rules.N; Nw(:,[2 4])=Nw(:,[2 4])*L; % rotation DOFs


  r1=x(:,1)*x(:,1)';  r2=x(:,2)*x(:,2)';  r3=x(:,3)*x(:,3)';

  

  b1=zeros(12,1);
  for j1=1:size(Nu,1) % loop on integration points

    d=zeros(12,3);
    d(1:3,:)=r1*Nu(j1,1);% node 1 translation
    d(7:9,:)=r1*Nu(j1,2);% node 2 translation
    % Work associated with transverse displacement
    d(1:3,:)=d(1:3,:)+r2*Nw(j1,1); d(4:6,:)  =r2*Nw(j1,3);
    d(7:9,:)=d(7:9,:)+r2*Nw(j1,2); d(10:12,:)=r2*Nw(j1,4);

    d(1:3,:)=d(1:3,:)+r3*Nw(j1,1); d(4:6,:)  =r3*Nw(j1,3);
    d(7:9,:)=d(7:9,:)+r3*Nw(j1,2); d(10:12,:)=r3*Nw(j1,4);

    F = dire*Nu(j1,1:2)';
    b1=b1+d*(F*dxds*constit(point(7)+10)*rules.w(j1,4));
  end % loop on integration points

  out=b1;

end


