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