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