% FEM
% ------ generate element matrices, assemble in global matrices 
nd = n_nodes*n_dof;
clear K M C
K = zeros(nd); M = zeros(nd); C = zeros(nd); 

eFile = fopen ('Elements.out','w');
s_output (eFile,'Elements')
pFile = fopen ('Elem_Properties.out','w');
s_output (pFile,'Properties')

nr_elem = 0;

loc = 'Stiff elements';
for ie = eStiff
  eo(ie,:) = eo_beam (ie,Ex,Ey,Ez);
  [ke,me] = beam3d (Ex(ie,:),Ey(ie,:),Ez(ie,:),eo(ie,:),Ep_stiff);
  K = assem(Edof(ie,:),K,ke);
  M = assem(Edof(ie,:),M,me);
  nr_elem = nr_elem+1;
  m_elem(ie) = L_elem(ie) * Ep_stiff(7);
  fprintf (eFile,' %3d  %4d %4d  %7.4f %7.4f %7.4f %7.3f    %s\n',...
     ie,Elem(ie,:),eo(ie,:),m_elem(ie),loc);
end
fprintf (pFile,'%15s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n',...
         loc,Ep_stiff);
           
loc = 'Flex beam elements';
for ie = eFlex
  eo(ie,:) = [0 0 1];  
  [ke,me] = beam3d (Ex(ie,:),Ey(ie,:),Ez(ie,:),[0 0 1],Ep_lame);
  K = assem(Edof(ie,:),K,ke);
  M = assem(Edof(ie,:),M,me);
  nr_elem = nr_elem+1;
  m_elem(ie) = L_elem(ie) * Ep_lame(7);
  fprintf (eFile,' %3d  %4d %4d  %7.4f %7.4f %7.4f %7.3f    %s\n',...
     ie,Elem(ie,:),eo(ie,:),m_elem(ie),loc);
end
fprintf (pFile,'%15s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n',...
   loc,Ep_lame);

loc = 'Mass elements';
for ie = eMass
  eo(ie,:) = [0 0 1];
  [ke,me] = beam3d (Ex(ie,:),Ey(ie,:),Ez(ie,:),eo(ie,:),Ep_block);
  K = assem(Edof(ie,:),K,ke);
  M = assem(Edof(ie,:),M,me);
  nr_elem = nr_elem+1;
  m_elem(ie) = L_elem(ie) * Ep_block(7);
  fprintf (eFile,' %3d  %4d %4d  %7.4f %7.4f %7.4f %7.3f    %s\n',...
     ie,Elem(ie,:),eo(ie,:),m_elem(ie),loc);
end
fprintf (pFile,'%15s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n',...
   loc,Ep_block);

loc = 'Bar elements';
for ie = eAct
  eo(ie,:) = [0 0 1];
  [ke,me] = beam3d (Ex(ie,:),Ey(ie,:),Ez(ie,:),eo(ie,:),Ep_bar);
  K = assem(Edof(ie,:),K,ke);
  M = assem(Edof(ie,:),M,me);
  nr_elem = nr_elem+1;
  m_elem(ie) = L_elem(ie) * Ep_bar(3);
  fprintf (eFile,' %3d  %4d %4d  %7.4f %7.4f %7.4f %7.3f    %s\n',...
     ie,Elem(ie,:),eo(ie,:),m_elem(ie),loc);
end
fprintf (pFile,'%15s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n',...
   loc,Ep_bar);

loc = 'Flex hinge elements';
for ie = eH
  eo(ie,:) = [0 0 1];
  [ke,me] = beam3d (Ex(ie,:),Ey(ie,:),Ez(ie,:),eo(ie,:),Ep_H);
  K = assem(Edof(ie,:),K,ke);
  M = assem(Edof(ie,:),M,me);
  nr_elem = nr_elem+1;
  m_elem(ie) = L_elem(ie) * Ep_H(3);
  fprintf (eFile,' %3d  %4d %4d  %7.4f %7.4f %7.4f %7.3f    %s\n',...
     ie,Elem(ie,:),eo(ie,:),m_elem(ie),loc);
end
fprintf (pFile,'%15s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n',...
   loc,Ep_H);

Check = nr_elem
m_elem = m_elem';
fclose (eFile);
fclose (pFile);