ng_mk_fwd_model

 NG_MK_FWD_MODEL: create a fwd_model object from a netgen vol file
 [fwd_mdl]= ...
      ng_mk_fwd_model( ng_vol_filename, centres, ...
                       name, stim_pattern, z_contact)

  ng_vol_filename:   filename output from netgen
  name:              name for object (if [] use ng_vol_filename)
  centres:           matrix of N x [x,y,z] electrode centres
                     centres can also be a Nx1 cell matrix of
                     functions which are 1 inside the electrode and 0 outside
  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
  fwd_mdl.mat_idx:   cell array of indices into netgen material properties
 (C) 2006 Andy Adler. (C) 2013 Alistair Boyle. License: GPL version 2 or version 3
 $Id: ng_mk_fwd_model.m 6355 2022-05-01 16:26:22Z aadler $

0001 function [fwd_mdl]= ...
0002              ng_mk_fwd_model( ng_vol_filename, centres, ...
0003                               name, stim_pattern, z_contact, postprocmesh)
0004 % NG_MK_FWD_MODEL: create a fwd_model object from a netgen vol file
0005 % [fwd_mdl]= ...
0006 %      ng_mk_fwd_model( ng_vol_filename, centres, ...
0007 %                       name, stim_pattern, z_contact)
0008 %
0009 %  ng_vol_filename:   filename output from netgen
0010 %  name:              name for object (if [] use ng_vol_filename)
0011 %  centres:           matrix of N x [x,y,z] electrode centres
0012 %                     centres can also be a Nx1 cell matrix of
0013 %                     functions which are 1 inside the electrode and 0 outside
0014 %  stim_pattern:      a stimulation pattern structure
0015 %                     empty ([]) if stim_pattern is not available
0016 %  z_contact:         vector or scalar electrode contact impedance
0017 %
0018 %  fwd_mdl:           eidors format fwd_model
0019 %  fwd_mdl.mat_idx:   cell array of indices into netgen material properties
0020 % (C) 2006 Andy Adler. (C) 2013 Alistair Boyle. License: GPL version 2 or version 3
0021 % $Id: ng_mk_fwd_model.m 6355 2022-05-01 16:26:22Z aadler $
0022 
0023 if nargin==1 && strcmp(ng_vol_filename,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 if isempty(name);
0026    name = ['MDL from', ng_vol_filename];
0027 end
0028 
0029 if nargin<5
0030    z_contact=0.01; % singular if z_contact=0
0031 end
0032 
0033 % Model Geometry
0034 [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_read_mesh(ng_vol_filename);
0035 if isempty(vtx);
0036    error('EIDORS: ng_mk_fwd_model: Netgen meshing failed. Stopping');
0037 end
0038 if nargin>=6
0039     N_elec = size(centres,1);
0040     [srf,vtx,fc,bc,simp,edg,mat_ind] = feval(postprocmesh,...
0041         srf,vtx,fc,bc,simp,edg,mat_ind, N_elec);
0042 end
0043 fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0044                              stim_pattern, centres, z_contact,fc);
0045 
0046 [fwd_mdl.mat_idx] = mk_mat_indices(mat_ind);
0047 
0048 if ~isfield(fwd_mdl,'normalize_measurements')
0049    fwd_mdl.normalize_measurements = 0;
0050 end
0051 
0052 % build fwd_model structure
0053 function fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0054                        stim_pattern, centres, z_contact,fc)
0055 % maybe there were bugs from an unsorted boundary. Just in case
0056 
0057 % get rid of duplicates, if there are any
0058 [~,idx] = unique(sort(srf,2),'rows','stable');
0059 if length(idx)~=size(srf,1)
0060    error('Duplicate elements on boundary. Netgen completed, but EIDORS can''t interpret the result. Try modifying the netgen parameters.')
0061 end
0062 %fc = fc(idx);
0063 %bc = bc(idx);
0064 
0065 mdl.nodes    = vtx;
0066 mdl.elems    = simp;
0067 mdl.boundary = srf;
0068 mdl.boundary_numbers=fc;
0069 mdl.gnd_node=    find_centre_node(vtx);
0070 mdl.np_fwd_solve.perm_sym =     '{n}';
0071 mdl.name = name;
0072 
0073 % Model Stimulation
0074 if ~isempty(stim_pattern)
0075    mdl.stimulation= stim_pattern;
0076 end
0077 
0078 nelec= size(centres,1);
0079 if nelec>0
0080    % Electrodes
0081    [elec,sels,electrodes] = ng_tank_find_elec(srf,vtx,bc,centres);
0082    if size(elec,1) ~= nelec
0083       error('Failed to find all the electrodes')
0084    end
0085 
0086    % set the z_contact
0087    z_contact= z_contact.*ones(nelec,1);
0088    for i=1:nelec
0089       electrodes(i).z_contact= z_contact(i);
0090    end
0091 
0092    mdl.electrode =     electrodes;
0093 end
0094 
0095 mdl.solve=      'eidors_default';
0096 mdl.jacobian=   'eidors_default';
0097 mdl.system_mat= 'eidors_default';
0098 
0099 fwd_mdl= eidors_obj('fwd_model', mdl);
0100 
0101 % Output mat_idx cell array of indices into each material
0102 % typei. Array order is the order of the specified material
0103 % (netgen 'tlo' statements in the .geo file).
0104 function [mat_idx] = mk_mat_indices( mat_ind);
0105   % find length of mat_indices
0106   % test example: mat_ind=[10 12 14 14 12 12 14 12];
0107 
0108   if isempty(mat_ind)
0109      mat_idx = [];
0110      return
0111   end
0112   mat_indices = unique( mat_ind );
0113   for i= 1:length(mat_indices);
0114      mat_idx{i}= find(mat_ind == mat_indices(i));
0115   end
0116 
0117 function gnd_node=    find_centre_node(vtx);
0118   %distance from zero
0119   d = sum( vtx.^2, 2);
0120   [jnk,gnd_node] = min(d);
0121   gnd_node= gnd_node(1);
0122 
0123 function do_unit_test
0124   % Test we can make one electrode
0125   fmdl = ng_mk_cyl_models(1,[],[]);
0126   unit_test_cmp('No elecs',num_elecs(fmdl),0);
0127   fmdl = ng_mk_cyl_models(1,[0,0.5;10,0.2],[0.1]);
0128   unit_test_cmp('Two elecs',num_elecs(fmdl),2);