%% Import the nodes,edges, and triangles to MATLAB from Gmsh

%{
 points: the coordinates of vertices of each triangle
 edges : the index of points which are located on the edges
 triangles: the connectivity matrix is constructed by the triangles
 index_points: it is required to plot the the solution in Gmsh
 Note  : when the lower boundary is periodic, we create some points on its
       boundary. They are in the list of points but it's possible that the
       nodes on the lower boundary has nothing to do with others triangles.
       That means, they do not exist in the connectivity matrix. Then, we
       have to delete them from points and also the edges.
%}

function [mesh] = Read_Mesh(file)

% Number of points can be found in the fourth line of file.msh
number_points = dlmread(file,'',[4, 0, 4, 0]);


% To read the coordinates of each vertex of triangle, we need to read ...
% the block below the numpoints from line 5 until line number_points + 4.
points = dlmread(file,'',[5, 1, number_points + 4, 2]);


% Number of elements can be found in the (number_points + 7)th line of file.msh
number_elements = dlmread(file,'',[number_points + 7, 0, number_points + 7, 0]);


% To read the lines and triangles, we read the block below the number_element ...
% from line (number_points + 8)th until line (number_points + number_elements + 7)th in the file.msh
elements = dlmread(file,'',[number_points + 8, 1, number_points + number_elements + 7, 7]);

% In the mentioned block, we delete all lines whose the second array is 15...
% because they related to the points which already has been specified.
elements(find(elements(:,1)==15),:) = [];


% In the mentioned block, the edges are found by the line whose the first line...
% is 1 and the index of triangles is found by the line whose the first line is 2.
index_edges = find(elements(:,1)==1);
index_triangles = find(elements(:,1)==2);

% For edges, the useful data is from the 4th to 6th column
% For triangles, the useful data is from the 5th to end.
edges = elements(index_edges,4:6);
triangles = elements(index_triangles,5:end);


points = points.';
edges  = edges.';
triangles = triangles.';

% Determine which point is not a vertex of the triangles (denoted by 0)
Vector = false(1,number_points);
Vector(triangles) = true;

% Find the points which are denoted by zero and one in Vector
index_extpoints = find(~Vector);
index_intpoints = find( Vector);


% Delete the points which are denoted by zero

if(isempty(index_extpoints) == 0) %if the index_extpoints is not empty:
    
    s = zeros(1,number_points);
    
    s(index_intpoints) = 1:length(index_intpoints);
    
    % Delete the additional points on the lower boundary
    points(:,index_extpoints) = [];
    
    % Just for adjusting the indexes
    edges(2:3,:) = s(edges(2:3,:));
    triangles = s(triangles);
end

% === Determinant of the Jacoubian Matrix

Jac =  points(1,triangles(1,:)) .* ( points(2,triangles(2,:)) - points(2,triangles(3,:)) ) + ...
       points(1,triangles(2,:)) .* ( points(2,triangles(3,:)) - points(2,triangles(1,:)) ) + ...
       points(1,triangles(3,:)) .* ( points(2,triangles(1,:)) - points(2,triangles(2,:)) );

  
mesh.numberpoints = length(points);
mesh.points = points;
mesh.edges  = edges;
mesh.triangles = triangles;
mesh.index_intpoints = index_intpoints;
mesh.Area = Jac;


%% Specify the index of the boundary conditions:

[boundary] = Boundary_Points(mesh.edges);

mesh.bnd.Left  = boundary.Left;
mesh.bnd.Right = boundary.Right;
mesh.bnd.Upper = boundary.Upper;
mesh.bnd.Lower = boundary.Lower;

end