function [mesh, info, group] = import_gmsh_mesh(fname)

fid = fopen(fname);
ftext = textscan(fid, '%s', 'Delimiter', '\n');
fclose(fid);

% Commands start with $ sign
rxp_command = '^\$';

r = regexp(ftext{1}, rxp_command, 'start');

k = 1;
commands = cell(0,2);


for i = 1 : numel(r)
    if ~isempty(r{i})
        commands{k, 1} = ftext{1}{i};
        commands{k, 2} = i;
        k = k+1;
    end
end

% Obtain info
info = read_info(ftext{1}, commands);

if ~strcmpi(info.file_type, 'ASCII')
    error('import_gmsh_mesh:unsupported_format', ...
        'Only ASCII format is supported, format found: %s', info.file_type);
end

if ~ismember(round(info.version), [2 4])
    error('import_gmsh_mesh:unsuppoerted_version', ...
        'Only versions 2.x and 4.x are supported, version found: %g', info.version);
end

% Read physical
group.physical = read_physical(ftext{1}, commands);

% Read entities
group.entities = read_entities(ftext{1}, commands);

% Read nodes
[mesh.Nodes, group.node2blk, group.nblk2ent] = read_nodes(ftext{1}, info, commands);

% Read elements
[mesh.Elements, group.elem2blk, group.eblk2ent] = read_elements(ftext{1}, info, commands);

% Create default materials and properties
[mesh.Materials, mesh.Properties] = default_mat_prop();


end % of function import_gmsh_mesh

function [start_ind, stop_ind] = find_start_stop(start_exp, stop_exp, commands)
    if ~isempty(start_exp)
        start_ind = regexp(commands(:,1), start_exp, 'start', 'once');
        start_ind = find(~cellfun(@isempty, start_ind), 1, 'first');
    else
        start_ind = [];
    end
    
    if ~isempty(stop_exp)
        stop_ind = regexp(commands(:,1), stop_exp, 'start', 'once');
        stop_ind = find(~cellfun(@isempty, stop_ind), 1, 'first');
    else
        stop_ind = [];
    end
end

function info = read_info(ftext, commands)
    
    start_exp = '^\$MeshFormat';    
    start_ind = find_start_stop(start_exp, [], commands);
    
    % If no info found
    if isempty(start_ind)
        info = [];
        return
    end
    
    % Starting line of information
    start_line = ftext{commands{start_ind,2}+1};
    
    % Process the structure
    data = sscanf(start_line, '%f', 3);
    
    info.version = data(1);
    if data(2) == 0
        info.file_type = 'ASCII';
    else
        info.file_type = data(2);
    end
    info.data_size = data(3);
end

function phys = read_physical(ftext, commands)
    start_exp = '^\$PhysicalNames';
    stop_exp = '^\$EndPhysicalNames';
    
    [start_ind, stop_ind] = find_start_stop(start_exp, stop_exp, commands);
    
    if isempty(start_ind) || isempty(stop_ind)
        phys = cell(0, 3);
        return
    end
    
    % Read number of physical entities
    start_line_ind = commands{start_ind, 2} + 1;
    n_phys = sscanf(ftext{start_line_ind}, '%d', 1);
    
    % Check if reading was successful
    if numel(n_phys) ~= 1 || n_phys < 0 
        error('import_gmsh_mesh:physical_num_error', ...
            'Error reading number of physical names in line %d', start_line_ind);
    end 
    
    % Preallocate
    phys = cell(n_phys, 3);
    
    % Go through physical names
    for i_phys = 1 : n_phys
        % Process line
        data = sscanf(ftext{start_line_ind + i_phys}, '%d %d \"%[^\"]\"', [1, inf]);
        
        % Check if success
        if numel(data) < 2
            error('import_gmsh_mesh:physical_format_error', ...
                'Error reading physical name #%d in line %d', i_phys, start_line_ind + i_phys);
        end
        
        % Store data
        phys{i_phys, 1} = data(1);
        phys{i_phys, 2} = data(2);
        if numel(data) >= 3
            phys{i_phys, 3} = char(data(3:end));
        else
            phys(i_phys, 3) = '';
        end
    end
end % of function read_physical

function entities = read_entities(ftext, commands)
    % Entity types
    types = {'points', 3; 'curves', 6; 'surfaces', 6; 'volumes', 6};
    n_types = size(types, 1);

    start_exp = '^\$Entities';
    stop_exp = '^\$EndEntities';
    
    [start_ind, stop_ind] = find_start_stop(start_exp, stop_exp, commands);

    args = [reshape(types(:,1), 1, []); reshape(cell(n_types, 1), 1, [])];
    entities = struct(args{:});
    
    % If no entities found
    if isempty(start_ind) || isempty(stop_ind)
        return
    end
   
    types = {'points', 3; 'curves', 6; 'surfaces', 6; 'volumes', 6};
    n_types = size(types, 1);
    start_line_ind = commands{start_ind, 2} + 1;
    
    % Read number of entities for each type
    n_per_type = sscanf(ftext{start_line_ind}, '%d', [1 n_types]);
    
    if numel(n_per_type) ~= n_types || any(n_per_type < 0)
        error('import_gmsh_mesh:entity_num_per_type', ...
            'Error reading number of entities in line %d', start_line_ind);
    end
    
    % Cycle through types
    i_line = 1;
    for i_type = 1 : n_types
        type_name = types{i_type, 1};
        type_num_data = types{i_type, 2};
        entities.(type_name) = zeros(n_per_type(i_type), type_num_data + 3);
        for i_ent = 1 : n_per_type(i_type)
            % data: [tag, type_num_data x double, num_phys, phys, num_node, nodes]
            data = sscanf(ftext{start_line_ind + i_line}, '%g', [1 inf]);
            
            if numel(data) < type_num_data + 2
                error('import_gms_mesh:entity_format', ...
                    'Error reading ''%s'' entity #%d in line %d', types{i_type,1}, i_ent, start_line_ind+i_line);
            end
            % Keep what is needed from data
            n = type_num_data + 2 + data(type_num_data + 2);
            % Write into matrix
            entities.(type_name)(i_ent, 1 : n) = data(1 : n);
            
            % Next line
            i_line = i_line + 1;
        end
    end
end

function [nodes, node2blk, nblk2ent] = read_nodes(ftext, info, commands)
    start_exp = '^\$Nodes';
    stop_exp = '^\$EndNodes';
    
    [start_ind, stop_ind] = find_start_stop(start_exp, stop_exp, commands);
    
    % If no nodes found
    if isempty(start_ind) || isempty(stop_ind)
        nodes = [];
        node2blk = [];
        nblk2ent = [];
        return
    end
    
    switch round(info.version)
        case 2
            start_line_ind = commands{start_ind,2}+2;
            stop_line_ind = commands{stop_ind,2}-1;
            
            nNodes = stop_line_ind - start_line_ind + 1;
            nodes = zeros(nNodes, 4);

            % ElemID, x, y, z format
            for iNode = 1 : nNodes
                nodes(iNode, :) = sscanf(ftext{start_line_ind + iNode -1}, '%f', 4);
            end
        case 4
            % Read number of entities
            % Format of line : numEntityBlocks(size_t) numNodes(size_t)
            %                  minNodeTag(size_t) maxNodeTag(size_t)
            start_line_ind = commands{start_ind,2} + 1;
            line_cnt = 0;
            node_cnt = 0;
            data = sscanf(ftext{start_line_ind + line_cnt}, '%d %d %d %d');
            
            if numel(data) ~= 4
                error('import_gmsh_mesh:node_block_info_error', ...
                    'Could not read node block information from line %d', start_line_ind + line_cnt);
            end
            line_cnt = line_cnt + 1;
            
            num_blocks = data(1);
            num_nodes = data(2);
            
            % Preallocate
            nodes = zeros(num_nodes, 4);
            node2blk = zeros(num_nodes, 1);
            nblk2ent = zeros(num_blocks, 2);
            
            % Read each block
            for i_block = 1 : num_blocks
                % Read block information
                % Format of line:
                %   entityDim(int) entityTag(int) parametric(int; 0 or 1)
                %   numNodesInBlock(size_t)
                data = sscanf(ftext{start_line_ind + line_cnt}, '%d %d %d %d');
                
                if numel(data) ~= 4
                    error('import_gmsh_mesh:node_entity_info_error', ...
                        'Could not read node entity information from line %d', start_line_ind + line_cnt);
                end
                line_cnt = line_cnt + 1;
                
                % Store the block data
                nblk2ent(i_block, 1:2) = data(1:2);
                n_nodes_in_block = data(4);
                
                % Read the node IDs
                node_ids = str2double(ftext(start_line_ind + line_cnt + (0 : n_nodes_in_block-1)));
                if any(isnan(node_ids))
                    err_idx = find(isnan(node_ids), 1, 'first');
                    error('import_gmsh_mesh:invalid_node_ID', ...
                        'Invalid node ID found for node #%d in block %d, line %d', ...
                            err_idx, i_block, start_line_ind + line_cnt + err_idx -1);
                end
                line_cnt = line_cnt + n_nodes_in_block;
                nodes(node_cnt + (1 : n_nodes_in_block), 1) = node_ids;
                node2blk(node_cnt + (1 : n_nodes_in_block)) = i_block;
                
                % Read the node coordinates
                coords = cellfun(@(x)sscanf(x, '%g %g %g'), ...
                    ftext(start_line_ind+line_cnt+(0:n_nodes_in_block-1)), 'UniformOutput', false);
                if ~all(cellfun(@isnumeric, coords)) || any(cellfun(@numel, coords) ~= 3)
                    err_idx = find(~cellfun(@isnumeric, coords) | ~(cellfun(@numel, coords) ~= 3), 1, 'first');
                    error('import_gmsh_mesh:invalid_node_coords', ...
                        'Invalid node coordinates found for node #%d in block %d, line %d', ...
                            err_idx, i_block, start_line_ind + line_cnt + err_idx-1);
                end
                line_cnt = line_cnt + n_nodes_in_block;
                coords = reshape(cell2mat(coords), 3, []).';
                nodes(node_cnt + (1 : n_nodes_in_block), 2:4) = coords;
                
                node_cnt = node_cnt + n_nodes_in_block;
            end
    end
    
end % of function read_nodes

function [elements, elem2blk, eblk2ent] = read_elements(ftext, info, commands)
    start_exp = '^\$Elements';
    stop_exp = '^\$EndElements';
    
    start_ind = regexp(commands(:,1), start_exp, 'start', 'once');
    start_ind = find(~cellfun(@isempty, start_ind), 1, 'first');
    
    stop_ind = regexp(commands(:,1), stop_exp, 'start', 'once');
    stop_ind = find(~cellfun(@isempty, stop_ind), 1, 'first');
    
    % If no info found
    if isempty(start_ind) || isempty(stop_ind)
        elements = [];
        elem2blk = [];
        eblk2ent = [];
        return
    end
    
    switch round(info.version)
        case 2
            % Version 2, ASCII
            start_line_ind = commands{start_ind,2}+2;
            stop_line_ind = commands{stop_ind,2}-1;
    
            nElem = stop_line_ind - start_line_ind + 1;
            elements = zeros(nElem, 12);
            jElem = 0;
            for iElem = 1 : nElem
                data = sscanf(ftext{start_line_ind + iElem - 1}, '%f', Inf);
                if (numel(data) < 2)
                    continue
                end
                % Select elem type
                switch data(2)
                    % 3-node tria
                    case 2
                        elements(jElem+1, 1:4) = [data(1) 23 1 1];
                        elements(jElem+1, 4+(1:3)) = data(end-2 : end);
                        jElem = jElem+1;
                    % 4-node quadrangle
                    % TODO: check node ordering 
                    case 3
                        elements(jElem+1, 1:4) = [data(1) 24 1 1];
                        elements(jElem+1, 4+(1:4)) = data(end-3 : end);
                        jElem = jElem+1;
                    % 4-node tetrahedron
                    case 4
                        elements(jElem+1, 1:4) = [data(1) 34 1 1];
                        elements(jElem+1, 4+(1:4)) = data(end-3 : end);
                        jElem = jElem+1;
                    % 8-node hexahedron 
                    % TODO: check node ordering
                    case 5
                        elements(jElem+1, 1:4) = [data(1) 38 1 1];
                        elements(jElem+1, 4+(1:8)) = data(end-7 : end);
                        jElem = jElem+1;
                    % 6-node prism
                    % TODO: check node ordering
                    case 6
                        elements(jElem+1, 1:4) = [data(1) 36 1 1];
                        elements(jElem+1, 4+(1:6)) = data(end-5 : end);
                        jElem = jElem+1;
                end
            end
            % Drop unread elements
            elements = elements(1 : jElem, :);
        
        case 4
            % Version 4, ASCII
            % Read number of entities
            % Format of line : numEntityBlocks(size_t) numElements(size_t)
            %                  minElementTag(size_t) maxElementTag(size_t)
            start_line_ind = commands{start_ind,2} + 1;
            line_cnt = 0;
            elem_cnt = 0;
            data = sscanf(ftext{start_line_ind + line_cnt}, '%d %d %d %d');
            
            if numel(data) ~= 4
                error('import_gmsh_mesh:elem_block_info_error', ...
                    'Could not read element block information from line %d', start_line_ind + line_cnt);
            end
            line_cnt = line_cnt + 1;
            
            n_blocks = data(1);
            n_elems = data(2);
            
            gmsh_to_nihu = [
                1 12 % 2-node line 
                2 23 % 3-node triangle 
                3 24 % 4-node quadrangle
                4 34 % 4-node tetrahedron
                5 38 % 8-node hexahedron
                6 36 % 6-node prism. 
                7 35 % -node pyramid
            ];
            
            % Pre-allocate
            elements = zeros(n_elems, 12);
            elem2blk = zeros(n_elems, 1);
            eblk2ent = zeros(n_blocks, 2);
        
            for i_block = 1 : n_blocks
                % Read block information
                % Format of line:
                %   entityDim(int) entityTag(int) elementType(int; see below)
                %   numElementsInBlock(size_t)
                data = sscanf(ftext{start_line_ind + line_cnt}, '%d %d %d %d');
                
                if numel(data) ~= 4
                    error('import_gmsh_mesh:element_entity_info_error', ...
                        'Could not read element entity information from line %d', start_line_ind + line_cnt);
                end
                line_cnt = line_cnt + 1;
                
                eblk2ent(i_block, 1:2) = data(1:2);
                elem_type = data(3);
                n_elems_in_block = data(4);
                
                [is, type_idx] = ismember(elem_type, gmsh_to_nihu(:,1));
                if is
                    % Number of nodes per element
                    n_node_per_elem = mod(gmsh_to_nihu(type_idx, 2), 10);
                    % Format: tag, node IDs
                    format = ['%d', repmat(' %d', 1, n_node_per_elem)];
                    
                    % Read indices
                    indices = cellfun(@(x)sscanf(x, format), ...
                        ftext(start_line_ind+line_cnt+(0:n_elems_in_block-1)), 'UniformOutput', false);
                    
                    % Check for invalid format
                    if ~all(cellfun(@isnumeric, indices)) || any(cellfun(@numel, indices) ~= n_node_per_elem+1)
                        err_idx = find(~cellfun(@isnumeric, indices) | ~(cellfun(@numel, indices) ~= n_node_per_elem+1), 1, 'first');
                        error('import_gmsh_mesh:invalid_elem_indices', ...
                            'Invalid element indices for elem #%d in block %d, line %d', ...
                                err_idx, i_block, start_line_ind + line_cnt + err_idx-1);
                    end
                    
                    indices = reshape(cell2mat(indices), n_node_per_elem+1, []).';
                    
                    % Store the elements
                    elements(elem_cnt + (1 : n_elems_in_block), [1 4+(1:n_node_per_elem)]) = indices;
                    elements(elem_cnt + (1 : n_elems_in_block), 2) = gmsh_to_nihu(type_idx, 2);
                    elements(elem_cnt + (1 : n_elems_in_block), [3 4]) = 1;
                    
                    elem2blk(elem_cnt + (1 : n_elems_in_block)) = i_block;
                    
                    elem_cnt = elem_cnt + n_elems_in_block;
                end
                    
                line_cnt = line_cnt + n_elems_in_block;
            end % of cycle over blocks
            % Only keep elements that were actually read
            elements = elements(1 : elem_cnt, :);
            elem2blk = elem2blk(1 : elem_cnt);
    end
end % of function read_elements
