Description of gmsh_write

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 
0014 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 is_img = (isfield(mdl, 'type') && strcmp(mdl.type, 'image')) || isfield(mdl, 'elem_data');
0017 is_imdl = (isfield(mdl, 'type') && strcmp(mdl.type, 'inv_model')) || isfield(mdl, 'fwd_model');
0018 
0019 if nargin < 3
0020     outfile = data;
0021     data = [];
0022 end
0023 
0024 if is_img && isempty(data)
0025     data = mdl.elem_data(:,1); % TODO only handles first frame of data at the moment
0026 end
0027 assert(isempty(data) || size(data,2) == 1, 'error: expect a single frame of data in gmsh_write_mesh');
0028 
0029 if is_imdl || is_img
0030     if isfield(mdl, 'rec_model')
0031         mdl = mdl.rec_model;
0032     else
0033         mdl = mdl.fwd_model;
0034     end
0035 end
0036 assert(isempty(data) || size(data,1) == size(mdl.elems,1), 'error: expect a data to match number of elements in gmsh_write_mesh');
0037 
0038 nn = size(mdl.nodes,1); % number of nodes
0039 ne = size(mdl.elems,1); % number of elements
0040 ndim = size(mdl.nodes,2); % number of dimensions: 2=2D or 3=3D
0041 assert((ndim >= 2) && (ndim <= 3), 'not 2D or 3D?!');
0042 
0043 fid = fopen(outfile, 'w');
0044 fprintf(fid, '$MeshFormat\n');
0045 fprintf(fid, '4.1 0 8\n'); % version file-type(0 for ASCII, 1 for binary) data-size=sizeof(size_t)
0046 fprintf(fid, '$EndMeshFormat\n');
0047 fprintf(fid, '$Nodes\n');
0048 fprintf(fid, '%d %d %d %d\n', 1, nn, 1, nn); % numEntityBlocks numNodes minNodeTag maxNodeTag
0049 fprintf(fid, '%d %d %d %d\n', ndim, 0, 0, nn); % entityDim entityTag parametric(0 or 1) numNodesInBlock
0050 fprintf(fid, '%d\n', 1:nn); % nodeTag
0051 nodes = [ mdl.nodes zeros(nn, 3 - ndim) ]; % always nn x 3, with z=0.0 if 2D
0052 fprintf(fid, '%f %f %f\n', nodes'); % x y z
0053 fprintf(fid, '$EndNodes\n');
0054 fprintf(fid, '$Elements\n');
0055 if ndim == 2
0056     type = 2; % triangle
0057 else
0058     type = 4; % tetrahedra
0059 end
0060 fprintf(fid, '%d %d %d %d\n', 1, ne, 1, ne); % numEntityBlocks numElements minElementTag maxElementTag
0061 fprintf(fid, '%d %d %d %d\n', ndim, 0, type, ne); % entityDim entityTag elementType(4=tet) numElementsInBlock
0062 elems = [[1:ne]' mdl.elems];
0063 if ndim == 2
0064     fprintf(fid, '%d %d %d %d\n', elems'); % elementTag nodeTag ...
0065 else
0066     fprintf(fid, '%d %d %d %d %d\n', elems'); % elementTag nodeTag ...
0067 end
0068 fprintf(fid, '$EndElements\n');
0069 if ~isempty(data)
0070     fprintf(fid, '$ElementData\n');
0071     fprintf(fid, '%d\n', 1); % numStringTags
0072     name = 'gmsh_write_mesh output';
0073     if isfield(mdl, 'name')
0074         name = mdl.name;
0075     end
0076     fprintf(fid, '%s\n', name); % stringTag
0077     fprintf(fid, '%d\n', 1); % numRealTags
0078     fprintf(fid, '%f\n', 0); % realTag(time=0)
0079     fprintf(fid, '%d\n', 3); % numIntegerTags
0080     fprintf(fid, '%d\n', 0); % integerTag(time-step=0)
0081     fprintf(fid, '%d\n', 1); % integerTag(1-component/scalar field)
0082     fprintf(fid, '%d\n', ne); % integerTag(associated element values)
0083     data = [[1:ne]' data];
0084     fprintf(fid, '%d %f\n', data'); % integerTag(associated element values)
0085     fprintf(fid, '$EndElementData\n');
0086 end
0087 fclose(fid);
0088 
0089 function do_unit_test()
0090 % TODO cannot check 'data' because gmsh_read_mesh and gmsh_mk_fwd_model
0091 %      don't know about ElementData regions in the .msh file
0092 outfile = tempname();
0093 selfdir = fileparts(which('gmsh_read_mesh'));
0094 mdl1 = gmsh_mk_fwd_model(fullfile(selfdir, 'cube-4.1.msh'));
0095 mdl2 = gmsh_mk_fwd_model(fullfile(selfdir, 'box-4.1.msh'));
0096 data1 = rand(size(mdl1.elems,1), 1);
0097 data2 = rand(size(mdl1.elems,1), 1);
0098 data3 = rand(size(mdl2.elems,1), 1);
0099 
0100 gmsh_write_mesh(mdl1, data1, outfile);
0101 mdl = gmsh_mk_fwd_model(outfile);
0102 delete(outfile);
0103 unit_test_cmp('fmdl+data (nn)', mdl1.nodes, mdl.nodes);
0104 unit_test_cmp('fmdl+data (el)', mdl1.elems, mdl.elems);
0105 % TODO check 'data' matches data1
0106 
0107 gmsh_write_mesh(mdl1, outfile);
0108 mdl = gmsh_mk_fwd_model(outfile);
0109 delete(outfile);
0110 unit_test_cmp('fmdl (nn)', mdl1.nodes, mdl.nodes);
0111 unit_test_cmp('fmdl (ee)', mdl1.elems, mdl.elems);
0112 % TODO check 'data' does not exist
0113 
0114 img = mk_image(mdl1, data2);
0115 gmsh_write_mesh(img, data1, outfile);
0116 mdl = gmsh_mk_fwd_model(outfile);
0117 delete(outfile);
0118 unit_test_cmp('img+data (nn)', mdl1.nodes, mdl.nodes);
0119 unit_test_cmp('img+data (ee)', mdl1.elems, mdl.elems);
0120 % TODO check 'data' matches data1
0121 
0122 img = mk_image(mdl1, data2);
0123 gmsh_write_mesh(img, outfile);
0124 mdl = gmsh_mk_fwd_model(outfile);
0125 delete(outfile);
0126 unit_test_cmp('img (nn)', mdl1.nodes, mdl.nodes);
0127 unit_test_cmp('img (ee)', mdl1.elems, mdl.elems);
0128 % TODO check 'data' matches data2
0129 
0130 imdl = select_imdl(mdl1, {'Basic GN dif'});
0131 
0132 gmsh_write_mesh(imdl, outfile);
0133 mdl = gmsh_mk_fwd_model(outfile);
0134 delete(outfile);
0135 unit_test_cmp('imdl+fmdl (nn)', mdl1.nodes, mdl.nodes);
0136 unit_test_cmp('imdl+fmdl (ee)', mdl1.elems, mdl.elems);
0137 % TODO check 'data' does not exist
0138 
0139 gmsh_write_mesh(imdl, data1, outfile);
0140 mdl = gmsh_mk_fwd_model(outfile);
0141 delete(outfile);
0142 unit_test_cmp('imdl+fmdl+data (nn)', mdl1.nodes, mdl.nodes);
0143 unit_test_cmp('imdl+fmdl+data (ee)', mdl1.elems, mdl.elems);
0144 % TODO check 'data' matches data1
0145 
0146 imdl.rec_model = mdl2;
0147 
0148 gmsh_write_mesh(imdl, outfile);
0149 mdl = gmsh_mk_fwd_model(outfile);
0150 delete(outfile);
0151 unit_test_cmp('imdl+rmdl (nn)', mdl2.nodes, mdl.nodes);
0152 unit_test_cmp('imdl+rmdl (ee)', mdl2.elems, mdl.elems);
0153 % TODO check 'data' does not exist
0154 
0155 gmsh_write_mesh(imdl, data3, outfile);
0156 mdl = gmsh_mk_fwd_model(outfile);
0157 delete(outfile);
0158 unit_test_cmp('imdl+rmdl+data (nn)', mdl2.nodes, mdl.nodes);
0159 unit_test_cmp('imdl+rmdl+data (ee)', mdl2.elems, mdl.elems);
0160 % TODO check 'data' matches data3
gmsh_write_mesh

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
