gmsh_write_mesh

gmsh_read_mesh(mdl, [data,] outfile)
 Write a GMSH .msh file (v4.1 format)

 mdl      An EIDORS model (fwd_model or rec_model) or an image (img)
 data     Per element data (optional)
 outfile  Output filename, should end in '.msh'

 If 'mdl' is an EIDORS image, then img.elem_data will be used if 'data' is not
 provided.

0001 function gmsh_write_mesh(mdl, data, outfile)
0002 %gmsh_read_mesh(mdl, [data,] outfile)
0003 % Write a GMSH .msh file (v4.1 format)
0004 %
0005 % mdl      An EIDORS model (fwd_model or rec_model) or an image (img)
0006 % data     Per element data (optional)
0007 % outfile  Output filename, should end in '.msh'
0008 %
0009 % If 'mdl' is an EIDORS image, then img.elem_data will be used if 'data' is not
0010 % provided.
0011 
0012 % (C) 2021 Alistair Boyle. Licensed under GPL v2 or v3
0013 % Modified by Bartek Grychtol
0014 % $Id: gmsh_write_mesh.m 7102 2024-12-26 12:12:18Z aadler $
0015 
0016 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0017 
0018 is_img = (isfield(mdl, 'type') && strcmp(mdl.type, 'image')) || isfield(mdl, 'elem_data');
0019 is_imdl = (isfield(mdl, 'type') && strcmp(mdl.type, 'inv_model')) || isfield(mdl, 'fwd_model');
0020 
0021 if nargin < 3
0022     outfile = data;
0023     data = [];
0024 end
0025 
0026 if is_img && isempty(data)
0027     data = mdl.elem_data(:,1); % TODO only handles first frame of data at the moment
0028 end
0029 assert(isempty(data) || size(data,2) == 1, 'error: expect a single frame of data in gmsh_write_mesh');
0030 
0031 if is_imdl || is_img
0032     if isfield(mdl, 'rec_model')
0033         mdl = mdl.rec_model;
0034     else
0035         mdl = mdl.fwd_model;
0036     end
0037 end
0038 assert(isempty(data) || size(data,1) == size(mdl.elems,1), 'error: expect a data to match number of elements in gmsh_write_mesh');
0039 
0040 nn = size(mdl.nodes,1); % number of nodes
0041 ne = size(mdl.elems,1); % number of elements
0042 ndim = size(mdl.nodes,2); % number of dimensions: 2=2D or 3=3D
0043 assert((ndim >= 2) && (ndim <= 3), 'not 2D or 3D?!');
0044 
0045 fid = fopen(outfile, 'w');
0046 fprintf(fid, '$MeshFormat\n');
0047 fprintf(fid, '4.1 0 8\n'); % version file-type(0 for ASCII, 1 for binary) data-size=sizeof(size_t)
0048 fprintf(fid, '$EndMeshFormat\n');
0049 
0050 nodes = [ mdl.nodes zeros(nn, 3 - ndim) ]; % always nn x 3, with z=0.0 if 2D
0051 
0052 
0053 [mdl, mat_vec] = prepare_materials(mdl, ne);
0054 numEntityBlocks = 1;
0055 try numEntityBlocks = numel(mdl.mat_idx); end
0056 numPhysicalNames = 0;
0057 try numPhysicalNames = numel(mdl.mat_idx); end
0058 
0059 if numPhysicalNames
0060    fprintf(fid, '$PhysicalNames\n');
0061    fprintf(fid, '%d\n',numPhysicalNames);
0062    for i = 1:numPhysicalNames
0063       fprintf(fid, '%d %d "%s"\n', ndim, i, mdl.mat_names{i}); % dimension physicalTag "name"
0064    end
0065    fprintf(fid, '$EndPhysicalNames\n');
0066 end
0067 
0068 fprintf(fid, '$Entities\n');
0069 if ndim == 2
0070    fprintf(fid, '%d %d %d %d\n', 0, 0, numEntityBlocks, 0); % numPoints numCurves numSurfaces numVolumes
0071 else
0072    fprintf(fid, '%d %d %d %d\n', 0, 0, 0, numEntityBlocks); % numPoints numCurves numSurfaces numVolumes
0073 end
0074 for i = 1:numEntityBlocks
0075    physicalTag = 0; 
0076    if numPhysicalNames, physicalTag = i; end
0077    if numEntityBlocks == 1
0078       idx = (1:ne)';
0079    else 
0080       idx = reshape(mdl.mat_idx{i},[],1);
0081    end
0082    minlim = min(nodes(unique(mdl.elems(idx,:)), :));
0083    maxlim = max(nodes(unique(mdl.elems(idx,:)), :));
0084    fprintf(fid, '%d %f %f %f %f %f %f %d %d %d\n', ...  
0085             i, ... % volumeTag
0086             minlim, ... % minX minY minZ 
0087             maxlim, ... % maxX maxY maxZ
0088             1, physicalTag, 0); % numPhysicalTags physicalTag numBoundingSurfaces
0089 end
0090 fprintf(fid, '$EndEntities\n');
0091 
0092 fprintf(fid, '$Nodes\n');
0093 fprintf(fid, '%d %d %d %d\n', 1, nn, 1, nn); % numEntityBlocks numNodes minNodeTag maxNodeTag
0094 fprintf(fid, '%d %d %d %d\n', ndim, 0, 0, nn); % entityDim entityTag parametric(0 or 1) numNodesInBlock
0095 fprintf(fid, '%d\n', 1:nn); % nodeTag
0096 fprintf(fid, '%f %f %f\n', nodes'); % x y z
0097 fprintf(fid, '$EndNodes\n');
0098 
0099 fprintf(fid, '$Elements\n');
0100 if ndim == 2
0101     type = 2; % triangle
0102 else
0103     type = 4; % tetrahedra
0104 end
0105 fprintf(fid, '%d %d %d %d\n', numEntityBlocks, ne, 1, ne); % numEntityBlocks numElements minElementTag maxElementTag
0106 
0107 for i = 1:numEntityBlocks
0108    if isempty(mat_vec)
0109       neb = ne;
0110       idx = (1:ne)';
0111    else
0112       neb = numel(mdl.mat_idx{i});
0113       idx = mdl.mat_idx{i};
0114    end
0115    fprintf(fid, '%d %d %d %d\n', ndim, i, type, neb); % entityDim entityTag elementType(4=tet) numElementsInBlock
0116    elems = [idx mdl.elems(idx,:)];
0117    if ndim == 2
0118       fprintf(fid, '%d %d %d %d\n', elems'); % elementTag nodeTag ...
0119    else
0120       fprintf(fid, '%d %d %d %d %d\n', elems'); % elementTag nodeTag ...
0121    end
0122 end
0123 fprintf(fid, '$EndElements\n');
0124 
0125 if ~isempty(data)
0126     fprintf(fid, '$ElementData\n');
0127     fprintf(fid, '%d\n', 1); % numStringTags
0128     name = 'gmsh_write_mesh output';
0129     if isfield(mdl, 'name')
0130         name = mdl.name;
0131     end
0132     fprintf(fid, '%s\n', name); % stringTag
0133     fprintf(fid, '%d\n', 1); % numRealTags
0134     fprintf(fid, '%f\n', 0); % realTag(time=0)
0135     fprintf(fid, '%d\n', 3); % numIntegerTags
0136     fprintf(fid, '%d\n', 0); % integerTag(time-step=0)
0137     fprintf(fid, '%d\n', 1); % integerTag(1-component/scalar field)
0138     fprintf(fid, '%d\n', ne); % integerTag(associated element values)
0139     data = [[1:ne]' data];
0140     fprintf(fid, '%d %f\n', data'); % integerTag(associated element values)
0141     fprintf(fid, '$EndElementData\n');
0142 end
0143 fclose(fid);
0144 
0145 function [mdl, mat_vec] = prepare_materials(mdl, ne)
0146 mat_vec = [];
0147 if isfield(mdl, 'mat_idx') && ~isempty(mdl.mat_idx)
0148    mat_vec = zeros(ne,1);
0149    for i = 1:numel(mdl.mat_idx)
0150       idx = mdl.mat_idx{i};
0151       if any(mat_vec(idx))
0152          error('Materials are not unique');
0153       end
0154       mat_vec(idx) = i;
0155    end
0156    if isfield(mdl,'mat_names') && ~isempty(mdl.mat_names)
0157       assert(numel(mdl.mat_names) == numel(mdl.mat_idx),'Every mat_idx must have a name!')
0158    end
0159    if any(mat_vec==0)
0160       mdl.mat_idx(end+1) = {find(mat_vec==1)};
0161       if isfield(mdl,'mat_names')
0162          mdl.mat_names(end+1) = {'Background'};
0163       end
0164    end
0165 end
0166 
0167 
0168 function do_unit_test()
0169 % TODO cannot check 'data' because gmsh_read_mesh and gmsh_mk_fwd_model
0170 %      don't know about ElementData regions in the .msh file
0171 outfile = tempname();
0172 selfdir = fileparts(which('gmsh_read_mesh'));
0173 mdl1 = gmsh_mk_fwd_model(fullfile(selfdir, 'tests/cube-4.1.msh'));
0174 mdl2 = gmsh_mk_fwd_model(fullfile(selfdir, 'tests/box-4.1.msh'));
0175 data1 = rand(size(mdl1.elems,1), 1);
0176 data2 = rand(size(mdl1.elems,1), 1);
0177 data3 = rand(size(mdl2.elems,1), 1);
0178 
0179 gmsh_write_mesh(mdl1, data1, outfile);
0180 mdl = gmsh_mk_fwd_model(outfile);
0181 delete(outfile);
0182 unit_test_cmp('fmdl+data (nn)', mdl1.nodes, mdl.nodes);
0183 unit_test_cmp('fmdl+data (el)', mdl1.elems, mdl.elems);
0184 % TODO check 'data' matches data1
0185 
0186 gmsh_write_mesh(mdl1, outfile);
0187 mdl = gmsh_mk_fwd_model(outfile);
0188 delete(outfile);
0189 unit_test_cmp('fmdl (nn)', mdl1.nodes, mdl.nodes);
0190 unit_test_cmp('fmdl (ee)', mdl1.elems, mdl.elems);
0191 % TODO check 'data' does not exist
0192 
0193 img = mk_image(mdl1, data2);
0194 gmsh_write_mesh(img, data1, outfile);
0195 mdl = gmsh_mk_fwd_model(outfile);
0196 delete(outfile);
0197 unit_test_cmp('img+data (nn)', mdl1.nodes, mdl.nodes);
0198 unit_test_cmp('img+data (ee)', mdl1.elems, mdl.elems);
0199 % TODO check 'data' matches data1
0200 
0201 img = mk_image(mdl1, data2);
0202 gmsh_write_mesh(img, outfile);
0203 mdl = gmsh_mk_fwd_model(outfile);
0204 delete(outfile);
0205 unit_test_cmp('img (nn)', mdl1.nodes, mdl.nodes);
0206 unit_test_cmp('img (ee)', mdl1.elems, mdl.elems);
0207 % TODO check 'data' matches data2
0208 
0209 imdl = select_imdl(mdl1, {'Basic GN dif'});
0210 
0211 gmsh_write_mesh(imdl, outfile);
0212 mdl = gmsh_mk_fwd_model(outfile);
0213 delete(outfile);
0214 unit_test_cmp('imdl+fmdl (nn)', mdl1.nodes, mdl.nodes);
0215 unit_test_cmp('imdl+fmdl (ee)', mdl1.elems, mdl.elems);
0216 % TODO check 'data' does not exist
0217 
0218 gmsh_write_mesh(imdl, data1, outfile);
0219 mdl = gmsh_mk_fwd_model(outfile);
0220 delete(outfile);
0221 unit_test_cmp('imdl+fmdl+data (nn)', mdl1.nodes, mdl.nodes);
0222 unit_test_cmp('imdl+fmdl+data (ee)', mdl1.elems, mdl.elems);
0223 % TODO check 'data' matches data1
0224 
0225 imdl.rec_model = mdl2;
0226 
0227 gmsh_write_mesh(imdl, outfile);
0228 mdl = gmsh_mk_fwd_model(outfile);
0229 delete(outfile);
0230 unit_test_cmp('imdl+rmdl (nn)', mdl2.nodes, mdl.nodes);
0231 unit_test_cmp('imdl+rmdl (ee)', mdl2.elems, mdl.elems);
0232 % TODO check 'data' does not exist
0233 
0234 gmsh_write_mesh(imdl, data3, outfile);
0235 mdl = gmsh_mk_fwd_model(outfile);
0236 delete(outfile);
0237 unit_test_cmp('imdl+rmdl+data (nn)', mdl2.nodes, mdl.nodes);
0238 unit_test_cmp('imdl+rmdl+data (ee)', mdl2.elems, mdl.elems);
0239 % TODO check 'data' matches data3