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 material indices from eidors

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 material indices from eidors
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
0022 
0023 if isempty(name); 
0024    name = ['fwd_mdl based on ', vol_filename];
0025 end
0026 
0027 if nargin < 4
0028     stim_pattern=[];
0029 end
0030 
0031 if nargin<3 || isempty(eprefix); 
0032    eprefix = 'electrode-';
0033 end
0034 
0035 if nargin<5
0036    z_contact=0.01; % singular if z_contact=0
0037 end
0038 
0039 
0040 % Model Geometry
0041 [srf,vtx,fc,bc,simp,edg,mat_ind,phys_names] = gmsh_read_mesh(vol_filename);
0042 fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0043                              stim_pattern, eprefix, z_contact, fc,phys_names);
0044 mat_indices= mk_mat_indices( mat_ind);
0045 if isempty(srf)
0046    fwd_mdl.boundary = find_boundary(fwd_mdl);
0047 end
0048 
0049 fwd_mdl.mat_idx = mat_indices;
0050 
0051 % build fwd_model structure
0052 function fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0053     stim_pattern, eprefix, z_contact, fc, phys_names)
0054     mdl.nodes    = vtx;
0055     mdl.elems    = simp;
0056     mdl.boundary = srf;
0057     mdl.boundary_numbers=fc; 
0058     mdl.gnd_node = 1;
0059     mdl.name = name;
0060     
0061     % Model Stimulation
0062     if ~isempty(stim_pattern)
0063         mdl.stimulation= stim_pattern;
0064     end
0065     
0066     % Electrodes
0067     electrodes = find_elec(phys_names,eprefix,z_contact);
0068     if ~isempty(fieldnames(electrodes));
0069         mdl.electrode =     electrodes;
0070     end
0071     mdl.solve=          'eidors_default';
0072     mdl.jacobian=       'eidors_default';
0073     mdl.system_mat=     'eidors_default';
0074 
0075     fwd_mdl= eidors_obj('fwd_model', mdl);
0076 
0077 % Output cell array of indices into each material type
0078 %   array order is sorted by length of material type
0079 function mat_indices= mk_mat_indices( mat_ind);
0080     % find length of mat_indices
0081     % test example: mat_ind=[10 12 14 14 12 12 14 12];
0082     sort_mi= sort(mat_ind(:));
0083     find_mi= find( diff([-1e8;sort_mi]) );
0084     len_mi = diff([find_mi;length(sort_mi)+1]);
0085     [jnk,idxs]= sort(-len_mi); %reverse sort
0086     l_idxs= length(idxs);
0087     mat_indices= cell(1, l_idxs);
0088     for i= 1:l_idxs;
0089         mat_idx_i= sort_mi(find_mi(idxs(i)));
0090         mat_indices{i}= find(mat_ind == mat_idx_i);
0091     end
0092 
0093 % Assumes that electrodes are numbered starting at 1, with prefix provided
0094 function electrodes = find_elec(phys_names,prefix,z_contact)
0095 electrodes = struct();
0096 phys_elecs = find(arrayfun(@(x)strncmp(x.name,prefix,length(prefix)),phys_names));
0097 for i = 1:length(phys_elecs)
0098     cur_elec = arrayfun(@(x)strcmp(sprintf('%s%d',prefix,i),x.name),phys_names(phys_elecs));
0099     electrodes(i).nodes = unique(phys_names(phys_elecs(cur_elec)).nodes(:));
0100     electrodes(i).z_contact = z_contact;
0101 end