function _mesh = load_gmsh(filename, init_mem)

if (~isdef('init_mem','local')) then
  init_mem = %F;
end

// Reads a mesh in msh format, version 1
_mesh     = [];    
_mesh.MIN = zeros(3, 1);
_mesh.MAX = zeros(3, 1);
fid = mopen(filename, 'r');
while %T
  endoffile = %F;
  while %T
    tline = mgetl(fid,1);
    if meof(fid) then endoffile=%T; break; end
    if part(tline,1) == '$' then break; end  
  end
  if endoffile then break; end
  if part(tline,2) == 'N' & part(tline,3) == 'O' then
    printf('reading nodes\n')
    _mesh.nbNod = mfscanf(1, fid, '%d');
    _mesh.POS   = zeros(_mesh.nbNod, 3);
    for I=1:_mesh.nbNod
      iNod      = mfscanf(1, fid, '%d');
      X         = mfscanf(3, fid, '%g');
      IDS(iNod) = I;
      if (I == 1) then
	       _mesh.MIN = X;
	       _mesh.MAX = X;
      else
	       if(_mesh.MAX(1) < X(1)) then _mesh.MAX(1) = X(1); end
	       if(_mesh.MAX(2) < X(2)) then _mesh.MAX(2) = X(2); end
	       if(_mesh.MAX(3) < X(3)) then _mesh.MAX(3) = X(3); end
	       if(_mesh.MIN(1) > X(1)) then _mesh.MIN(1) = X(1); end
	       if(_mesh.MIN(2) > X(2)) then _mesh.MIN(2) = X(2); end
	       if(_mesh.MIN(3) > X(3)) then _mesh.MIN(3) = X(3); end
      end
      _mesh.POS(I, 1) = X(1);
      _mesh.POS(I, 2) = X(2);
      _mesh.POS(I, 3) = X(3);
    end
    tline = mgetl(fid,1);
    printf('nodes have been read\n')
  elseif part(tline,2) == 'E' & part(tline,3) == 'L' then
    printf('reading elements\n')
    _mesh.nbElm       = mfscanf(1, fid, '%d');
    _mesh.ELE_INFOS   = zeros(_mesh.nbElm, 5);
    _mesh.nbPoints    = 0;
    _mesh.nbLines     = 0;
    _mesh.nbTriangles = 0;
    _mesh.nbQuads     = 0;
    _mesh.nbTets      = 0;
    _mesh.nbHexas     = 0;
    _mesh.nbPrisms    = 0;
    _mesh.nbPyramids  = 0;
    // comment next 8 lines to get "tight" arrays (will slow down reading)
    if init_mem then
      _mesh.POINTS    = zeros(_mesh.nbElm, 2);
      _mesh.LINES     = zeros(_mesh.nbElm, 3);
      _mesh.TRIANGLES = zeros(_mesh.nbElm, 4);
      _mesh.QUADS     = zeros(_mesh.nbElm, 5);
      _mesh.TETS      = zeros(_mesh.nbElm, 5);
      _mesh.HEXAS     = zeros(_mesh.nbElm, 9);
      _mesh.PRISMS    = zeros(_mesh.nbElm, 7);
      _mesh.PYRAMIDS  = zeros(_mesh.nbElm, 6);
    end
    for I = 1:_mesh.nbElm
      _mesh.ELE_INFOS(I, :) = (mfscanf(5, fid, '%d'))';
      NODES_ELEM = mfscanf(_mesh.ELE_INFOS(I, 5), fid, '%d');
      if(_mesh.ELE_INFOS(I, 2) == 15) then // point
	       _mesh.nbPoints = _mesh.nbPoints + 1;
        	_mesh.POINTS(_mesh.nbPoints, 1) = IDS(NODES_ELEM(1));
        	_mesh.POINTS(_mesh.nbPoints, 2) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 1) then // line
        	_mesh.nbLines = _mesh.nbLines + 1;
        	_mesh.LINES(_mesh.nbLines, 1) = IDS(NODES_ELEM(1));
        	_mesh.LINES(_mesh.nbLines, 2) = IDS(NODES_ELEM(2));
        	_mesh.LINES(_mesh.nbLines, 3) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 2) then // triangle
        	_mesh.nbTriangles = _mesh.nbTriangles + 1;
        	_mesh.TRIANGLES(_mesh.nbTriangles, 1) = IDS(NODES_ELEM(1));
        	_mesh.TRIANGLES(_mesh.nbTriangles, 2) = IDS(NODES_ELEM(2));
        	_mesh.TRIANGLES(_mesh.nbTriangles, 3) = IDS(NODES_ELEM(3));
        	_mesh.TRIANGLES(_mesh.nbTriangles, 4) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 3) then // quadrangle
        	_mesh.nbQuads = _mesh.nbQuads + 1;
        	_mesh.QUADS(_mesh.nbQuads, 1) = IDS(NODES_ELEM(1));
        	_mesh.QUADS(_mesh.nbQuads, 2) = IDS(NODES_ELEM(2));
        	_mesh.QUADS(_mesh.nbQuads, 3) = IDS(NODES_ELEM(3));
        	_mesh.QUADS(_mesh.nbQuads, 4) = IDS(NODES_ELEM(4));
        	_mesh.QUADS(_mesh.nbQuads, 5) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 4) then // tetrahedron
        	_mesh.nbTets = _mesh.nbTets + 1;
        	_mesh.TETS(_mesh.nbTets, 1) = IDS(NODES_ELEM(1));
        	_mesh.TETS(_mesh.nbTets, 2) = IDS(NODES_ELEM(2));
        	_mesh.TETS(_mesh.nbTets, 3) = IDS(NODES_ELEM(3));
        	_mesh.TETS(_mesh.nbTets, 4) = IDS(NODES_ELEM(4));
        	_mesh.TETS(_mesh.nbTets, 5) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 5) then // hexahedron
        	_mesh.nbHexas = _mesh.nbHexas + 1;
        	_mesh.HEXAS(_mesh.nbHexas, 1) = IDS(NODES_ELEM(1));
        	_mesh.HEXAS(_mesh.nbHexas, 2) = IDS(NODES_ELEM(2));
        	_mesh.HEXAS(_mesh.nbHexas, 3) = IDS(NODES_ELEM(3));
        	_mesh.HEXAS(_mesh.nbHexas, 4) = IDS(NODES_ELEM(4));
        	_mesh.HEXAS(_mesh.nbHexas, 5) = IDS(NODES_ELEM(5));
        	_mesh.HEXAS(_mesh.nbHexas, 6) = IDS(NODES_ELEM(6));
        	_mesh.HEXAS(_mesh.nbHexas, 7) = IDS(NODES_ELEM(7));
        	_mesh.HEXAS(_mesh.nbHexas, 8) = IDS(NODES_ELEM(8));
        	_mesh.HEXAS(_mesh.nbHexas, 9) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 6) then // prism
        	_mesh.nbPrisms = _mesh.nbPrisms + 1;
        	_mesh.PRISMS(_mesh.nbPrisms, 1) = IDS(NODES_ELEM(1));
        	_mesh.PRISMS(_mesh.nbPrisms, 2) = IDS(NODES_ELEM(2));
        	_mesh.PRISMS(_mesh.nbPrisms, 3) = IDS(NODES_ELEM(3));
        	_mesh.PRISMS(_mesh.nbPrisms, 4) = IDS(NODES_ELEM(4));
        	_mesh.PRISMS(_mesh.nbPrisms, 5) = IDS(NODES_ELEM(5));
        	_mesh.PRISMS(_mesh.nbPrisms, 6) = IDS(NODES_ELEM(6));
        	_mesh.PRISMS(_mesh.nbPrisms, 7) = _mesh.ELE_INFOS(I, 3);
      end
      if(_mesh.ELE_INFOS(I, 2) == 7) then // pyramid
        	_mesh.nbPyramids = _mesh.nbPyramids + 1;
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 1) = IDS(NODES_ELEM(1));
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 2) = IDS(NODES_ELEM(2));
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 3) = IDS(NODES_ELEM(3));
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 4) = IDS(NODES_ELEM(4));
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 5) = IDS(NODES_ELEM(5));
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 6) = IDS(NODES_ELEM(6));
        	_mesh.PYRAMIDS(_mesh.nbPyramids, 7) = _mesh.ELE_INFOS(I, 3);
      end
    end
    tline = mgetl(fid,1);
    printf('elements have been read\n')
  end
end
mclose(fid);
endfunction
