read_msh.m

function [p, t] = read_msh(filename)
%function [p, t] = read_msh(filename)
%
% IM 2014-08-14

fid = fopen(filename, 'r');

% read number of nodes
C = textscan(fid, '%d', 1, 'Headerlines', 4);
nNodes = C{1};

% read nodes coordinates
C = textscan(fid, '%d%f%f%f', nNodes);

p = [C{2}(:), C{3}(:)]';

% read elements
C = textscan(fid, '%d', 1, 'CommentStyle', '$');
textscan(fid, '%s', 1, 'Delimiter', '\n');
nElements = C{1};
t = [];

for i = 1:nElements
    
    C = textscan(fid, '%s', 1, 'Delimiter', '\n');
    C = textscan(C{1}{1}, '%f');
    if C{1}(2) == 2
        
        t = [t, C{1}(C{1}(3) + (4:6))];
    end
end

fclose(fid);

% make orientation consistent
X1 = p(1, t(2, :)) - p(1, t(1, :));
Y1 = p(2, t(2, :)) - p(2, t(1, :));
X2 = p(1, t(3, :)) - p(1, t(1, :));
Y2 = p(2, t(3, :)) - p(2, t(1, :));

reverse = X1.*Y2 - X2.* Y1 < 0;

if any(reverse)    
    t([2 3], reverse) = t([3 2], reverse);
end

end