gmsh_read_mesh

[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(filename)
 Function to read in a mesh model from Gmsh and saves it in
 five arrays; surface (srf), veritices (vtx), face no. (fc)
 volume (simp) and edges (edg)

 srf        = The surfaces indices into vtx
 simp       = The volume indices into vtx
 vtx        = The vertices matrix
 fc         = A one column matrix containing the face numbers
 edg        = Edge segment information
 filename   = Name of file containing NetGen .vol information
 mat_ind    = Material index
 phys_names = Structure of "Physical" entities in the mesh
              .dim   = dimension
              .name  = name (string)
              .tag   = physical tag
              .nodes = N-x-dim array of indices into vtx
 
 This mostly works on GMSH v2. I very basic (and probably
   buggy GMSH v4 reader is now included)

0001 function[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(filename)
0002 %[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(filename)
0003 % Function to read in a mesh model from Gmsh and saves it in
0004 % five arrays; surface (srf), veritices (vtx), face no. (fc)
0005 % volume (simp) and edges (edg)
0006 %
0007 % srf        = The surfaces indices into vtx
0008 % simp       = The volume indices into vtx
0009 % vtx        = The vertices matrix
0010 % fc         = A one column matrix containing the face numbers
0011 % edg        = Edge segment information
0012 % filename   = Name of file containing NetGen .vol information
0013 % mat_ind    = Material index
0014 % phys_names = Structure of "Physical" entities in the mesh
0015 %              .dim   = dimension
0016 %              .name  = name (string)
0017 %              .tag   = physical tag
0018 %              .nodes = N-x-dim array of indices into vtx
0019 %
0020 % This mostly works on GMSH v2. I very basic (and probably
0021 %   buggy GMSH v4 reader is now included)
0022 
0023 % $Id: gmsh_read_mesh.m 6043 2019-12-31 17:27:26Z aadler $
0024 % (C) 2009 Bartosz Sawicki. Licensed under GPL V2
0025 % Modified by James Snyder <jbsnyder@fanplastic.org>
0026 
0027 if ischar(filename) && strcmp(filename,'UNIT_TEST'); do_unit_test; return; end
0028 
0029 fid = fopen(filename,'r');
0030 phys_names = [];
0031 while 1
0032     tline = fgetl(fid);
0033     if ~ischar(tline); fclose(fid); break; end
0034 
0035     if strcmp(tline,'$MeshFormat')
0036        gmshformat = parse_format( fid);
0037     elseif strcmp(tline,'$Elements')
0038        elements= parse_elements( fid, gmshformat );
0039     elseif strcmp(tline,'$Nodes')
0040        nodes= get_lines_with_nodes( fid, gmshformat );
0041     elseif strcmp(tline,'$PhysicalNames')
0042        phys_names= parse_names( fid, gmshformat );
0043     end
0044 end
0045 
0046 if ~isempty(phys_names)
0047     for i = 1:length(phys_names)
0048         tmpelements = find(arrayfun(@(x)x.phys_tag==phys_names(i).tag,elements));
0049         phys_names(i).nodes = cat(1,elements(tmpelements).simp);
0050     end
0051 end
0052 
0053 edg = [];
0054 bc = [];
0055 mat_ind = [];
0056 
0057 % Select 2d vs 3d model by checking if Z is all the same
0058 if length( unique( nodes(:,4) ) ) > 1 
0059     vtx = nodes(:,2:4);
0060     % Type 2: 3-node triangle
0061 keyboard
0062     tri = find(arrayfun(@(x)x.type==2,elements));
0063     % Type 4: 4-node tetrahedron
0064     tet = find(arrayfun(@(x)x.type==4,elements));
0065     simp = cat(1,elements(tet).simp);
0066     srf = cat(1,elements(tri).simp);
0067 else
0068     vtx = nodes(:,2:3);
0069     tri = find(arrayfun(@(x)x.type==2,elements));
0070     simp = cat(1,elements(tri).simp);
0071     srf = [];
0072 end
0073 
0074 elemtags = cat(1,elements.phys_tag);
0075 fc = elemtags(tri,1);
0076 end
0077 
0078 
0079 function mat = get_lines_with_nodes( fid, gmshformat )
0080    tline = fgetl(fid);
0081    n_rows = parse_rows(tline,gmshformat);
0082    switch floor(gmshformat)
0083 % Line Format:
0084 % node-number x-coord y-coord z-coord
0085      case 2; mat= fscanf(fid,'%f',[4,n_rows])';
0086      case 4; mat= zeros(n_rows,4);
0087         while (true)
0088           tline = fgetl(fid);
0089           n = sscanf(tline, '%d')';
0090           n_block = n(4);
0091           for b = 1:n_block;
0092              tline = fgetl(fid);
0093              el= sscanf(tline, '%f')';
0094              mat(el(1),:) = el(1:end);
0095           end
0096           if (el(1) == n_rows); break; end % got them all
0097       end
0098       tline = fgetl(fid); % get the EndElements
0099 
0100      otherwise; error('cant parse gmsh file of this format');
0101    end
0102 end
0103 
0104 function gmshformat = parse_format( fid);
0105    tline = fgetl(fid);
0106    rawformat = sscanf(tline,'%f');
0107    tline = fgetl(fid); % should be EndMeshFormat
0108    gmshformat = rawformat(1);
0109 end
0110 
0111 function n_rows = parse_rows(tline, gmshformat);
0112    n_rows = sscanf(tline,'%d');
0113    switch floor(gmshformat)
0114      case 2; n_rows = n_rows(1);
0115      case 4; n_rows = n_rows(2);
0116      otherwise; error('cant parse gmsh file of this format');
0117    end
0118 end
0119 
0120 function names = parse_names( fid, gmshformat )
0121 % Line Format:
0122 % physical-dimension physical-number "physical-name"
0123 tline = fgetl(fid);
0124 n_rows = parse_rows(tline,gmshformat);
0125 names = struct('tag',{},'dim',{},'name',{});
0126 for i = 1:n_rows
0127     tline = fgetl(fid);
0128     if exist('OCTAVE_VERSION')
0129         parts = strsplit(tline,' ');
0130     else
0131         parts = regexp(tline,' ','split');
0132     end
0133     nsz = size(names,2)+1;
0134     names(nsz).dim = str2double( parts(1) );
0135     names(nsz).tag = str2double( parts(2) );
0136     tname = parts(3);
0137     names(nsz).name = strrep(tname{1},'"','');
0138 end
0139 end
0140 
0141 function elements = parse_elements( fid, gmshformat )
0142    tline = fgetl(fid);
0143    n_rows = parse_rows(tline,gmshformat);
0144    switch floor(gmshformat)
0145      case 2; elements = parse_v2_elements(fid,n_rows);
0146      case 4; elements = parse_v4_elements(fid,n_rows);
0147      otherwise error('cant parse this file type');
0148    end
0149 end
0150 
0151 function elements = parse_v4_elements(fid,n_rows);
0152 % http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
0153 % Partial implementation ... look here
0154    elements = struct('simp',{},'type',{},'phys_tag',{});
0155    while (true)
0156        tline = fgetl(fid);
0157        n = sscanf(tline, '%d')';
0158        n_block = n(4);
0159        n_type = n(2);
0160        for b = 1:n_block;
0161           tline = fgetl(fid);
0162           el= sscanf(tline, '%d')';
0163           if n_type>1; %not interested in points
0164              elements(end+1) = struct( ...
0165               'simp',el(2:end),'type',n_type, ...
0166               'phys_tag',0);
0167           end
0168        end
0169        if (el(1) == n_rows); break; end % got them all
0170    end
0171    tline = fgetl(fid); % get the EndElements
0172 end
0173 
0174 function elements = parse_v2_elements(fid,n_rows);
0175 % Line Format:
0176 % elm-number elm-type number-of-tags < tag > ... node-number-list
0177 elements(n_rows).simp = [];
0178 elements(n_rows).phys_tag = [];
0179 elements(n_rows).geom_tag = [];
0180 elements(n_rows).type = [];
0181 
0182 for i = 1:n_rows
0183     tline = fgetl(fid);
0184     n = sscanf(tline, '%d')';
0185 %     nsz = size(elements,2)+1;
0186     elements(i).simp = n(n(3) + 4:end);
0187     %
0188     elements(i).type = n(2);
0189     if n(3) > 0 % get tags if they exist
0190         % By default, first tag is number of parent physical entity
0191         % second is parent elementary geometrical entity
0192         % third is number of parent mesh partitions followed by
0193         % partition ids
0194         tags = n(4:3+n(3));
0195         if length(tags) >= 1
0196             elements(i).phys_tag = tags(1);
0197             if length(tags) >= 2
0198                 elements(i).geom_tag = tags(2);
0199             end
0200         end
0201     end
0202 end
0203 end
0204 
0205 function do_unit_test
0206    tmpnam = tempname;
0207 
0208    fid = fopen(tmpnam,'w');
0209    fprintf(fid,gmshv4file);
0210    fclose(fid);
0211    [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(tmpnam);
0212     unit_test_cmp('v4 vtx ',vtx(2:3,:),[1,0;-1,0])
0213     unit_test_cmp('v4 simp',simp(2:3,:),[3,7,12; 12, 7,14]);
0214 
0215    fid = fopen(tmpnam,'w');
0216    fprintf(fid,gmshv2file);
0217    fclose(fid);
0218    [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(tmpnam);
0219     unit_test_cmp('v2 vtx ',vtx(2:3,:),[1,0;-1,0])
0220     unit_test_cmp('v2 simp',simp(2:3,:),[2,4,15; 14,17,19]);
0221     
0222 end
0223 
0224 % Example of gmsh v4 file
0225 function t = gmshv4file; t=[ ...
0226 '$MeshFormat\n' ...
0227 '4 0 8\n' ...
0228 '$EndMeshFormat\n' ...
0229 '$Entities\n' ...
0230 '3 2 1 0\n' ...
0231 '1 0 0 0 0 0 0 0 \n' ...
0232 '2 1 0 0 1 0 0 0 \n' ...
0233 '3 -1 0 0 -1 0 0 0 \n' ...
0234 '1 -1 0 0 1 0.984807753012208 0 0 2 2 -3 \n' ...
0235 '2 -1 -0.984807753012208 0 1 0 0 0 2 3 -2 \n' ...
0236 '4 -1 -0.984807753012208 0 1 0.984807753012208 0 0 2 1 2 \n' ...
0237 '$EndEntities\n' ...
0238 '$Nodes\n' ...
0239 '6 14\n' ...
0240 '1 0 0 1\n' ...
0241 '1 0 0 0\n' ...
0242 '2 0 0 1\n' ...
0243 '2 1 0 0\n' ...
0244 '3 0 0 1\n' ...
0245 '3 -1 0 0\n' ...
0246 '1 1 0 3\n' ...
0247 '4 0.7071067795891069 0.7071067827839882 0\n' ...
0248 '5 -4.62506049552981e-09 1 0\n' ...
0249 '6 -0.7071067828954029 0.7071067794776923 0\n' ...
0250 '2 1 0 3\n' ...
0251 '7 -0.7071067795891069 -0.7071067827839882 0\n' ...
0252 '8 4.62506049552981e-09 -1 0\n' ...
0253 '9 0.7071067828954029 -0.7071067794776923 0\n' ...
0254 '4 2 0 5\n' ...
0255 '10 -7.225004064876332e-18 5.806922560102616e-17 0\n' ...
0256 '11 0.4208610471654206 -0.174326352879787 0\n' ...
0257 '12 -0.4208610471654206 0.1743263528797871 0\n' ...
0258 '13 0.1601886191568594 0.3867295406713929 0\n' ...
0259 '14 -0.1601886191568594 -0.3867295406713928 0\n' ...
0260 '$EndNodes\n' ...
0261 '$Elements\n' ...
0262 '6 27\n' ...
0263 '1 0 15 1\n' ...
0264 '1 1 \n' ...
0265 '2 0 15 1\n' ...
0266 '2 2 \n' ...
0267 '3 0 15 1\n' ...
0268 '3 3 \n' ...
0269 '1 1 1 4\n' ...
0270 '4 2 4 \n' ...
0271 '5 4 5 \n' ...
0272 '6 5 6 \n' ...
0273 '7 6 3 \n' ...
0274 '2 1 1 4\n' ...
0275 '8 3 7 \n' ...
0276 '9 7 8 \n' ...
0277 '10 8 9 \n' ...
0278 '11 9 2 \n' ...
0279 '4 2 2 16\n' ...
0280 '12 2 4 11 \n' ...
0281 '13 3 7 12 \n' ...
0282 '14 12 7 14 \n' ...
0283 '15 11 4 13 \n' ...
0284 '16 5 6 12 \n' ...
0285 '17 8 9 11 \n' ...
0286 '18 8 11 14 \n' ...
0287 '19 5 12 13 \n' ...
0288 '20 10 11 13 \n' ...
0289 '21 10 12 14 \n' ...
0290 '22 4 5 13 \n' ...
0291 '23 7 8 14 \n' ...
0292 '24 12 10 13 \n' ...
0293 '25 11 10 14 \n' ...
0294 '26 6 3 12 \n' ...
0295 '27 9 2 11 \n' ...
0296 '$EndElements\n'];
0297 end
0298 
0299 % Example of gmsh v2 file
0300 function t = gmshv2file; t=[ ...
0301 '$MeshFormat\n' ...
0302 '2.2 0 8\n' ...
0303 '$EndMeshFormat\n' ...
0304 '$Nodes\n' ...
0305 '26\n' ...
0306 '1 0 0 0\n' ...
0307 '2 1 0 0\n' ...
0308 '3 -1 0 0\n' ...
0309 '4 0.8660254037855228 0.4999999999981222 0\n' ...
0310 '5 0.500000000002617 0.8660254037829277 0\n' ...
0311 '6 3.95800631512864e-12 1 0\n' ...
0312 '7 -0.4999999999980611 0.866025403785558 0\n' ...
0313 '8 -0.8660254037843073 0.5000000000002276 0\n' ...
0314 '9 -0.8660254037855228 -0.4999999999981222 0\n' ...
0315 '10 -0.500000000002617 -0.8660254037829277 0\n' ...
0316 '11 -3.95800631512864e-12 -1 0\n' ...
0317 '12 0.4999999999980611 -0.866025403785558 0\n' ...
0318 '13 0.8660254037843073 -0.5000000000002276 0\n' ...
0319 '14 5.368222989669768e-17 1.262108679519644e-17 0\n' ...
0320 '15 0.5000000000000001 0.1339745962149803 0\n' ...
0321 '16 -0.5 -0.1339745962149803 0\n' ...
0322 '17 0.133974596217502 0.5000000000002625 0\n' ...
0323 '18 -0.133974596217502 -0.5000000000002626 0\n' ...
0324 '19 -0.3660254037841135 0.3660254037850291 0\n' ...
0325 '20 0.3660254037841134 -0.3660254037850291 0\n' ...
0326 '21 0.6830127018922291 -0.1830127018925211 0\n' ...
0327 '22 -0.683012701892229 0.1830127018925211 0\n' ...
0328 '23 0.5000000000014104 0.4999999999990732 0\n' ...
0329 '24 -0.5000000000014104 -0.4999999999990732 0\n' ...
0330 '25 -0.1830127018901805 0.6830127018929458 0\n' ...
0331 '26 0.1830127018901805 -0.6830127018929458 0\n' ...
0332 '$EndNodes\n' ...
0333 '$Elements\n' ...
0334 '51\n' ...
0335 '1 15 2 0 1 1\n' ...
0336 '2 15 2 0 2 2\n' ...
0337 '3 15 2 0 3 3\n' ...
0338 '4 1 2 0 1 2 4\n' ...
0339 '5 1 2 0 1 4 5\n' ...
0340 '6 1 2 0 1 5 6\n' ...
0341 '7 1 2 0 1 6 7\n' ...
0342 '8 1 2 0 1 7 8\n' ...
0343 '9 1 2 0 1 8 3\n' ...
0344 '10 1 2 0 2 3 9\n' ...
0345 '11 1 2 0 2 9 10\n' ...
0346 '12 1 2 0 2 10 11\n' ...
0347 '13 1 2 0 2 11 12\n' ...
0348 '14 1 2 0 2 12 13\n' ...
0349 '15 1 2 0 2 13 2\n' ...
0350 '16 2 2 0 4 3 9 16\n' ...
0351 '17 2 2 0 4 2 4 15\n' ...
0352 '18 2 2 0 4 14 17 19\n' ...
0353 '19 2 2 0 4 14 18 20\n' ...
0354 '20 2 2 0 4 5 6 17\n' ...
0355 '21 2 2 0 4 10 11 18\n' ...
0356 '22 2 2 0 4 7 8 19\n' ...
0357 '23 2 2 0 4 12 13 20\n' ...
0358 '24 2 2 0 4 14 16 18\n' ...
0359 '25 2 2 0 4 14 15 17\n' ...
0360 '26 2 2 0 4 14 19 16\n' ...
0361 '27 2 2 0 4 14 20 15\n' ...
0362 '28 2 2 0 4 6 7 25\n' ...
0363 '29 2 2 0 4 11 12 26\n' ...
0364 '30 2 2 0 4 2 21 13\n' ...
0365 '31 2 2 0 4 3 22 8\n' ...
0366 '32 2 2 0 4 9 24 16\n' ...
0367 '33 2 2 0 4 4 23 15\n' ...
0368 '34 2 2 0 4 17 25 19\n' ...
0369 '35 2 2 0 4 18 26 20\n' ...
0370 '36 2 2 0 4 3 16 22\n' ...
0371 '37 2 2 0 4 2 15 21\n' ...
0372 '38 2 2 0 4 9 10 24\n' ...
0373 '39 2 2 0 4 4 5 23\n' ...
0374 '40 2 2 0 4 15 20 21\n' ...
0375 '41 2 2 0 4 16 19 22\n' ...
0376 '42 2 2 0 4 5 17 23\n' ...
0377 '43 2 2 0 4 10 18 24\n' ...
0378 '44 2 2 0 4 8 22 19\n' ...
0379 '45 2 2 0 4 13 21 20\n' ...
0380 '46 2 2 0 4 15 23 17\n' ...
0381 '47 2 2 0 4 16 24 18\n' ...
0382 '48 2 2 0 4 6 25 17\n' ...
0383 '49 2 2 0 4 11 26 18\n' ...
0384 '50 2 2 0 4 12 20 26\n' ...
0385 '51 2 2 0 4 7 19 25\n' ...
0386 '$EndElements\n'];
0387 end
0388