1%GMSH2PDETOOLBOX Reads a mesh in msh format, version 1 or 2 and returns the
2%arrays p, e and t according to the MATLAB pde toolbox syntax. Filename
3%refers to the .msh file. Note that only triangular 2D mesh are processed
4%since these are the only elements allowed in pde toolbox.
5%
6%SEE ALSO load_gmsh load_gmsh4 initmesh
7%
8% Copyright (C) 2014 Arthur Levy. https://github.com/arthurlevy/
9%
10% This fucntion uses the routine load_gmsh4 in load_gmsh2.m: Copyright (C)
11% 2007  JP Moitinho de Almeida (moitinho@civil.ist.utl.pt) and  R
12% Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) It is based on
13% load_gmsh.m supplied with gmsh-2.0
14function [p, e, t] = gmsh2pdetoolbox(filename)
15
16%loads the gmsh file using JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
17% and  R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) code.
18mesh_structure = load_gmsh4(filename,-1);
19
20
21%% NODES
22%nodes positions
23disp('importing nodes');
24p = mesh_structure.POS(:,1:2)';
25
26
27%% EDGES
28disp('importing edges');
29%find the edges (ie dimension 1)
30is_edge = (mesh_structure.ELE_INFOS(:,2)==1);
31%add edge data
32e = mesh_structure.ELE_NODES(is_edge,1:2)';
33e(3,:) = 0;
34e(4,:) = 1;
35%tag is important for applying boundary conditions
36e(5,:) = mesh_structure.ELE_TAGS(is_edge,1)';
37e(6,:) = 1;
38e(7,:) = 0;
39
40%% TRIANGLES
41disp('importing triangles')
42is_triangle = (mesh_structure.ELE_INFOS(:,2)==2);
43t = mesh_structure.ELE_NODES(is_triangle,1:3)';
44t(4,:) = 2;
45
46
47
48
49
50function msh = load_gmsh4(filename, which)
51
52%% Reads a mesh in msh format, version 1 or 2
53%
54% Usage:
55% To define all variables m.LINES, M.TRIANGLES, etc
56% (Warning: This creates a very large structure. Do you really need it?)
57%            m = load_gmsh4('a.msh')
58%
59% To define only certain variables (for example TETS and HEXS)
60%            m = load_gmsh4('a.msh', [ 5 6])
61%
62% To define no variables (i.e. if you prefer to use m.ELE_INFOS(i,2))
63%            m = load_gmsh4('a.msh', -1)
64%            m = load_gmsh4('a.msh', [])
65%
66% Copyright (C) 2007  JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
67% and  R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be)
68%
69% based on load_gmsh.m supplied with gmsh-2.0
70%
71% Structure msh always has the following elements:
72%
73% msh.MIN, msh.MAX - Bounding box
74% msh.nbNod - Number of nodes
75% msh.nbElm - Total number of elements
76% msh.nbType(i) - Number of elements of type i (i in 1:19)
77% msh.POS(i,j) - j'th coordinate of node i (j in 1:3, i in 1:msh.nbNod)
78% msh.ELE_INFOS(i,1) - id of element i (i in 1:msh.nbElm)
79% msh.ELE_INFOS(i,2) - type of element i
80% msh.ELE_INFOS(i,3) - number of tags for element i
81% msh.ELE_INFOS(i,4) - Dimension (0D, 1D, 2D or 3D) of element i
82% msh.ELE_TAGS(i,j) - Tags of element i (j in 1:msh.ELE_INFOS(i,3))
83% msh.ELE_NODES(i,j) - Nodes of element i (j in 1:k, with
84%                       k = msh.NODES_PER_TYPE_OF_ELEMENT(msh.ELE_INFOS(i,2)))
85%
86% These elements are created if requested:
87%
88% msh.nbLines - number of 2 node lines
89% msh.LINES(i,1:2) - nodes of line i
90% msh.LINES(i,3) - tag (WHICH ?????) of line i
91%
92% msh.nbTriangles - number of 2 node triangles
93% msh.TRIANGLES(i,1:3) - nodes of triangle i
94% msh.TRIANGLES(i,4) - tag (WHICH ?????) of triangle i
95%
96% Etc
97
98% These definitions need to be updated when new elemens types are added to gmsh
99%
100% msh.Types{i}{1} Number of an element of type i
101% msh.Types{i}{2} Dimension (0D, 1D, 2D or 3D) of element of type i
102% msh.Types{i}{3} Name to add to the structure with elements of type i
103% msh.Types{i}{4} Name to add to the structure with the number of elements of type i
104%
105
106nargchk(1, 2, nargin);
107
108msh.Types = { ...
109    { 2, 1, 'LINES', 'nbLines'}, ... % 1
110    { 3,  2, 'TRIANGLES', 'nbTriangles'}, ...
111    { 4,  2, 'QUADS', 'nbQuads'}, ...
112    { 4,  3, 'TETS', 'nbTets'}, ...
113    { 8,  3, 'HEXAS', 'nbHexas'}, ... %5
114    { 6,  3, 'PRISMS', 'nbPrisms'}, ...
115    { 5,  3, 'PYRAMIDS', 'nbPyramids'}, ...
116    { 3,  1, 'LINES3', 'nbLines3'}, ...
117    { 6,  2, 'TRIANGLES6', 'nbTriangles6'}, ...
118    { 9,  2, 'QUADS9', 'nbQuads9'}, ... % 10
119    { 10,  3, 'TETS10', 'nbTets10'}, ...
120    { 27,  3, 'HEXAS27', 'nbHexas27'}, ...
121    { 18,  3, 'PRISMS18', 'nbPrisms18'}, ...
122    { 14,  3, 'PYRAMIDS14', 'nbPyramids14'}, ...
123    { 1,  0, 'POINTS', 'nbPoints'}, ... % 15
124    { 8,  3, 'QUADS8', 'nbQuads8'}, ...
125    { 20,  3, 'HEXAS20', 'nbHexas20'}, ...
126    { 15,  3, 'PRISMS15', 'nbPrisms15'}, ...
127    { 13,  3, 'PYRAMIDS13', 'nbPyramids13'}, ...
128};
129
130ntypes = length(msh.Types);
131
132if (nargin==1)
133    which = 1:ntypes;
134else
135    if isscalar(which) && which == -1
136        which = [];
137    end
138end
139
140% Could check that "which" is properlly defined....
141
142fid = fopen(filename, 'r');
143fileformat = 0; % Assume wrong file
144
145tline = fgetl(fid);
146if (feof(fid))
147    disp (sprintf('Empty file: %s',  filename));
148    msh = [];
149    return;
150end
151
152%% Read mesh type
153if (strcmp(tline, '$NOD'))
154    fileformat = 1;
155elseif (strcmp(tline, '$MeshFormat'))
156    fileformat = 2;
157    tline = fgetl(fid);
158    if (feof(fid))
159        disp (sprintf('Syntax error (no version) in: %s',  filename));
160        fileformat = 0;
161    end
162    [ form ] = sscanf( tline, '%f %d %d');
163%    if (form(1) ~= 2)
164%        disp (sprintf('Unknown mesh format: %s', tline));
165%        fileformat = 0;
166%    end
167    if (form(2) ~= 0)
168        disp ('Sorry, this program can only read ASCII format');
169        fileformat = 0;
170    end
171    fgetl(fid);    % this should be $EndMeshFormat
172    if (feof(fid))
173        disp (sprintf('Syntax error (no $EndMeshFormat) in: %s',  filename));
174        fileformat = 0;
175    end
176    tline = fgetl(fid);    % this should be $Nodes
177    if (feof(fid))
178        disp (sprintf('Syntax error (no $Nodes) in: %s',  filename));
179        fileformat = 0;
180    end
181end
182
183if (~fileformat)
184    msh = [];
185    return
186end
187
188%% Read nodes
189
190if strcmp(tline, '$NOD') || strcmp(tline, '$Nodes')
191    msh.nbNod = fscanf(fid, '%d', 1);
192    aux = fscanf(fid, '%g', [4 msh.nbNod]);
193    msh.POS = aux(2:4,:)';
194    numids = max(aux(1,:));
195    IDS = zeros(1, numids);
196    IDS(aux(1,:)) = 1:msh.nbNod; % This may not be an identity
197    msh.MAX = max(msh.POS);
198    msh.MIN = min(msh.POS);
199    fgetl(fid); % End previous line
200    fgetl(fid); % Has to be $ENDNOD $EndNodes
201else
202    disp (sprintf('Syntax error (no $Nodes/$NOD) in: %s',  filename));
203    fileformat = 0;
204end
205
206%% Read elements
207tline = fgetl(fid);
208if strcmp(tline,'$ELM') || strcmp(tline, '$Elements')
209    msh.nbElm = fscanf(fid, '%d', 1);
210    % read all info about elements into aux (it is faster!)
211    aux = fscanf(fid, '%g', inf);
212    start = 1;
213    msh.ELE_INFOS = zeros(msh.nbElm, 4); % 1 - id, 2 - type, 3 - ntags, 4 - Dimension
214    msh.ELE_NODES = zeros(msh.nbElm,6); % i - Element, j - ElNodes
215    if (fileformat == 1)
216        ntags = 2;
217    else
218        ntags = 3; % just a prediction
219    end
220    msh.ELE_TAGS = zeros(msh.nbElm, ntags); % format 1: 1 - physical number, 2 - geometrical number
221                                            % format 2: 1 - physical number, 2 - geometrical number, 3 - mesh partition number
222    msh.nbType = zeros(ntypes,1);
223    for I = 1:msh.nbElm
224        if (fileformat == 2)
225            finnish = start + 2;
226            msh.ELE_INFOS(I, 1:3) = aux(start:finnish);
227            ntags = aux(finnish);
228            start = finnish + 1;
229            finnish = start + ntags -1;
230            msh.ELE_TAGS(I, 1:ntags) = aux(start:finnish);
231            start = finnish + 1;
232        else
233            finnish = start + 1;
234            msh.ELE_INFOS(I, 1:2) = aux(start:finnish);
235            start = finnish + 1; % the third element is nnodes, which we get from the type
236            msh.ELE_INFOS(I, 3) = 2;
237            finnish = start + 1;
238            msh.ELE_TAGS(I, 1:2) = aux(start:finnish);
239            start = finnish + 2;
240        end
241        type = msh.ELE_INFOS(I, 2);
242        msh.nbType(type) = msh.nbType(type) + 1;
243        msh.ELE_INFOS(I, 4) = msh.Types{type}{2};
244        nnodes = msh.Types{type}{1};
245        finnish = start + nnodes - 1;
246        msh.ELE_NODES(I, 1:nnodes) = IDS(aux(start:finnish));
247        start=finnish + 1;
248    end
249    fgetl(fid); % Has to be $ENDELM or $EndElements
250else
251    disp (sprintf('Syntax error (no $Elements/$ELM) in: %s',  filename));
252    fileformat = 0;
253end
254
255if (fileformat == 0)
256    msh = [];
257end
258
259fclose(fid);
260
261%% This is used to create the explicit lists for types of elements
262
263for i = which
264    if (~isempty(msh.Types{i}{3}))
265        cmd = sprintf('msh.%s=msh.nbType(%d);', msh.Types{i}{4}, i);
266        eval(cmd);
267        % Dimension
268        cmd = sprintf('msh.%s=zeros(%d,%d);', msh.Types{i}{3}, msh.nbType(i), msh.Types{i}{1}+1);
269        eval(cmd);
270        % Clear nbType for counting, next loop will recompute it
271        msh.nbType(i) = 0;
272    end
273end
274
275for i = 1:msh.nbElm
276    type = msh.ELE_INFOS(i,2);
277    if (find(which == type))
278        if (~isempty(msh.Types{type}{3}))
279            msh.nbType(type) = msh.nbType(type)+1;
280            aux=[ msh.ELE_NODES(i,1:msh.Types{type}{1}), msh.ELE_TAGS(i,1) ];
281            cmd = sprintf('msh.%s(%d,:)=aux;', msh.Types{type}{3}, msh.nbType(type));
282            eval(cmd);
283        end
284    end
285end
286
287return;
288