function Fe = readFEMfile(filename)
%  Fe = readFEMfile(filename)
%  Reads a file with a particular format that contains a description of 
%  an FEM project. Total.FEM is at the time of writing this software
% the only viable file. 
  
  fid = fopen(filename,'r');
  if fid<0, 
    disp([filename ,' could not be opened']); 
    return;
  end;
  
  Fe=0;
% looking for 'NNodes', 'NEleme'  'NGauss' 'NFiber' 'Muscle' 'NDirec' 'EOFile'

eof = 0;

cr=sprintf('\n');

while eof == 0,
  validstring=0;
  while validstring==0,
    cc = fscanf(fid,'%c',1);
    switch cc,
     case cr,
      disp('CR encountered');
      
     case '#',
      disp('Comment line');
      while cc~=cr,
	cc=fscanf(fid,'%c',1);
      end;
     case 'N',
      str = fscanf(fid,'%5s',1);
      str = [cc,str];
      validstring=1;
     case 'E',
      str = fscanf(fid,'%5s',1);
      str = [cc,str];
      validstring=1;
    end;
  end;
  
  disp(['Encountered String: ''',str,'''']);
  
  switch str,
   case 'NNodes',
    disp('Nodes found');
    [Fe.xyz,Fe.constr] = read_nodecoords(fid); 
    
   case 'NEleme',
    Fe.elem = read_elements(fid);
    
   case 'NGauss',
    ignore_gausspoints(fid);  % since you can compute them otherwise.
    
   case 'NFiber',
    nfibers=fscanf(fid,'%i',1);
    disp(sprintf('There are %i fiber types in the model\n',nfibers));
    for kfib = 1:nfibers,
      [musclename, indexfield, directionfield, densityfield] = read_a_fiber(fid);
      fprintf('Muscle name: %s  with %i Gauss points affected\n', ...
	      musclename,length(indexfield));
      Fe.musclenames{kfib}=musclename;
      Fe.fiberindx{kfib}=indexfield;
      Fe.fiberdir{kfib}=directionfield;
      Fe.fiberdensity{kfib}=densityfield;
    end;
    
    
    
   case 'EOFile',
    eof=1;
   
   otherwise,
    disp('File corrupted');
    eof=1;
    
  end;
end;

fclose(fid);
return;

function [XYZ,constraints]= read_nodecoords(fid); 
    nnodes = fscanf(fid,'%i',1);
    
    if nnodes>100000,
      disp('too many nodes - something is funny I bail out');
      return;
    end;
    
    for k=1:nnodes,
      no = fscanf(fid,'%i',1);
      xyz = fscanf(fid,'%f',3);
      con = fscanf(fid,'%i',3);
      XYZ(k,1:3)=xyz(:)';
      constraints(k,1:3)=con(:)';
    end;
    disp('Done reading node coordinates');
    return;
    
      
function Elements=read_elements(fid)
  nelem = fscanf(fid,'%i',1);
  if nelem>100000,
    disp('too many elements - something is funny . I bail out');
    return;
  end;
  
  for ne=1:nelem,
    no = fscanf(fid,'%i',1);
    
    ele=fscanf(fid,'%i',27);
    
    Elements(ne,:) = ele(:)';
  end;
  disp('Done reading elements');
  
  return;
  
function ignore_gausspoints(fid)
  ngauss=fscanf(fid,'%i',1);
  
  disp(sprintf('%i Gauss points are going to be skipped\n',ngauss));
  
  for k=1:ngauss,
    gg=fscanf(fid,'%f',3);
  end;
  
  disp('done reading and ignoring Gauss points');
  return;
  
function [name, indx, direct, dense] = read_a_fiber(fid)
  musc = fscanf(fid,'%6s',1);
  if ~isequal(musc,'Muscle'),
    disp('Oh-oh, file corrupted. Expecting the word ''Muscle''');
    indx=-1;
    return;
  end;
  
  musc = fscanf(fid,'%16s',1);
  str=splitaname(musc);
  name = str{1};
  
  musc = fscanf(fid,'%6s',1);    % expecting 'NDirec'
  
  if ~isequal(musc,'NDirec'),
    disp('Oh-oh, file corrupted. Expecting the word ''NDirec''');
    indx=-1;
    return;
  end;
  
  nfib = fscanf(fid,'%i',1);
  
  for k=1:nfib,
    no = fscanf(fid,'%i',1);
    xyzd = fscanf(fid,'%f',4);
    indx(k)=no;
    direct(k,1:3)=xyzd(1:3)';
    dense(k)=xyzd(4);
  end;
  
  return;