gmsh_read_mesh

[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names,entities] = 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. A very basic GMSH v4 reader is now
 included.

0001 function [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names,entities] = gmsh_read_mesh(filename)
0002 %[srf,vtx,fc,bc,simp,edg,mat_ind,phys_names,entities] = 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. A very basic GMSH v4 reader is now
0021 % included.
0022 
0023 % $Id: gmsh_read_mesh.m 6848 2024-05-04 22:20:11Z bgrychtol $
0024 % (C) 2009 Bartosz Sawicki. Licensed under GPL V2
0025 % Modified by James Snyder, Mark Campbell, Symon Stowe, Alistair Boyle,
0026 % Bartek Grychtol
0027 
0028 if ischar(filename) && strcmp(filename,'UNIT_TEST'); do_unit_test; return; end
0029 
0030 fid = fopen(filename,'r');
0031 content = textscan(fid,'%s','Delimiter','\n');
0032 content = content{:};
0033 fclose(fid);
0034 phys_names = [];
0035 entities = [];
0036 idx = find(startsWith(content,'$'));
0037 for i = 1:2:numel(idx)
0038    block = content( (idx(i)+1) : idx(i+1)-1 );
0039    switch content{idx(i)}
0040       case '$MeshFormat'
0041          gmshformat = parse_format(block);
0042       case '$Entities'
0043          entities= parse_entities( block, gmshformat );
0044       case '$Elements'
0045          elements= parse_elements( block, gmshformat );
0046       case '$Nodes'
0047          nodes= get_lines_with_nodes( block, gmshformat );
0048       case '$PhysicalNames'
0049          phys_names= parse_names( block, gmshformat );
0050    end
0051 end
0052 
0053 % look up nodes for each of phys_names
0054 if ~isempty(phys_names)
0055     if (gmshformat >= 4.0)
0056         assert(~isempty(entities), 'expected $Entities section (GMSH format 4)');
0057         for i = 1:length(phys_names)
0058             phys_tag = phys_names(i).tag;
0059             dim = phys_names(i).dim;
0060             tmpentities = find(arrayfun(@(x) any(x.phys_tag == phys_tag), entities));
0061             tags = cat(1,entities(tmpentities).entity_tag);
0062             idx = ismember(elements.entity_tag, tags) & elements.dim == dim;
0063             phys_names(i).nodes = elements.simp(idx,1:dim+1);
0064             elements.phys_tag(idx) = phys_tag;
0065         end
0066     else % GMSH 2.0 doesn't have Entities
0067         assert(isempty(entities), 'unexpected $Entities section (GMSH format 2)');
0068         for i = 1:length(phys_names)
0069             dim = phys_names(i).dim;
0070             tag = phys_names(i).tag;
0071             idx = elements.phys_tag == tag & elements.dim == dim;
0072             phys_names(i).nodes = elements.simp(idx,1:dim+1);
0073         end
0074     end
0075 end
0076 
0077 edg = [];
0078 bc = [];
0079 
0080 % Type 2: 3-node triangle
0081 tri = elements.type==2;
0082 % Type 4: 4-node tetrahedron
0083 tet = elements.type==4;
0084 
0085 % Select 2d vs 3d model by checking if Z is all the same
0086 if any(nodes(:,4) ~= nodes(1,4))
0087     vtx = nodes(:,2:4);
0088     simp = elements.simp(tet,:);
0089     srf = elements.simp(tri,1:3);
0090     mat_ind = elements.phys_tag(tet);
0091 else
0092     vtx = nodes(:,2:3);
0093     simp = elements.simp(tri,1:3);
0094     srf = [];
0095     mat_ind = elements.phys_tag(tri);
0096 end
0097 
0098 fc = elements.phys_tag(tri); % for 4.1, these are all zero, why?
0099 end
0100 
0101 function mat = get_lines_with_nodes( lines, gmshformat )
0102     
0103     
0104     switch floor(gmshformat)
0105     % Version 2 Line Format:
0106     % node-number x-coord y-coord z-coord
0107     % Version 4 Line Format: (not always like this)
0108     % node-number
0109     % x-coord y-coord z-coord
0110     case 2 
0111       n_rows = parse_rows(lines{1},gmshformat);
0112       %mat= sscanf(sprintf('%s ',lines{2:end}),'%f',[4,n_rows])';
0113       mat= textscan(sprintf('%s ',lines{2:end}),'%f');
0114       mat = reshape(mat{1},[4,n_rows])';
0115     case 4;
0116         n = sscanf(lines{1}, '%d')';
0117         n_block = n(1);
0118         n_nodes = n(2);
0119         mat = zeros(n_nodes,4);
0120         count = 2;
0121         while n_block > 0
0122             n_block = n_block - 1;
0123             tline = lines{count}; count = count + 1; % fgetl(fid);
0124             blk = sscanf(tline, '%d')';
0125             blk_nodes = blk(end); % n(2) for v4.0, n(4) for v4.1
0126             if blk_nodes == 0,  continue, end
0127             if gmshformat == 4.1 % v4.1: node tags first as a block, then node coordinates in a block
0128                node_tag =  textscan(sprintf('%s ',lines{count:count+blk_nodes-1}),'%f');
0129                node_tag = node_tag{1};
0130                count = count + blk_nodes;
0131 %                el = sscanf(sprintf('%s ',lines{count:count+blk_nodes-1}),'%f',[3 blk_nodes])';
0132                el = textscan(sprintf('%s ',lines{count:count+blk_nodes-1}),'%f');
0133                el = reshape(el{1},[3 blk_nodes])';
0134                count = count + blk_nodes;
0135                mat(node_tag,:) = [node_tag(:), el(:,1:end)];
0136             else % v4.0: node tag and coordinates on the same line
0137 %                data = sscanf(sprintf('%s ',lines{count:count+blk_nodes-1}),'%f',[4 blk_nodes])';
0138                data = textscan(sprintf('%s ',lines{count:count+blk_nodes-1}),'%f');
0139                data = reshape(data{1},[4 blk_nodes])';
0140                mat(data(:,1),:) = data;
0141                count = count + blk_nodes;               
0142             end
0143             
0144         end
0145         assert(n_block == 0, 'failed to consume all $Nodes blocks');
0146     otherwise; error('cant parse gmsh file of this format');
0147     end
0148 end
0149 
0150 function gmshformat = parse_format(lines)
0151    rawformat = sscanf(lines{1},'%f');
0152    gmshformat = rawformat(1);
0153 end
0154 
0155 function n_rows = parse_rows(tline, gmshformat)
0156    n_rows = sscanf(tline,'%d');
0157    switch floor(gmshformat)
0158      case 2; n_rows = n_rows(1);
0159      case 4; n_rows = n_rows(2);
0160      otherwise; error('cant parse gmsh file of this format');
0161    end
0162 end
0163 
0164 function names = parse_names( lines, version )
0165     % Line Format:
0166     % physical-dimension physical-number "physical-name"
0167     n_rows = sscanf(lines{1},'%d');
0168     names = struct('tag',{},'dim',{},'name',{});
0169     for i = 1:n_rows
0170         tline = lines{1+i};
0171         parts = regexp(tline,' ','split');
0172         nsz = size(names,2)+1;
0173         names(nsz).dim = str2double( parts(1) );
0174         names(nsz).tag = str2double( parts(2) );
0175         tname = parts(3);
0176         names(nsz).name = strrep(tname{1},'"','');
0177     end
0178 end % end function
0179 
0180 function elements = parse_entities( lines, gmshformat )
0181     elements = [];
0182     switch floor(gmshformat)
0183       case 2; warning('ignoring $Entities$ for GMSH v2 file');
0184       case 4;
0185           elements = parse_v4_entities(lines, gmshformat);
0186       otherwise error('cant parse this file type');
0187     end
0188  end
0189 
0190  function entities = parse_v4_entities(lines, gmshformat)
0191 % http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
0192 % $Entities
0193 %  numPoints(size_t) numCurves(size_t)
0194 %    numSurfaces(size_t) numVolumes(size_t)
0195 %  pointTag(int) X(double) Y(double) Z(double)
0196 %    numPhysicalTags(size_t) physicalTag(int) ...
0197 %  ...
0198 %  curveTag(int) minX(double) minY(double) minZ(double)
0199 %    maxX(double) maxY(double) maxZ(double)
0200 %    numPhysicalTags(size_t) physicalTag(int) ...
0201 %    numBoundingPoints(size_t) pointTag(int) ...
0202 %  ...
0203 %  surfaceTag(int) minX(double) minY(double) minZ(double)
0204 %    maxX(double) maxY(double) maxZ(double)
0205 %    numPhysicalTags(size_t) physicalTag(int) ...
0206 %    numBoundingCurves(size_t) curveTag(int) ...
0207 %  ...
0208 %  volumeTag(int) minX(double) minY(double) minZ(double)
0209 %    maxX(double) maxY(double) maxZ(double)
0210 %    numPhysicalTags(size_t) physicalTag(int) ...
0211 %    numBoundngSurfaces(size_t) surfaceTag(int) ...
0212 %  ...
0213 % $EndEntities
0214 % Entities are necessary to map '$PhysicalNames' to 'entityDim' & 'entityTag'
0215 % in $Elements and $Nodes.
0216 % Note that the entityTag can repeat for each type of entityDim, so when we
0217 % look up the elements associated with a PhysicalName we need to then use the
0218 % entityDim AND entityTag to find matching Nodes/Elements.
0219     entities = struct('phys_tag',{},'dim',{},'entity_tag',{});
0220 
0221     bl = sscanf(lines{1}, '%d')'; % Get the line info
0222     n_points = bl(1);
0223     n_curves = bl(2);
0224     n_surf = bl(3);
0225     n_vol = bl(4);
0226     %fprintf('entities: %d pts, %d curves, %d surf, %d vol\n', n_points, n_curves, n_surf, n_vol);
0227     for i = 2:numel(lines)
0228         tline = lines{i};
0229         bl = sscanf(tline, '%f')'; % Get the line info
0230         off = 8; % offset for 'numPhysicalTags'
0231         if n_points > 0
0232             n_points = n_points - 1;
0233             dim = 0; % point
0234             off = 5; % offset for 'numPhysicalTags'
0235         elseif n_curves > 0
0236             n_curves = n_curves - 1;
0237             dim = 1; % line/curve
0238         elseif n_surf > 0
0239             n_surf = n_surf - 1;
0240             dim = 2; % surface
0241         elseif n_vol > 0
0242             n_vol = n_vol - 1;
0243             dim = 3; % volume
0244         else
0245             error('extra $Entities found in v4.1 GMSH file');
0246         end
0247         if bl(off) > 0
0248             entities(end+1).phys_tag = bl((off+1):uint32(off+bl(off)));
0249             entities(end).entity_tag = uint32(bl(1));
0250             entities(end).dim = dim;
0251         end
0252     end
0253     assert(n_points == 0, 'missing Entity/points from GMSH 4.1 file');
0254     assert(n_curves == 0, 'missing Entity/curves from GMSH 4.1 file');
0255     assert(n_surf == 0, 'missing Entity/surfaces from GMSH 4.1 file');
0256     assert(n_vol == 0, 'missing Entity/volumes from GMSH 4.1 file');
0257 end
0258 
0259 
0260 function elements = parse_elements( lines, gmshformat )
0261    n_rows = parse_rows(lines{1},gmshformat);
0262    switch floor(gmshformat)
0263      case 2; elements = parse_v2_elements(lines(2:end),n_rows);
0264      case 4; elements = parse_v4_elements(lines,gmshformat);
0265      otherwise error('cant parse this file type');
0266    end
0267 end
0268 
0269 function elements = parse_v4_elements(lines,gmshformat)
0270 % http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
0271 % $Elements
0272 % numEntityBlocks numElements minElementTag maxElementTag
0273 % ...
0274 % entityDim entityTag elementType numElementsInBlock
0275 % elementTag nodeTag ... nodeTag
0276 % ...
0277 % $EndElements
0278 % An entity is a node, curve, surface or volume
0279 % Each entity has a list of contained elements and corresponding node tags
0280 % 0D - 1 node tag (ignore these)
0281 % 1D - 2 node tags
0282 % 2D - 3 node tags
0283 % 3D - 4 node tags
0284  
0285     bl = sscanf(lines{1}, '%d')'; % Get the line info
0286     n_blocks = bl(1);
0287     n_elems = bl(2);
0288     e_block = 0;
0289     simp = zeros([n_elems,4],'uint32');
0290     type = zeros([n_elems,1],'uint8');
0291     entity_tag = zeros([n_elems,1],'uint32');
0292     type = zeros([n_elems,1], 'uint8');
0293     dim = zeros([n_elems,1], 'uint8');
0294     phys_tag = zeros([n_elems,1], 'uint32'); 
0295     count = 2;
0296     for i = 1:n_blocks
0297         bl = sscanf(lines{count}, '%d')'; % Get the line info
0298         count = count + 1;
0299         if gmshformat >= 4.1
0300            e_dim = bl(1); % entityDim
0301            e_tag = bl(2); % entityTag
0302         else % 4.0
0303            e_tag = bl(1); % entityTag
0304            e_dim = bl(2); % entityDim
0305         end
0306         e_type = bl(3); % elementType
0307         e_block = bl(4); % numElementsInBlock: Size of entity block
0308 %         data = sscanf(sprintf('%s ',lines{count:count+e_block-1}),'%d',[e_dim+2 e_block])';
0309         data = textscan(sprintf('%s ',lines{count:count+e_block-1}),'%d');
0310         data = reshape(data{1},[e_dim+2 e_block])';
0311         count = count + e_block;
0312         n_elems = n_elems - size(data,1);
0313         idx = data(:,1);
0314         simp(idx,1:e_dim+1) = data(:,2:end);
0315         type(idx) = e_type;
0316         entity_tag(idx) = e_tag;
0317         dim(idx) = e_dim;
0318     end
0319     assert(n_elems == 0, 'missing Elements from GMSH 4 file');
0320 %     assert(n_blocks == 0, 'missing Element/blocks from GMSH 4 file');
0321     elements = struct('simp',simp,'type',type,'entity_tag',entity_tag,...
0322                       'dim',dim,'phys_tag',phys_tag);
0323 
0324 % Copied from somewhere else in the code, to reintroduce once we deal
0325 % with PhysTags again
0326 %     assert(length(elements(j).phys_tag) > 0, ...
0327 %        'GMSH format v4 mesh volume elements can only have one $PhysicalName when imported into EIDORS');
0328 end
0329 
0330 function elements = parse_v2_elements(lines,n_rows)
0331 % Line Format:
0332 % elm-number elm-type number-of-tags < tag > ... node-number-list
0333 elements.simp = zeros([n_rows,4],'uint32');;
0334 elements.phys_tag = zeros([n_rows,1], 'uint32');
0335 elements.geom_tag = zeros([n_rows,1],'uint32');;
0336 elements.type = zeros([n_rows,1],'uint8');
0337 elements.dim = zeros([n_rows,1], 'uint8');
0338 
0339 for i = 1:n_rows % can we get rid of this loop?
0340     tline = lines{i};
0341     n = sscanf(tline, '%d')';
0342     ndim = numel(n) - n(3) - 4;
0343     elements.simp(i,1:ndim+1) = n(n(3) + 4:end);
0344     elements.type(i) = n(2);
0345     elements.dim(i) = ndim;
0346     if n(3) > 0 % get tags if they exist
0347         % By default, first tag is number of parent physical entity
0348         % second is parent elementary geometrical entity
0349         % third is number of parent mesh partitions followed by
0350         % partition ids
0351         tags = n(4:3+n(3));
0352         if length(tags) >= 1
0353             elements.phys_tag(i) = tags(1);
0354             if length(tags) >= 2
0355                 elements.geom_tag(i) = tags(2);
0356             end
0357         end
0358     end
0359 end
0360 end
0361 
0362 function do_unit_test
0363    selfdir = fileparts(which('gmsh_read_mesh'));
0364    [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(fullfile(selfdir, 'tests/test-4.0.msh'));
0365     unit_test_cmp('v4 vtx ',vtx(2:3,:),[1,0;-1,0])
0366     unit_test_cmp('v4 simp',simp(2:3,:),[3,7,12; 12, 7,14]);
0367 
0368    [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(fullfile(selfdir, 'tests/test-2.2.msh'));
0369     unit_test_cmp('v2 vtx ',vtx(2:3,:),[1,0;-1,0])
0370     unit_test_cmp('v2 simp',simp(2:3,:),[2,4,15; 14,17,19]);
0371 
0372    vers = {'2.2', '4.0', '4.1'};
0373    for ver = vers(:)'
0374        ver = ver{1};
0375        [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = ...
0376            gmsh_read_mesh( fullfile(selfdir, ['tests/box-' ver '.msh']) );
0377        unit_test_cmp(['2d v' ver ' vtx '],vtx,[0,0;1,0;1,1;0,1;0.5,0.5])
0378        unit_test_cmp(['2d v' ver ' simp'],simp,[2,5,1;1,5,4;3,5,2;,4,5,3]);
0379        unit_test_cmp(['2d v' ver ' phys'],{phys_names(:).name},{'elec#1','elec#2','main'});
0380 
0381        [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = ...
0382            gmsh_read_mesh( fullfile(selfdir, ['tests/cube-' ver '.msh']) );
0383        unit_test_cmp(['3d v' ver ' vtx '],vtx([1,2,13,end],:),[0,0,1;0,0,0;,0.5,0.5,0;0.5,0.5,1])
0384        unit_test_cmp(['3d v' ver ' simp'],simp([1,2,23,end],:), ...
0385            [10,11,12,13;9,12,14,11;12,14,10,7;13,10,11,6]);
0386        unit_test_cmp(['3d v' ver ' phys'],{phys_names(:).name},{'elec#1','elec#2','main'});
0387    end
0388 
0389 end