1function gmesh = load_gmsh ( filename )
2   [fid,ierr] = mopen ( filename, 'r');
3	if ierr <> 0 then;
4  		mclose(fid);
5		error('no such file: '+filename)
6	end
7    gmesh = [];
8    gmesh.MIN = zeros(3,1);
9    gmesh.MAX = zeros(3,1);
10
11   while 1
12        endoffile = 0;
13
14        while (1)
15            // jump over irrelevant text
16            tline = mgetl(fid,1);
17            if meof(fid), endoffile=1,break, end
18            if (part(tline,1) == '$'),  break, end
19        end
20        if (endoffile == 1), break, end
21        if (part(tline,1:4) == '$NOD' | part(tline,1:4) == '$Nod')
22            disp('reading nodes')
23            gmesh.nbNod = mfscanf (1, fid, '%d');
24	         gmesh.POS = zeros (gmesh.nbNod, 3);
25            for(I=1:gmesh.nbNod)
26                iNod = mfscanf(1,fid,'%d');
27                X = mfscanf(3,fid,'%g');
28                IDS(iNod) = I;
29                if (I == 1)
30                    gmesh.MIN = X;
31                    gmesh.MAX = X;
32                else
33                    if (gmesh.MAX(1) < X(1)), gmesh.MAX(1)=X(1);end
34                    if (gmesh.MAX(2) < X(2)), gmesh.MAX(2)=X(2);end
35                    if (gmesh.MAX(3) < X(3)), gmesh.MAX(3)=X(3);end
36                    if (gmesh.MIN(1) > X(1)), gmesh.MIN(1)=X(1);end
37                    if (gmesh.MIN(2) > X(2)), gmesh.MIN(2)=X(2);end
38                    if (gmesh.MIN(3) > X(3)), gmesh.MIN(3)=X(3);end
39                    end
40                gmesh.POS(I,1) = X(1);
41                gmesh.POS(I,2) = X(2);
42                gmesh.POS(I,3) = X(3);
43            end
44            // Here we have to read two lines, instead of one, in order to continue.
45            // The reason may be due to an unread end-of-line character from the previous
46            // scan.
47            tline = mgetl(fid,2); // read 2 dummy lines
48            disp('nodes have been read')
49
50        elseif(part(tline,1:4)=='$ELM' | part(tline,1:4)=='$Ele')
51            disp('reading elements')
52            gmesh.nbElm = mfscanf (1,fid, '%d');
53            gmesh.ELE_INFOS = zeros (gmesh.nbElm,5);
54            gmesh.nbLines = 0;
55            gmesh.nbPoints = 0;
56            gmesh.nbTriangles = 0;
57            gmesh.nbQuads = 0;
58            gmesh.nbPrism=0;
59            gmesh.nbPyr=0;
60            gmesh.nbTet = 0;
61            gmesh.nbHex = 0;
62            gmesh.nbQTriangles = 0;
63            gmesh.nbQQuads = 0;
64            gmesh.POINTS=[];
65            gmesh.LINES=[];
66            gmesh.TRIANGLES=[];
67            gmesh.QUADS=[];
68            gmesh.TETS=[];
69            gmesh.HEXA=[];
70            gmesh.PRISM=[];
71            gmesh.PYRAMID=[];
72            gmesh.QTRIANGLES=[];
73            gmesh.QQUADS=[];
74            for (I=1:gmesh.nbElm)
75                if (part(tline,1:4)=='$Ele')
76                   // file is in GMSH version 2.0 format
77                   gmesh.ELE_INFOS(I,1:3) = mfscanf(3,fid,'%d')';
78                   n_of_tags = gmesh.ELE_INFOS(I,3);
79                   gmesh.ELE_INFOS(I,4+(1:n_of_tags)) = mfscanf(n_of_tags,fid,'%d')';
80                   // the number of points
81                   select gmesh.ELE_INFOS(I,2)
82                   case 15 then,  np=1;
83                   case  1 then,  np=2;
84                   case  2 then,  np=3;
85                   case  3 then,  np=4;
86                   case  4 then,  np=4;
87                   case  5 then,  np=8;
88                   case  6 then,  np=6;
89                   case  7 then,  np=5;
90                   case  9 then,  np=6;
91                   case 10 then,  np=9;
92                   end
93                   NODES_ELEM = mfscanf(np,fid,'%d')';
94                else
95                   // version 1.0 format
96                   gmesh.ELE_INFOS(I,:) = mfscanf(5,fid,'%d')';
97                   NODES_ELEM = mfscanf(gmesh.ELE_INFOS(I,5),fid,'%d')';
98                end
99                select gmesh.ELE_INFOS(I,2)
100                case 15 then  // point
101                    gmesh.nbPoints = gmesh.nbPoints + 1;
102                    gmesh.POINTS =[gmesh.POINTS; IDS(NODES_ELEM (1))' I];
103                case 1 then // beam
104                    gmesh.nbLines = gmesh.nbLines + 1;
105                    gmesh.LINES =[gmesh.LINES; IDS(NODES_ELEM (1:2))' I];
106                case 2 then // triangle
107                    gmesh.nbTriangles = gmesh.nbTriangles + 1;
108                    gmesh.TRIANGLES = [gmesh.TRIANGLES; IDS(NODES_ELEM (1:3))' I];
109                case 3 then // quadrangle
110                    gmesh.nbQuads = gmesh.nbQuads + 1;
111                    gmesh.QUADS =[gmesh.QUADS; IDS(NODES_ELEM (1:4))' I];
112                case 4 then // tetrahedron (4 node)
113                    gmesh.nbTet = gmesh.nbTet + 1;
114                    gmesh.TETS = [gmesh.TETS; IDS(NODES_ELEM (1:4))' I];
115                case 5 then // hexahedron (8 nodes)
116                    gmesh.nbHex = gmesh.nbHex + 1;
117                    gmesh.HEXA =[gmesh.HEXA; IDS(NODES_ELEM (1:8))' I];
118                case 6 then // prism (6 nodes)
119                    gmesh.nbPrism = gmesh.nbPrism + 1;
120                    gmesh.PRISM =[gmesh.PRISM; IDS(NODES_ELEM (1:6))' I];
121                case 7 then // pyramid (5 nodes)
122                    gmesh.nbPyr = gmesh.nbPyr + 1;
123                    gmesh.PYRAMID =[gmesh.PYRAMID; IDS(NODES_ELEM (1:5))' I];
124                case 9 then // second order triangle (6 nodes)
125                    gmesh.nbQTriangles = gmesh.nbQTriangles + 1;
126                    gmesh.QTRIANGLES = [gmesh.QTRIANGLES; IDS(NODES_ELEM (1:6))' I];
127                case 10 then // second order quadrangle (9 nodes)
128                    gmesh.nbQQuads = gmesh.nbQQuads + 1;
129                    gmesh.QQUADS = [gmesh.QQUADS; IDS(NODES_ELEM (1:9))' I];
130                else
131                    disp(' ')
132                    warning('Unknown element type ' + string(gmesh.ELE_INFOS(I,2)) + ' !')
133                end
134            end
135            disp('elements have been read')
136            tline = mgetl(fid,1);
137        end  //
138    end  // while
139mclose (fid);
140endfunction
141