function [mesh, info, ids] = import_medit_mesh(fname)
%IMPORT_INRIA_MESH Summary of this function goes here
%   Detailed explanation goes here

fid = fopen(fname, 'r');

if fid < 0
    error('import_inria_mesh:cannot_open', ...
        'Cannot open mesh file for reading: %s', fname);
end

info = struct('version', [], 'dimension', []);

% Read version
[v, count] = fscanf(fid, ' MeshVersionFormatted %d');
if count ~= 1
    error('import_inria_mesh:format_not_found', ...
        'Could not find format string');
end
info.version = v;

% Read dimension
[d, count] = fscanf(fid, ' Dimension %d');
if count ~= 1
    error('import_inria_mesh:dimension_not_found', ...
        'Could not find dimension string');
end
info.dimension = d;

tokens = {'Vertices', 'Edges', 'Triangles', 'Tetrahedra', 'End'};

n_nodes = 0;
n_edges = 0;
n_tria = 0;
n_tetra = 0;

while ~feof(fid)
    line = strtrim(fgetl(fid));
    if isempty(line) || line(1) == '#';
        continue
    end
    
    i_tok = find(strcmp(line, tokens));
    if isempty(i_tok)
        error('import_inria_mesh:unrecognized_token', ...
            'Unrecognized token found: %s', line);
    end
    
    switch tokens{i_tok}
        case 'Vertices'
            [n_nodes, count] = fscanf(fid, ' %d', 1);
            if count ~= 1
                error('import_inria_mesh:number_of_nodes_not_found', ...
                    'Could not find number of nodes');
            end
            
            % Read the nodes
            nodes = fscanf(fid, '%g', [4 n_nodes]);
        case 'Edges'
            [n_edges, count] = fscanf(fid, ' %d', 1);
            if count ~= 1
               error('import_inria_mesh:number_of_edges_not_found', ...
                    'Could not find number of edges');
            end
            
            edges = fscanf(fid, '%d', [3 n_edges]);
        case 'Triangles'
            [n_tria, count] = fscanf(fid, ' %d', 1);
            if count ~= 1
               error('import_inria_mesh:number_of_triangles_not_found', ...
                    'Could not find number of triangles');
            end
            tria = fscanf(fid, '%d', [4 n_tria]);
        case 'Tetrahedra'
            [n_tetra, count] = fscanf(fid, ' %d', 1);
            if count ~= 1
               error('import_inria_mesh:number_of_tetrahedra_not_found', ...
                    'Could not find number of tetrahedra');
            end
            tetra = fscanf(fid, '%d', [5 n_tetra]);
        case 'End'
            break
    end
end

fclose(fid);

% Assemble the mesh
mesh.Nodes = [(1:n_nodes).', nodes(1:3,:).'];
mesh.Elements = zeros(0, 12);
n_elems = 0;
% Add triangles
if n_tria ~= 0
    mesh.Elements = [mesh.Elements; ...
        [n_elems+(1:n_tria).', 23*ones(n_tria, 1), ones(n_tria, 2), ... 
         tria(1:3, :).', zeros(n_tria, 5)]];
    n_elems = n_elems + n_tria;
end
% Add tetrahedra
if n_tetra ~= 0
    mesh.Elements = [mesh.Elements; ...
        [n_elems+(1:n_tetra).', 34*ones(n_tetra, 1), ones(n_tetra, 2), ... 
         tetra(1:4, :).', zeros(n_tetra, 4)]];
     n_elems = n_elems + n_tetra;
end

[mesh.Materials, mesh.Properties] = default_mat_prop();

% Export the IDs
if nargout > 2
    ids.vertices = nodes(4,:).';
    ids.edges = edges(3,:).';
    ids.elems = zeros(n_elems, 1);
    n_elems = 0;
    if n_tria ~= 0
        ids.elems(n_elems + (1 : n_tria)) = tria(4,:).';
        n_elems = n_elems + n_tria;
    end
    if n_tetra ~= 0
        ids.elems(n_elems + (1 : n_tetra)) = tetra(5,:).';
        n_elems = n_elems + n_tetra;
    end
end

end % of function import_inria_mesh

