gmsh_mk_fwd_model

 GMSH_MK_FWD_MODEL: create a fwd_model object from a Gmsh file
 [fwd_mdl, mat_indices]= ...
      gmsh_mk_fwd_model( vol_filename, centres, ...
                       name, stim_pattern, z_contact)

  vol_filename:      filename output from Gmsh
  name:              name for object (if [] use vol_filename)
  eprefix:           prefix used for names of electrodes
                     (if [] or omitted use 'electrode-')
  stim_pattern:      a stimulation pattern structure
                     empty ([]) if stim_pattern is not available
  z_contact:         vector or scalar electrode contact impedance

  fwd_mdl:           eidors format fwd_model
  mat_indices:       cell array of element indices, per material

0001 function [fwd_mdl, mat_indices]= ...
0002              gmsh_mk_fwd_model( vol_filename, name, eprefix,stim_pattern, z_contact)
0003 % GMSH_MK_FWD_MODEL: create a fwd_model object from a Gmsh file
0004 % [fwd_mdl, mat_indices]= ...
0005 %      gmsh_mk_fwd_model( vol_filename, centres, ...
0006 %                       name, stim_pattern, z_contact)
0007 %
0008 %  vol_filename:      filename output from Gmsh
0009 %  name:              name for object (if [] use vol_filename)
0010 %  eprefix:           prefix used for names of electrodes
0011 %                     (if [] or omitted use 'electrode-')
0012 %  stim_pattern:      a stimulation pattern structure
0013 %                     empty ([]) if stim_pattern is not available
0014 %  z_contact:         vector or scalar electrode contact impedance
0015 %
0016 %  fwd_mdl:           eidors format fwd_model
0017 %  mat_indices:       cell array of element indices, per material
0018 
0019 % Gmsh mesher for EIDORS was based on Netgen interface.
0020 % (C) 2009 Bartosz Sawicki. License: GPL version 2 or version 3
0021 % Modified by James Snyder, Bartlomiej Grychtol, Alistair Boyle
0022 
0023 if ischar(vol_filename) && strcmp(vol_filename,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 if nargin<2;
0026    name = vol_filename;
0027 end
0028 
0029 if nargin<3 || isempty(eprefix);
0030    eprefix = 'electrode-';
0031 end
0032 
0033 if nargin < 4
0034     stim_pattern=[];
0035 end
0036 
0037 if nargin<5
0038    z_contact=0.01; % singular if z_contact=0
0039 end
0040 
0041 
0042 % Model Geometry
0043 [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(vol_filename);
0044 fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0045                              stim_pattern, eprefix, z_contact, fc,phys_names);
0046 [mat_indices,mat_names] = mk_mat_indices( mat_ind, phys_names);
0047 if isempty(srf)
0048    fwd_mdl.boundary = find_boundary(fwd_mdl);
0049 end
0050 
0051 fwd_mdl.mat_idx = mat_indices;
0052 if length(mat_names) > 0
0053     fwd_mdl.mat_names = mat_names;
0054 end
0055 
0056 % build fwd_model structure
0057 function fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0058     stim_pattern, eprefix, z_contact, fc, phys_names)
0059     mdl.nodes    = vtx;
0060     mdl.elems    = simp;
0061     mdl.boundary = srf;
0062     mdl.boundary_numbers=fc; % TODO this is not very useful without mapping to phys_names, like mat_idx/mat_names
0063     mdl.gnd_node = 1;
0064     mdl.name = name;
0065 
0066     % Model Stimulation
0067     if ~isempty(stim_pattern)
0068         mdl.stimulation= stim_pattern;
0069     end
0070 
0071     % Electrodes
0072     electrodes = find_elec(phys_names,eprefix,z_contact);
0073     if ~isempty(fieldnames(electrodes));
0074         mdl.electrode =     electrodes;
0075     end
0076     mdl.solve=          'eidors_default';
0077     mdl.jacobian=       'eidors_default';
0078     mdl.system_mat=     'eidors_default';
0079 
0080     fwd_mdl= eidors_obj('fwd_model', mdl);
0081 
0082 % Output cell array of indices into each material type
0083 %   array order is sorted by length of material type
0084 function [mat_indices,mat_names]= mk_mat_indices(mat_ind,phys_names);
0085     % find length of mat_indices
0086     % test example: mat_ind=[10 12 14 14 12 12 14 12];
0087     sort_mi= sort(mat_ind(:));
0088     find_mi= find( diff([-1e8;sort_mi]) );
0089     len_mi = diff([find_mi;length(sort_mi)+1]);
0090     [jnk,idxs]= sort(-len_mi); %reverse sort
0091     l_idxs= length(idxs);
0092     mat_indices= cell(1, l_idxs);
0093     for i= 1:l_idxs;
0094         mat_idx_i= sort_mi(find_mi(idxs(i)));
0095         mat_indices{i}= find(mat_ind == mat_idx_i);
0096     end
0097     mat_names = {};
0098     if length(phys_names) > 0
0099         phys_dims = max(cat(1,phys_names(:).dim));
0100         dim = max(phys_dims);
0101         mat_names = cell(1, l_idxs);
0102         for i = 1:l_idxs
0103             tag = sort_mi(find_mi(idxs(i)));
0104             idx = find(arrayfun(@(x) and((x.tag == tag), (x.dim == dim)), phys_names));
0105             assert(length(idx) == 1, 'missing physical name for tag');
0106             mat_names{i} = phys_names(idx).name;
0107         end
0108     end
0109 
0110 % Assumes that electrodes are numbered starting at 1, with prefix provided
0111 function electrodes = find_elec(phys_names,prefix,z_contact)
0112 electrodes = struct();
0113 phys_elecs = find(arrayfun(@(x)strncmp(x.name,prefix,length(prefix)),phys_names));
0114 n_prefix = length(prefix);
0115 for i = 1:length(phys_elecs)
0116     cur_elec = arrayfun(@(x) str2double(x.name((n_prefix+1):end)) == i,phys_names(phys_elecs));
0117     electrodes(i).nodes = unique(phys_names(phys_elecs(cur_elec)).nodes(:));
0118     electrodes(i).z_contact = z_contact;
0119 end
0120 
0121 function do_unit_test
0122 DEBUG = 0; % enable if tests fail to help diagnose the problem
0123 
0124 selfdir = fileparts(which('gmsh_read_mesh'));
0125 vers = {'2.2', '4.0', '4.1'};
0126 for ver = vers(:)'
0127     ver = ver{1};
0128 
0129     % Expected forward model:
0130     stim = 'asdf';
0131     elec.nodes = [];
0132     elec.z_contact = 0.11;
0133     fmdl2d.nodes = [0,0;1,0;1,1;0,1;0.5,0.5];
0134     fmdl2d.elems = [2,5,1;1,5,4;3,5,2;,4,5,3];
0135     fmdl2d.boundary = uint32([1,2;1,4;2,3;3,4]);
0136     fmdl2d.boundary_numbers = ones(4,1)*5;
0137     fmdl2d.gnd_node = 1;
0138     fmdl2d.name = 'test-name';
0139     fmdl2d.stimulation = stim;
0140     fmdl2d.solve = 'eidors_default';
0141     fmdl2d.jacobian = 'eidors_default';
0142     fmdl2d.system_mat = 'eidors_default';
0143     fmdl2d.electrode(1) = elec;
0144     fmdl2d.electrode(2) = elec;
0145     fmdl2d.electrode(1).nodes = [1,4]';
0146     fmdl2d.electrode(2).nodes = [2,3]';
0147     fmdl2d.type = 'fwd_model';
0148     fmdl2d.mat_idx = {{1:4}};
0149     fmdl2d.mat_names = {'main'};
0150 
0151     fmdl3d = fmdl2d;
0152     fmdl3d.electrode(1).nodes = [1,2,3,4,9]';
0153     fmdl3d.electrode(2).nodes = [5,6,7,8,10]';
0154     fmdl3d.mat_idx = {{1:24}};
0155     fmdl3d.boundary_numbers = [ones(4,1)*13; ones(4,1)*14];
0156     fmdl3d.boundary = [
0157          2     1     9
0158          1     3     9
0159          4     2     9
0160          3     4     9
0161          6    10     5
0162          5    10     7
0163          8    10     6
0164          7    10     8 ];
0165     fmdl3d.elems = [
0166         10    11    12    13
0167          9    12    14    11
0168         12    14    11    10
0169          9    12    11    13
0170          2     9     1    11
0171          1     9     3    14
0172         11    14     1     5
0173          4     9    12     3
0174          2     4     9    13
0175         12     3    14     7
0176          5    10    14     7
0177          7    10    12     8
0178         11    10     5     6
0179         12     4     8    13
0180         13     8    10     6
0181         13    11     2     6
0182         14     1     9    11
0183          9    12     3    14
0184         11     2     9    13
0185         12     9     4    13
0186          8    10    12    13
0187         14    11    10     5
0188         12    14    10     7
0189         13    10    11     6 ];
0190     fmdl3d.nodes = [
0191              0         0    1.0000
0192              0         0         0
0193              0    1.0000    1.0000
0194              0    1.0000         0
0195         1.0000         0    1.0000
0196         1.0000         0         0
0197         1.0000    1.0000    1.0000
0198         1.0000    1.0000         0
0199              0    0.5000    0.5000
0200         1.0000    0.5000    0.5000
0201         0.5000         0    0.5000
0202         0.5000    1.0000    0.5000
0203         0.5000    0.5000         0
0204         0.5000    0.5000    1.0000 ];
0205 
0206     % Test 2D and 3D parsing
0207     [fmdl,mat_ind] = gmsh_mk_fwd_model( ...
0208         fullfile(selfdir, ['box-' ver '.msh']), fmdl2d.name, 'elec#', stim, elec.z_contact );
0209     if DEBUG
0210         for x = fields(fmdl)'
0211             x = x{1};
0212             unit_test_cmp(['2d v' ver ' fmdl.' x],fmdl.(x),fmdl2d.(x))
0213         end
0214         fmdl
0215         fmdl2d
0216     end
0217     unit_test_cmp(['2d v' ver ' fmdl'],fmdl,fmdl2d)
0218 
0219     [fmdl,mat_ind] = gmsh_mk_fwd_model( ...
0220         fullfile(selfdir, ['cube-' ver '.msh']), fmdl3d.name, 'elec#', stim, elec.z_contact );
0221     if DEBUG
0222         for x = fields(fmdl)'
0223             x = x{1};
0224             unit_test_cmp(['3d v' ver ' fmdl.' x],fmdl.(x),fmdl3d.(x))
0225         end
0226         fmdl
0227         fmdl3d
0228     end
0229     unit_test_cmp(['3d v' ver ' fmdl'],fmdl,fmdl3d)
0230 end