mk_geophysics_model

 imdl = mk_geophysics_model(str, ne, [option])

 ne  - number of electrodes, 5 metre spacing (+5,+10,...)
       and 0.1 metre diameter
         OR
       a list of electrode locations in the x-dimension or a 2- or
       3-dimensional array, one electrode per row, missing columns
       will be set to zero
       ne = 16
       ne = [4 6 10 20] % 1d: (x)
       ne = [0 1; 2 1.5; 3 1.2; 7 2.5] % 2d: (x,y)
       ne = [0 0.1 1; 2 -0.1 1.5; 3 -0.15 1.2; 7 0 2.5] % 3d: (x,y,z)
 str - model, x = see hmax_rec
       h2x -   2D half-space, linear CEM array (2d fwd)
       h2p5x - 2.5D half-space, linear CEM array (2d fwd + Fourier y-dimension)
       h3x   - 3D half-space, linear CEM array (3d fwd)
       h22x  - 2D half-space, linear CEM array (2d fwd, 2d rec)
       h32x  - 2D half-space, linear CEM array (3d fwd, 2d rec)
       h33x  - 2D half-space, linear CEM array (3d fwd, 3d rec)
       H2x -   2D half-space, linear PEM array (2d fwd)
       H2p5x - 2.5D half-space, linear PEM array (2d fwd + Fourier y-dimension)
       H3x   - 3D half-space, linear PEM array (3d fwd)
       H22x  - 2D half-space, linear PEM array (2d fwd, 2d rec)
       H32x  - 2D half-space, linear PEM array (3d fwd, 2d rec)
       H33x  - 2D half-space, linear PEM array (3d fwd, 3d rec)
 opt - override default configuration options (optional cell array)
       'hmax_fwd' - fine reconstruction mesh density, given an
                    array width xw, and electrode spacing es
                    Note that for ne=16, 'A' and 'a' are equivalent.
                ['a' : hmax_fwd=xw/1;    ['A' : hmax_fwd=es*16;
                 'b' : hmax_fwd=xw/2;     'B' : hmax_fwd=es*8;
                 'c' : hmax_fwd=xw/4;     'C' : hmax_fwd=es*4;
                 ...                       ...
                 'z' : hmax_fwd=xw/2^25]  'Z' : hmax_fwd=es*2^-21]
       'hmax_fwd_inner'
                  - reconstruction model mesh density for the inner region
                    [default: 1/2 of the outer region density hmax_fwd]
       'hmax_rec' - reconstruction model mesh density [hmax_fwd*2]
       'elec_width' - electrode width [0.1 m]
                    width = 0 requests a PEM, rather than CEM
       'z_contact' - electrode contact impedance [0.01 \Ohm.m]
       'elec_spacing' - distance between electrode centers [5 m]
       'extend_x' - extra mesh in the principle axis of the
                    electrode array, multiple of array width [1]
       'extend_y' - extra mesh in the minor axis of the
                    electrode array, multiple of array width
                    (3D models only) [1]
       'extend_z' - extra depth of model, multiple of array width [1]
       'extend_inner_x'
                  - inner (denser) mesh, multiple of array width [3/5]
       'extend_inner_y'
                  - inner (denser) mesh, multiple of array width [3/5]
                    (3D models only)
       'extend_inner_z'
                  - inner (denser) mesh, multiple of array width [2/5]
       'skip_c2f' - skip building the rec_model to fwd_model mapping [0]
       'threshold' - threshold for electrode placement errors [1e-12]
                    if error exceeds threshold, mash nodes by cubic interp
                    maintaining the side and lower boundaries of the
                    model... the nodes will be mashed until the electrodes
                    conform to the requested electrode positions,
                    correcting for Netgen inaccuracies and lack of
                    vertical profile (Inf = disable node mashing)
       'build_stim' - use stim_pattern_geophys to build a standard geophysics
                      stim/meas sequence and add it to the model, based on the
                      number of electrodes, and assuming a co-linear array;
                      set to 'none' to skip this step
                      [default: 'Wenner']
       'extra_ng_code' - pass extra netgen code, any 'tlo' are subtracted from
                      the inner region, units are scaled, so left to right-most
                      electrodes are -1 to +1
                      [default: {}]

 The linear electrode array runs in the +X direction at Z=0. For
 the 3D model, the Y-axis is perpendicular to the electrode array.

 (C) 2015--2017 A. Boyle
 License: GPL version 2 or version 3

0001 function imdl = mk_geophysics_model(str, ne, opt);
0002 % imdl = mk_geophysics_model(str, ne, [option])
0003 %
0004 % ne  - number of electrodes, 5 metre spacing (+5,+10,...)
0005 %       and 0.1 metre diameter
0006 %         OR
0007 %       a list of electrode locations in the x-dimension or a 2- or
0008 %       3-dimensional array, one electrode per row, missing columns
0009 %       will be set to zero
0010 %       ne = 16
0011 %       ne = [4 6 10 20] % 1d: (x)
0012 %       ne = [0 1; 2 1.5; 3 1.2; 7 2.5] % 2d: (x,y)
0013 %       ne = [0 0.1 1; 2 -0.1 1.5; 3 -0.15 1.2; 7 0 2.5] % 3d: (x,y,z)
0014 % str - model, x = see hmax_rec
0015 %       h2x -   2D half-space, linear CEM array (2d fwd)
0016 %       h2p5x - 2.5D half-space, linear CEM array (2d fwd + Fourier y-dimension)
0017 %       h3x   - 3D half-space, linear CEM array (3d fwd)
0018 %       h22x  - 2D half-space, linear CEM array (2d fwd, 2d rec)
0019 %       h32x  - 2D half-space, linear CEM array (3d fwd, 2d rec)
0020 %       h33x  - 2D half-space, linear CEM array (3d fwd, 3d rec)
0021 %       H2x -   2D half-space, linear PEM array (2d fwd)
0022 %       H2p5x - 2.5D half-space, linear PEM array (2d fwd + Fourier y-dimension)
0023 %       H3x   - 3D half-space, linear PEM array (3d fwd)
0024 %       H22x  - 2D half-space, linear PEM array (2d fwd, 2d rec)
0025 %       H32x  - 2D half-space, linear PEM array (3d fwd, 2d rec)
0026 %       H33x  - 2D half-space, linear PEM array (3d fwd, 3d rec)
0027 % opt - override default configuration options (optional cell array)
0028 %       'hmax_fwd' - fine reconstruction mesh density, given an
0029 %                    array width xw, and electrode spacing es
0030 %                    Note that for ne=16, 'A' and 'a' are equivalent.
0031 %                ['a' : hmax_fwd=xw/1;    ['A' : hmax_fwd=es*16;
0032 %                 'b' : hmax_fwd=xw/2;     'B' : hmax_fwd=es*8;
0033 %                 'c' : hmax_fwd=xw/4;     'C' : hmax_fwd=es*4;
0034 %                 ...                       ...
0035 %                 'z' : hmax_fwd=xw/2^25]  'Z' : hmax_fwd=es*2^-21]
0036 %       'hmax_fwd_inner'
0037 %                  - reconstruction model mesh density for the inner region
0038 %                    [default: 1/2 of the outer region density hmax_fwd]
0039 %       'hmax_rec' - reconstruction model mesh density [hmax_fwd*2]
0040 %       'elec_width' - electrode width [0.1 m]
0041 %                    width = 0 requests a PEM, rather than CEM
0042 %       'z_contact' - electrode contact impedance [0.01 \Ohm.m]
0043 %       'elec_spacing' - distance between electrode centers [5 m]
0044 %       'extend_x' - extra mesh in the principle axis of the
0045 %                    electrode array, multiple of array width [1]
0046 %       'extend_y' - extra mesh in the minor axis of the
0047 %                    electrode array, multiple of array width
0048 %                    (3D models only) [1]
0049 %       'extend_z' - extra depth of model, multiple of array width [1]
0050 %       'extend_inner_x'
0051 %                  - inner (denser) mesh, multiple of array width [3/5]
0052 %       'extend_inner_y'
0053 %                  - inner (denser) mesh, multiple of array width [3/5]
0054 %                    (3D models only)
0055 %       'extend_inner_z'
0056 %                  - inner (denser) mesh, multiple of array width [2/5]
0057 %       'skip_c2f' - skip building the rec_model to fwd_model mapping [0]
0058 %       'threshold' - threshold for electrode placement errors [1e-12]
0059 %                    if error exceeds threshold, mash nodes by cubic interp
0060 %                    maintaining the side and lower boundaries of the
0061 %                    model... the nodes will be mashed until the electrodes
0062 %                    conform to the requested electrode positions,
0063 %                    correcting for Netgen inaccuracies and lack of
0064 %                    vertical profile (Inf = disable node mashing)
0065 %       'build_stim' - use stim_pattern_geophys to build a standard geophysics
0066 %                      stim/meas sequence and add it to the model, based on the
0067 %                      number of electrodes, and assuming a co-linear array;
0068 %                      set to 'none' to skip this step
0069 %                      [default: 'Wenner']
0070 %       'extra_ng_code' - pass extra netgen code, any 'tlo' are subtracted from
0071 %                      the inner region, units are scaled, so left to right-most
0072 %                      electrodes are -1 to +1
0073 %                      [default: {}]
0074 %
0075 % The linear electrode array runs in the +X direction at Z=0. For
0076 % the 3D model, the Y-axis is perpendicular to the electrode array.
0077 %
0078 % (C) 2015--2017 A. Boyle
0079 % License: GPL version 2 or version 3
0080 
0081 % model: 64 electrode, 2d half-space
0082 % Once upon a time, this code started out from the following tutorial.
0083 % model from http://eidors3d.sourceforge.net/tutorial/other_models/square_mesh.shtml
0084 
0085 if ischar(str) && strcmp(str,'UNIT_TEST'); do_unit_test; return; end
0086 copt.fstr = 'mk_geophysics_model';
0087 if nargin < 3
0088    opt = {};
0089 end
0090 SALT='z$Id: mk_geophysics_model.m 5903 2019-02-05 11:08:03Z aadler $'; % stick a key in the model 'save' file, so we can expire them when the model definitions age
0091 imdl = eidors_cache(@mk_model,{str, ne, opt, SALT}, copt);
0092 
0093 function imdl=mk_model(str,ne,opt,SALT);
0094 if str(1) ~= 'h' && str(1) ~= 'H'
0095    error([str ': I only know how to build linear half-space models: h***']);
0096 end
0097 
0098 MDL_2p5D_CONFIG = 0;
0099 switch str(2:end-1)
0100    case {'2', '3'} % simple meshes
0101       FMDL_DIM = str(2) - '0';
0102       CMDL_DIM = 0; % no cmdl
0103    case '2p5' % 2.5D Fourier transformed
0104       FMDL_DIM = 2;
0105       CMDL_DIM = 0; % no cmdl
0106       MDL_2p5D_CONFIG = 1;
0107    case {'22', '33', '32'} % dual meshes
0108       FMDL_DIM = str(2) - '0';
0109       CMDL_DIM = str(3) - '0';
0110    otherwise
0111       error([str ': unrecognized model type']);
0112 end
0113 assert(CMDL_DIM ~= 3, '3d rec_model not yet tested');
0114 
0115 skip_c2f = 0;
0116 if str(1) == 'h'
0117    elec_width = 0.1;
0118 else % str(1) == 'H'
0119    elec_width = 0;
0120 end
0121 z_contact= 0.01;
0122 nodes_per_elec= 3; %floor(elec_width/hmax_rec*10);
0123 elec_spacing= 5.0;
0124 threshold = 1e-12;
0125 save_model_to_disk = (length(ne) == 1) && (length(opt) == 0);
0126 
0127 extend_x = 1;
0128 extend_y = 1;
0129 extend_z = 1;
0130 extend_inner_x = 3/5;
0131 extend_inner_y = 3/5;
0132 extend_inner_z = 2/5;
0133 build_stim = 'Wenner';
0134 extra_ng_code = '';
0135 if length(opt) > 0 % allow overriding the default values
0136    assert(round(length(opt)/2)*2 == length(opt),'option missing value?');
0137    expect = {'hmax_rec','hmax_fwd', 'hmax_fwd_inner', ...
0138              'elec_width','z_contact','elec_spacing',...
0139              'extend_x', 'extend_y', 'extend_z', ...
0140              'extend_inner_x', 'extend_inner_y', 'extend_inner_z', ...
0141              'skip_c2f', 'threshold', 'build_stim', 'extra_ng_code'};
0142    opts = struct(opt{:})
0143    for i = fieldnames(opts)'
0144       assert(any(strcmp(i,expect)), ['unexpected option: ',i{:}]);
0145       eval([i{:} ' = opts.(i{:});']);
0146    end
0147    if (str(1) == 'H') && isfield(opts, 'elec_width')
0148       error('requested "H" PEM model but configured "elec_width" option');
0149    end
0150 end
0151 if length(ne) == 1 % ne: number of electrodes
0152    xw=(ne-1)*elec_spacing; % array width
0153    %xs=-(ne-1)*elec_spacing/2; % array centered
0154    xs=+5; % array at left-most at +5
0155    xyz = xs+([1:ne]'-1)*elec_spacing;
0156 else
0157    xyz = ne; % must be a set of coordinates for the electrodes...
0158 end
0159 if size(xyz,1) == 1
0160    xyz = xyz'; % flip to column
0161 end
0162 xyz = [xyz zeros(size(xyz,1),3-size(xyz,2))]; % [x 0 0] or [ x y 0 ] or [ x y z ]
0163 ne=size(xyz,1);
0164 [R, X] = rot_line_to_xaxis(xyz);
0165 % % rescale, centre electrodes so NetGen can be happy
0166 xyzc = (xyz - X)*R; % centre and scale electrodes: -1 to +1 y-axis
0167 % xyzc = xyz; % centre and scale electrodes: -1 to +1 y-axis
0168 xw=max(xyzc(:,1))-min(xyzc(:,1));
0169 xs=min(xyzc(:,1));
0170 elec_spacing = min(min(pdist(xyzc) + diag(inf*(1:size(xyzc,1))))); % min spacing btw elec
0171 
0172 if ~exist('hmax_fwd','var')
0173    if str(end)-'a' >= 0
0174       hmax_fwd = xw*2^-(str(end)-'a');
0175    else
0176       hmax_fwd = elec_spacing*2^-(str(end)-'A'-4);
0177    end
0178 end
0179 if ~exist('hmax_rec','var') % allow hmax_rec to depend on configured hmax_fwd
0180    hmax_rec=hmax_fwd*2.01; % avoid parametrization aliasing
0181 end
0182 if ~exist('hmax_fwd_inner','var') % allow hmax_rec to depend on configured hmax_fwd
0183    hmax_fwd_inner=hmax_fwd/2.0;
0184 end
0185 
0186 if save_model_to_disk
0187    isOctave = exist('OCTAVE_VERSION');
0188    if isOctave
0189       octavestr='-octave-';
0190    else
0191       octavestr='-';
0192    end
0193    % NOTE models are stored in the directory specified by eidors_cache('cache_path')
0194    filename=fullfile(eidors_cache('cache_path'),sprintf('mk_geophysics_model%simdl-%s-%03del.mat',octavestr,str,ne));
0195    if exist(filename, 'file') == 2
0196       tmp = whos('-file',filename,'SALT');
0197       if length(tmp) > 0
0198          tmp = load(filename,'SALT');
0199          fSALT = tmp.SALT;
0200       else
0201          fSALT = 'deadbeef';
0202       end
0203       if strcmp(fSALT, SALT)
0204          tmp=load(filename,'imdl');
0205          imdl = tmp.imdl;
0206          eidors_msg(sprintf('%s: %s, %d electrode model loaded from file',filename,str,ne));
0207          return
0208       end % hmm, the SALT doesn't match so we go back to generating a new model from scratch
0209    end
0210 end
0211 
0212 assert(extend_x>0,'extend_x must be > 0');
0213 assert(extend_y>0,'extend_y must be > 0');
0214 assert(extend_z>-1,'extend_z must be > -1');
0215 
0216 % Calculate cmdl, fmdl and inner fmdl (fmdlin) min/max coordinates.
0217 % After rescaling (normalizing), electrode major axis is along the x-axis, and
0218 % the electrode array will be 2 units long (-1,+1).
0219 fmdlbox =   [-(1+2*extend_x)        +(1+2*extend_x);
0220              -(1+2*extend_y)        +(1+2*extend_y);
0221              -(2+2*extend_z)        +3];
0222 fmdlinbox = [-(1+2*extend_inner_x)  +(1+2*extend_inner_x);
0223              -(1+2*extend_inner_y)  +(1+2*extend_inner_y);
0224              -(2+2*extend_inner_z)  +2];
0225 cmdlbox = fmdlbox;
0226 
0227 assert(all(fmdlbox(:,1) <= fmdlinbox(:,1)), 'oops, inner mesh must be smaller than outer mesh');
0228 assert(all(fmdlbox(:,2) >= fmdlinbox(:,2)), 'oops, inner mesh must be smaller than outer mesh');
0229 assert(all(fmdlbox(:,1) <= cmdlbox(:,1)), 'oops, reconstruction mesh must be smaller than forward mesh');
0230 assert(all(fmdlbox(:,2) >= cmdlbox(:,2)), 'oops, reconstruction mesh must be smaller than forward mesh');
0231 
0232 % build shape string for NetGen
0233 if elec_width ~= 0
0234    elec_shape = [elec_width*norm(R)/2, 0, elec_width*norm(R)/2/(nodes_per_elec-1)]; % CEM, circular, maxh
0235    elec_pos   = [ xyzc(:,1:FMDL_DIM), repmat([zeros(1,3-FMDL_DIM+2) 1],ne,1) ]; % p(x,y,z=0), n(0,0,1)
0236    cem2pem = 0;
0237 else
0238    if FMDL_DIM == 3
0239       elec_width = elec_spacing/2;
0240       elec_shape = [elec_width, elec_width, hmax_fwd_inner]; % PEM, sz, maxh
0241       elec_pos   = [ xyzc, repmat([0 0 1],ne,1) ]; % p(x,y,z=0), n(0,0,1)
0242       elec_pos(:,1) = elec_pos(:,1) + elec_width/2;
0243       elec_pos(:,2) = elec_pos(:,2) + elec_width/2;
0244       cem2pem = 1;
0245    else % FMDL_DIM == 2
0246       elec_width = elec_spacing/2;
0247       elec_shape = [elec_width, elec_width, hmax_fwd_inner]; % rectangular CEM->PEM, maxh
0248       elec_pos   = [ xyzc(:,1:2), repmat([0 0 0 1],ne,1) ]; % p(x,y,z=0), n(0,0,1)
0249       elec_pos(:,1) = elec_pos(:,1) + elec_width/2;
0250       cem2pem = 1;
0251    end
0252 end
0253 skip_tlo = '';
0254 for idx = strfind(extra_ng_code, 'tlo ')
0255    sc = strfind(extra_ng_code, ';'); % semicolons
0256    sc = min(sc + (sc<idx)*1e10);
0257    tlo = extra_ng_code(idx+4:sc-1);
0258    skip_tlo = [skip_tlo ' and (not ' tlo ')'];
0259 end
0260 elec_pos(:,FMDL_DIM) = 0; % flatten electrode positions onto the "ps" plane
0261 if FMDL_DIM == 3
0262    shape_str = [...
0263                 sprintf('solid ps = plane(%s);\n', a2s([0 0 0; 0 0 1])), ...
0264                 sprintf('solid bi = orthobrick(%s);\n', a2s(fmdlinbox')), ...
0265                 sprintf('solid bo = orthobrick(%s);\n', a2s(fmdlbox')), ...
0266                 extra_ng_code, ...
0267                 sprintf('solid ri = bi and ps%s -maxh=%f;\n', skip_tlo, hmax_fwd_inner), ...
0268                 sprintf('solid ro = bo and ps and (not bi) -maxh=%f;\n', hmax_fwd), ...
0269                 sprintf('solid mainobj = ri;\n')];
0270    % Note that ri must be the 'mainobj' so that it can intersect with the electrodes
0271    % additional top level objects for netgen
0272 else % netgen 2d model
0273    % to create the 2D slice we need to give NetGen something to work with
0274    %Need some width to let netgen work, but not too much so
0275    % that it meshes the entire region
0276    sw = range(fmdlinbox(1,:)) / 5; % initial estimate
0277    sw = min(sw,2*hmax_fwd); % coarse model slice width
0278    ri2d = fmdlinbox'; ri2d(:,2) = [-sw 0 ];
0279    ro2d = fmdlbox'; ro2d(:,2) = [-sw 0 ];
0280    shape_str = [...
0281                 sprintf('solid ps = plane(%s);\n', a2s([0 0 0; 0 0 1])), ...
0282                 sprintf('solid bi = orthobrick(%s);\n', a2s(ri2d)), ...
0283                 sprintf('solid bo = orthobrick(%s);\n', a2s(ro2d)), ...
0284                 extra_ng_code, ...
0285                 sprintf('solid ri = bi and ps%s -maxh=%f;\n', skip_tlo, hmax_fwd_inner), ...
0286                 sprintf('solid ro = bo and ps and (not bi) -maxh=%f;\n', hmax_fwd), ...
0287                 sprintf('solid mainobj = ri;\n')];
0288 end
0289 % fprintf('SHAPE_STR: %s', shape_str); elec_pos
0290 elec_obj = 'ps';
0291 [fmdl, mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj, 'tlo ro');
0292 if FMDL_DIM == 2 % 2D
0293    % now convert the roughly 2D slice into a true 2D plane
0294    [fmdl, mat_idx] = copy_mdl2d_from3d(fmdl, mat_idx, 'y');
0295 else % 3D
0296    if CMDL_DIM ~= 0
0297       % c2f
0298       cmdl.mk_coarse_fine_mapping.f2c_offset  = [0 0 0];
0299       cmdl.mk_coarse_fine_mapping.f2c_project = [1 0 0; 0 0 1; 0 1 0];
0300       % duplicate parameters since mk_analytic/approx_c2f have different names...
0301       cmdl.mk_analytic_c2f.f2c_offset  = cmdl.mk_coarse_fine_mapping.f2c_offset;
0302       cmdl.mk_analytic_c2f.f2c_project = cmdl.mk_coarse_fine_mapping.f2c_project;
0303    end
0304 end
0305 
0306 if cem2pem
0307    fmdl = convert_cem2pem(fmdl, xyzc);
0308 end
0309 
0310 % 2d cmdl
0311 xllim = fmdlbox(1,1);
0312 xrlim = fmdlbox(1,2);
0313 zdepth = fmdlbox(3,1);
0314 xr=max(floor((xrlim-xllim)/hmax_rec/2),1)*2+1; % odd number
0315 yr=max(floor(-zdepth/hmax_rec/2),1)*2+1; % odd number
0316 [x,y] = meshgrid( linspace(xllim,xrlim,xr), linspace(zdepth,0,yr));
0317 vtx= [x(:),y(:)];
0318 if CMDL_DIM ~= 0
0319    cmdl= mk_fmdl_from_nodes( vtx,{vtx(1,:)}, z_contact, 'sq_m2');
0320 end
0321 
0322 % stick electrode nodes into cmdl so that show_fem will plot them
0323 for i=1:ne
0324    n=fmdl.electrode(i).nodes;
0325    nn=length(n);
0326    nx=fmdl.nodes(n,:);
0327 
0328    fmdl.electrode(i).z_contact = z_contact;
0329    if CMDL_DIM ~= 0
0330       nnc = length(cmdl.nodes);
0331       cmdl.nodes = [cmdl.nodes; nx(:,[1 FMDL_DIM])];
0332       cmdl.electrode(i).nodes = (nnc+1):(nnc+nn);
0333       cmdl.electrode(i).z_contact = z_contact;
0334    end
0335 end
0336 
0337 % fix electrode locations if necessary
0338 elec_err = sqrt(sum(mdl_elec_err(fmdl, xyzc).^2,2));
0339 if max(elec_err) > threshold % put electrodes in the right place
0340    [fmdl, cf] = correct_electrode_positions(fmdl, xyzc);
0341 
0342    if CMDL_DIM ~= 0
0343       [cmdl, cc] = correct_electrode_positions(cmdl, xyzc);
0344    end
0345 end
0346 % save functions for later use
0347 fmdl.mv_elec = @shift_electrode_positions;
0348 
0349 % reverse the centre and scaling
0350 nn = size(fmdl.nodes,1);
0351 Xn = repmat(X(1,:), nn, 1);
0352 fmdl.nodes = ([fmdl.nodes zeros(nn,3-FMDL_DIM)]/ R) + Xn;
0353 fmdl.nodes = fmdl.nodes(:,1:FMDL_DIM);
0354 
0355 if CMDL_DIM ~= 0
0356    nn = size(cmdl.nodes,1);
0357    Xn = repmat(X(1,:), nn, 1);
0358    if CMDL_DIM ~= FMDL_DIM % 2D
0359       cmdl.nodes = ([cmdl.nodes(:,1) zeros(nn,1) cmdl.nodes(:,2:end)]/ R) + Xn;
0360       cmdl.nodes = cmdl.nodes(:,[1 3])/R([1 3],[1 3]);
0361    else % CMDL_DIM == FMDL_DIM
0362       cmdl.nodes = ([cmdl.nodes zeros(nn,3-CMDL_DIM)]/ R) + Xn;
0363       cmdl.nodes = cmdl.nodes(:,1:CMDL_DIM);
0364    end
0365 end
0366 
0367 % check that all electrodes were found
0368 for i = 1:length(fmdl.electrode)
0369    nn = fmdl.electrode(i).nodes;
0370    assert(length(nn(:)) > 0, sprintf('electrode#%d: failed to find nodes',i));
0371 end
0372 
0373 % show_fem(fmdl); xlabel('x'); ylabel('y'); zlabel('z');
0374 
0375 % TODO ... Actually we want a node in the middle: the boundary is bad too...
0376 [~, gn] = min(fmdl.nodes(:,end));
0377 fmdl.gnd_node = gn; % make sure the ground node is away from surface electrodes
0378 if CMDL_DIM ~= 0
0379    [~, gn] = min(cmdl.nodes(:,end));
0380    cmdl.gnd_node = gn; % make sure the ground node is away from surface electrodes
0381 end
0382 
0383 
0384 imdl= mk_common_model('a2d0c',4); % 2d model
0385 imdl.fwd_model = fmdl;
0386 imdl.name = ['EIDORS mk_geophysics_model ' str];
0387 imdl.fwd_model.name = ['EIDORS mk_geophysics_model fwd model ' str];
0388 if CMDL_DIM ~= 0
0389    imdl.rec_model.name = ['EIDORS mk_geophysics_model rec model ' str];
0390    imdl.rec_model = cmdl;
0391    % EIDORS "analytic_c2f" gets stuck, do an approximate one
0392    %eidors_default('set','mk_coarse_fine_mapping','mk_analytic_c2f');
0393    %eidors_default('set','mk_coarse_fine_mapping','mk_approx_c2f');
0394    %[c2f,out] = mk_coarse_fine_mapping(fmdl,cmdl);
0395    if ~skip_c2f
0396       [c2f,out] = mk_approx_c2f(fmdl,cmdl);
0397       imdl.fwd_model.coarse2fine = c2f;
0398       imdl.fwd_model.background = out;
0399    end
0400 end
0401 
0402 if MDL_2p5D_CONFIG
0403    imdl.fwd_model.jacobian   = @jacobian_adjoint_2p5d_1st_order;
0404    imdl.fwd_model.solve      = @fwd_solve_2p5d_1st_order;
0405    imdl.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0406    imdl.fwd_model.jacobian_adjoint_2p5d_1st_order.k = [0 3];
0407    imdl.fwd_model.fwd_solve_2p5d_1st_order.k = [0 3];
0408 end
0409 
0410 if ~strcmp(build_stim,'none')
0411    imdl.fwd_model.stimulation = stim_pattern_geophys(length(imdl.fwd_model.electrode), build_stim);
0412 end
0413 
0414 if save_model_to_disk
0415    save(filename,'imdl','SALT');
0416 end
0417 
0418 % convert 2x3 array to "x1,y1,z1;x2,y2,z2" string
0419 % convert 1x3 array to "x1,y1,z1" string
0420 % if 1x1 array, then use xy_ctr=x1,y1 and [0 0 1]=x2,y2,z2 (z+)
0421 function s = a2s(a)
0422 if length(a(:)) == 3
0423    s = sprintf('%f,%f,%f', ...
0424                a(1), a(2), a(3));
0425 else
0426    s = sprintf('%f,%f,%f;%f,%f,%f', ...
0427                 a(1,1), a(1,2), a(1,3), ...
0428                a(2,1), a(2,2), a(2,3));
0429 end
0430 
0431 % returns R rotation/scaling and X0 offset
0432 % xyz1 = R * xyz + X; % rotate and scale to +/- 1
0433 function [R,X,var] = rot_line_to_xaxis(xyz)
0434 x = xyz(:,1); y=xyz(:,2); z=xyz(:,3);
0435 
0436 % fit line to points
0437 % http://www.mathworks.com/matlabcentral/newsreader/view_thread/32502
0438 p = mean(xyz);
0439 [U,S,V] = svd([x-p(1), y-p(2), z-p(3)]);
0440 if V(end,1) ~= 0
0441    N=1/V(end,1)*V(:,1);
0442 else
0443    N=V(:,1);
0444 end
0445 A=p' + dot( xyz(1,  :) - p, N ) * N/norm(N)^2;
0446 B=p' + dot( xyz(end,:) - p, N ) * N/norm(N)^2;
0447 
0448 % rotate line to +y-axis
0449 a = N/norm(N);
0450 b = [ 1 0 0 ]'; % +x
0451 % http://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
0452 v = cross(a, b);
0453 s = norm(v); c = dot(a, b); % sin, cos
0454 xt = @(v) [   0  -v(3)  v(2); ... % skew symmetric cross-product of v
0455             v(3)    0  -v(1); ... % diagnoal is the scaling identity matrix
0456            -v(2)  v(1)    0];
0457 if abs(s) < eps*1e3
0458    R = eye(3);
0459 else
0460    R = (eye(3) + xt(v) + xt(v)^2 * (1-c)/s^2);
0461    R=R'; % for right multiply: xyz*R
0462 end
0463 
0464 X = repmat(p, size(xyz,1), 1);
0465 R = R / max(range((xyz-X)*R)) * 2;
0466 
0467 DEBUG=0;
0468 if DEBUG
0469    clf;
0470    subplot(211);
0471    plot3(x,y,z,'bo');
0472    hold on;
0473    plot3(x(1),y(1),z(1),'go');
0474    plot3(p(1),p(2),p(3),'ro');
0475    plot3([A(1),B(1)],[A(2),B(2)],[A(3),B(3)]);
0476    grid; axis equal; hold off;
0477    xlabel('x'); ylabel('y'); zlabel('z');
0478    title('pre-rotation');
0479    subplot(212);
0480    xx = (xyz-X)*R;
0481    plot3(xx(:,1), xx(:,2), xx(:,3),'bo');
0482    hold on;
0483    plot3(xx(1,1), xx(1,2), xx(1,3),'go');
0484    plot3(0,0,0,'ro');
0485    grid; axis equal; hold off;
0486    xlabel('x'); ylabel('y'); zlabel('z');
0487    title('post-rotation');
0488 end
0489 
0490 % the (max-min) range of a variable's values
0491 function r = range(a)
0492 r = max(a(:))-min(a(:));
0493 
0494 function [mdl2,idx2] = copy_mdl2d_from3d(mdl3,idx3,sel);
0495 % AB: taken from EIDORS function ng_mk_gen_models() subfunction of the same name
0496 % AB: NEW: sel = 'x', 'y' or 'z' -- default was Z, we want X
0497    if sel == 'x'
0498       % swap Z and X
0499       T = [ 0 0 1; 0 1 0; 1 0 0 ];
0500    elseif sel == 'y'
0501       % swap Z and Y
0502       T = [ 1 0 0; 0 0 1; 0 1 0 ];
0503    elseif sel == 'z'
0504       T = eye(3);
0505    else
0506       error('sel must be "x", "y" or "z"');
0507    end
0508    mdl3.nodes = mdl3.nodes * T; % AB: SWAP axes
0509 
0510    % set name
0511    mdl2 = eidors_obj('fwd_model',sprintf('%s 2D',mdl3.name));
0512 
0513    % set nodes
0514    [bdy,idx] = find_boundary(mdl3.elems);
0515    vtx = mdl3.nodes;
0516    z_vtx = reshape(vtx(bdy,3), size(bdy) );
0517    z_vtx_thres = max(z_vtx(:))-10*eps*range(z_vtx(:));
0518    lay0  = find( all(z_vtx >= z_vtx_thres, 2) );
0519    bdy0  = bdy( lay0, :);
0520 
0521    vtx0  = unique(bdy0(:));
0522    mdl2.nodes = vtx(vtx0,1:2);
0523 
0524    % set elems
0525    nmap  = zeros(size(vtx,1),1); nmap(vtx0) = 1:length(vtx0);
0526    bdy0  = reshape(nmap(bdy0), size(bdy0) ); % renumber to new scheme
0527    mdl2.elems = bdy0;
0528 
0529    % set boundary
0530    mdl2.boundary = find_boundary( mdl2.elems);
0531 
0532    % set gnd_node
0533    mdl2.gnd_node = nmap(mdl3.gnd_node);
0534 
0535    % set material indices
0536    % TODO: vectorize code
0537    idx2 = {};
0538    idx0  = idx( lay0, :);
0539    for i=1:size(idx3,2)
0540       idx2{i} = [];
0541       ii = 1;
0542       for j=1:size(idx3{i},1)
0543          idx_tmp = find( idx0==idx3{i}(j) );
0544          if not(isempty(idx_tmp))
0545             idx2{i}(ii,1) = idx_tmp(1,1);
0546             ii = ii + 1;
0547          end
0548       end
0549    end
0550 
0551    % set electrode
0552    if isfield(mdl3,'electrode')
0553       mdl2.electrode = mdl3.electrode;
0554       for i=1:length(mdl2.electrode);
0555          nn = mdl3.nodes(mdl3.electrode(i).nodes,:);
0556          enodes = nmap( mdl2.electrode(i).nodes );
0557          enodes(enodes==0) = []; % Remove 3D layers
0558          mdl2.electrode(i).nodes = enodes(:)';
0559       end
0560    end
0561 
0562    ignore = {'electrode', 'nodes', 'boundary', 'elems', 'gnd_node', 'boundary_numbers', 'mat_idx'};
0563    for n=fieldnames(mdl3)'
0564       if ~any(strcmp(n,ignore))
0565          mdl2.(n{:}) = mdl3.(n{:});
0566       end
0567    end
0568 
0569 function mdl = convert_cem2pem(mdl, xyzc)
0570    if ~isfield(mdl, 'electrode')
0571       return;
0572    end
0573    nd = size(mdl.nodes,2); % number of dimensions
0574    for i=1:length(mdl.electrode)
0575       en = mdl.electrode(i).nodes;
0576       nn = mdl.nodes(en,:); % all nodes for this electrode
0577       if nd == 2
0578          np = xyzc(i,[1 3]); % true elec location (2d)
0579       else
0580          np = xyzc(i,:); % true elec location (3d)
0581       end
0582       D = pdist([np; nn]);
0583       [~,idx] = min(D(1,2:end)); % closest CEM node to ideal PEM location
0584       mdl.electrode(i).nodes = en(idx);
0585       % NOTE: used to bump node's location but now we place a corner of the
0586       % square electrode at precisely the correct location
0587       %mdl.nodes(en,:) = np; % shift PEM node to true location
0588    end
0589 
0590 function D2 = pdist(X) % row vectors
0591    if nargin == 2; error('only supports Euclidean distances'); end
0592    %D2 = bsxfun(@plus, dot(X, X, 1)', dot(Y, Y, 1)) - 2*(X'*Y) % 1d
0593    D2 = bsxfun(@minus, X, permute(X,[3 2 1]));
0594    D2 = squeeze(sqrt(sum(D2.^2,2)));
0595 
0596 function [mdl, c] = shift_electrode_positions(mdl, dx)
0597    for i = 1:length(mdl.electrode)
0598       en = mdl.electrode(i).nodes;
0599       ex = mdl.nodes(en,:);
0600       ep(i,:) = (max(ex,[],1) + min(ex,[],1))/2; % mid-point
0601    end
0602    [mdl, c] = correct_electrode_positions(mdl, ep + dx);
0603 
0604 function [mdl, c] = correct_electrode_positions(mdl, xyzc)
0605    %eidors_msg('correct_electrode_positions');
0606    nd = size(mdl.nodes,2);
0607    c = 0; err = 1;
0608    while max(err) > eps
0609       for n = 1:nd
0610          switch(n)
0611             case 1
0612                mdl = mash_nodes(mdl, 'shift_all',     1, 1, xyzc); % X (downslope)
0613             case 2 % 2d: Y, 3d: Z
0614                mdl = mash_nodes(mdl, 'shift_surface', 1, nd, xyzc); % Y or Z (vertical)
0615             case 3 % note: we only do this for 3d
0616                mdl = mash_nodes(mdl, 'shift_middle',  1, 2, xyzc); % Y (cross-slope)
0617             otherwise
0618                error('duh!');
0619          end
0620       end
0621       err = sqrt(sum(mdl_elec_err(mdl, xyzc).^2,2));
0622       c=c+1;
0623       if c >= 100
0624          break;
0625       end
0626    end
0627 
0628 function   mdl = mash_nodes(mdl, method, idm, dim, elec_true)
0629    elec_err = mdl_elec_err(mdl, elec_true);
0630    err = elec_err(:,dim);
0631 
0632    % add borders for electrode positions at
0633    % the volume boundary and 50% of the edge to electrode distance
0634    xq = mdl.nodes(:,idm);
0635    x = [min(xq); ...
0636         mean([min(xq) min(elec_true(:,idm))]); ...
0637         elec_true(:,idm);
0638         mean([max(xq) max(elec_true(:,idm))]); ...
0639         max(xq)];
0640    v = [0; 0; err; 0; 0];
0641 
0642    % scale error to match the electrode locations
0643    interp_method = 'linear';
0644    if strcmp(method, 'shift_surface')
0645          interp_method = 'pchip';
0646    end
0647    vq = interp1(x, v, xq, interp_method, 'extrap');
0648 
0649    switch method
0650       case 'shift_all'
0651          vqs = 1;
0652       case 'shift_middle'
0653          yq = mdl.nodes(:,dim);
0654          yqr = max(yq)-min(yq); % range
0655          yqm = (max(yq) + min(yq))/2; % middle = (max + min)/2
0656          vqs = 1-abs(yq-yqm)./(yqr/2); % scale the shift depending on x's distance from midline
0657       case 'shift_surface' % assumes positive surface
0658          yq = mdl.nodes(:,dim);
0659          yqr = max(yq) - min(yq);
0660          yqm = min(yq); % min
0661          vqs = abs(yq - yqm)./yqr; % scale by distance from surface
0662       otherwise
0663          error(['unrecognized method: ',method]);
0664    end
0665    mdl.nodes(:,dim) = mdl.nodes(:,dim) + (vq .* vqs);
0666 
0667 % calculate the error in electrode position for the fwd_model
0668 function err = mdl_elec_err(mdl, xyzc)
0669    if ~isfield(mdl, 'electrode')
0670       error('electrodes not available on this model, must supply positional error');
0671    end
0672 
0673    nel=length(mdl.electrode); % number of electrodes
0674    nd=size(mdl.nodes,2); % number of dimensions
0675 
0676    eu = ones(nel,nd)*NaN; % init
0677    for i=1:length(mdl.electrode)
0678       nn = mdl.nodes(mdl.electrode(i).nodes,:); % nodes per electrode
0679       eu(i,:) = (max(nn,[],1) + min(nn,[],1))./2; % approx centre of each electrode
0680    end
0681    err = xyzc(:,1:nd) - eu;
0682 
0683 function do_unit_test
0684    ne = 16;
0685    imdl = mk_geophysics_model('h2p5a', ne);
0686    imdl.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0687    img = mk_image(imdl.fwd_model, 1);
0688    img.fwd_solve.get_all_meas = 1;
0689    vh = fwd_solve_halfspace(img);
0690    vd = fwd_solve(img);
0691 clf; h=plot([vh.meas vd.meas],'o--'); legend('analytic','FEM'); set(gca,'box','off'); set(h,'LineWidth',2);
0692 
0693    imdl1 = mk_geophysics_model('h2a',[1:6]);
0694    imdl2 = mk_geophysics_model('h2a',[1:6]');
0695    imdl3 = mk_geophysics_model('h2a',[1:6]'*[1 0]);
0696    imdl3Hnm2d = mk_geophysics_model('H2a',[1:6],{'threshold',Inf}); % try without mashing nodes, no veritcal geometry... electrodes should be precisely located if the electrodes were correctly placed
0697    imdl3Hnm3d = mk_geophysics_model('H3a',[1:6],{'threshold',Inf}); % try without mashing nodes, no veritcal geometry... electrodes should be precisely located if the electrodes were correctly placed
0698    imdl4 = mk_geophysics_model('h2a',[1:6]'*[1 0] + ([1:6]*0+2)'*[0 1]);
0699    R = @(x) [cosd(x) -sind(x); sind(x) cosd(x)]; % rotation matrix
0700    X = [0 2];
0701    imdl5 = mk_geophysics_model('h2a',([1:6]'*[1 0] + ([1:6]*0+1)'*X)*R(-135));
0702    elec_pos_2d = [1 1; 2 2; 3 1; 4 1.5];
0703    elec_pos_3d = [1 0 0; 2 0.5 1; 3 -0.5 2.5; 10 0 3];
0704    imdl2dc = mk_geophysics_model('h2a', elec_pos_2d); % 2D complete electrode model
0705    imdl3dc = mk_geophysics_model('h3a', elec_pos_3d); % 3D complete electrode model
0706    imdl2dp = mk_geophysics_model('H2a', elec_pos_2d); % 2D point electrode model
0707    imdl3dp = mk_geophysics_model('H3a', elec_pos_3d); % 3D point electrode model
0708    imdl2df = mk_geophysics_model('H2a', elec_pos_2d, {'threshold', Inf}); % 2D point electrode model... without mashing
0709    imdl3df = mk_geophysics_model('H3a', elec_pos_3d, {'threshold', Inf}); % 3D point electrode model... without mashing
0710 
0711    % dual meshes
0712    imdlh32 = mk_geophysics_model('h32a', elec_pos_3d);
0713    imdlh22 = mk_geophysics_model('h22a', elec_pos_2d);
0714    imdlH32 = mk_geophysics_model('H32a', elec_pos_3d);
0715    imdlH22 = mk_geophysics_model('H22a', elec_pos_2d);
0716 
0717    % std dual meshes w/ 16 elec
0718    imdlh32_16 = mk_geophysics_model('h32a', 16);
0719    imdlh22_16 = mk_geophysics_model('h22a', 16);
0720    imdlH32_16 = mk_geophysics_model('H32a', 16);
0721    imdlH22_16 = mk_geophysics_model('H22a', 16);
0722 
0723 if 0
0724    % Nolwenn's grid
0725    [yy,xx] = meshgrid(0:3:24, [0 4:2:20*2 20*2+4]); % 9 x 21 grid
0726    xyz = [xx(:) yy(:) zeros(length(yy(:)),1)];
0727 
0728    % test grid
0729    % 1. bad electrodes ... ? fixme?
0730    % 2. mash_nodes breaks for co-located nodes in any particular dimension... need 2d interp?
0731    [yy,xx]=meshgrid(1:3,1:3); xyz = [xx(:) yy(:) zeros(length(yy(:)),1)]; % 3 x 3 grid
0732    imdl = mk_geophysics_model('H3a',xyz);
0733    % TODO add 'g3a' and 'G3a' types?
0734 end
0735 
0736    unit_test_cmp('h2a halfspace vs default TEST', norm(vh.meas - vd.meas), 0, 4e-3);
0737    unit_test_cmp('1d elec list equivalence (row/col)',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl2));
0738    unit_test_cmp('1d vs. 2d elec list equivalence',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl3));
0739    unit_test_cmp('2D PEM w/o node mashing, no vertical relief',[1:6]'*[1 0], unit_test_elec_pos(imdl3Hnm2d), eps);
0740    unit_test_cmp('2D PEM *is* PEM',length(imdl3Hnm2d.fwd_model.electrode(1).nodes),1)
0741    unit_test_cmp('3D PEM w/o node mashing, no vertical relief',[1:6]'*[1 0 0], unit_test_elec_pos(imdl3Hnm3d), eps);
0742    unit_test_cmp('3D PEM *is* PEM',length(imdl3Hnm3d.fwd_model.electrode(1).nodes),1)
0743    unit_test_cmp('1d vs. 2d + y=2 elec list equivalence',unit_test_elec_pos(imdl1), unit_test_elec_pos(imdl4)-([1:6]*0+2)'*[0 1]);
0744    unit_test_cmp('1d vs. 2d + y=2 - 135 deg elec eq', ...
0745                  unit_test_elec_pos(imdl1), ...
0746                  unit_test_elec_pos(imdl5, R(135), -X), eps*10);
0747    unit_test_cmp('2d with vertical geometry (mash nodes) CEM', ...
0748                  elec_pos_2d, ...
0749                  unit_test_elec_pos(imdl2dc), 0.01);
0750    unit_test_cmp('3d with vertical geometry (mash nodes) CEM', ...
0751                  elec_pos_3d, ...
0752                  unit_test_elec_pos(imdl3dc), 0.001);
0753    unit_test_cmp('2d with vertical geometry (mash nodes) PEM', ...
0754                  elec_pos_2d, ...
0755                  unit_test_elec_pos(imdl2dp), 10*eps);
0756    unit_test_cmp('3d with vertical geometry (mash nodes) PEM', ...
0757                  elec_pos_3d, ...
0758                  unit_test_elec_pos(imdl3dp), 10*eps);
0759    unit_test_cmp('2d with vertical geometry (w/o mash nodes) PEM', ...
0760                  elec_pos_2d, ...
0761                  unit_test_elec_pos(imdl2df), -10*eps);
0762    unit_test_cmp('3d with vertical geometry (w/o mash nodes) PEM', ...
0763                  elec_pos_3d, ...
0764                  unit_test_elec_pos(imdl3df), -10*eps);
0765 
0766    unit_test_cmp('h32a dual 3d/2d', ...
0767                  elec_pos_3d, ...
0768                  unit_test_elec_pos(imdlh32), 0.001);
0769    unit_test_cmp('H32a dual 3d/2d', ...
0770                  elec_pos_3d, ...
0771                  unit_test_elec_pos(imdlH32), 10*eps);
0772    unit_test_cmp('h22a dual 2d/2d', ...
0773                  elec_pos_2d, ...
0774                  unit_test_elec_pos(imdlh22), 0.002);
0775    unit_test_cmp('H22a dual 2d/2d', ...
0776                  elec_pos_2d, ...
0777                  unit_test_elec_pos(imdlH22), 10*eps);
0778 
0779 clf; subplot(221); show_fem(imdl1.fwd_model); title('models match? A');
0780      subplot(222); show_fem(imdl5.fwd_model); title('models match? C');
0781      subplot(223); show_fem(imdl2dc.fwd_model); title('2d deformations');
0782      subplot(224); show_fem(imdl3dc.fwd_model); title('3d deformations'); view([0 -1 0.01]);
0783 
0784 if 1
0785    clf; subplot(221); show_fem(imdlh22.fwd_model);
0786         subplot(222); show_fem(imdlh32.fwd_model);
0787         subplot(223); show_fem(imdlH22.fwd_model);
0788         subplot(224); show_fem(imdlH32.fwd_model);
0789 end
0790 
0791 if 0
0792    clf; subplot(221); show_fem(imdlh22.rec_model);
0793         subplot(222); show_fem(imdlh32.rec_model);
0794         subplot(223); show_fem(imdlH22.rec_model);
0795         subplot(224); show_fem(imdlH32.rec_model);
0796 end
0797 
0798    imdl = mk_geophysics_model('h2p5a', ne, {'extra_ng_code', 'solid tt = orthobrick(-1,-1,-1;-0.5,0,-0.5);\ntlo tt;\n'});
0799    clf; show_fem(imdl.fwd_model);
0800    % Note: this just has to be swallowed safely...
0801 
0802 function xyz = unit_test_elec_pos(imdl, R, X)
0803    if nargin < 2; R = 1; end
0804    if nargin < 3; X = 0; end
0805    fmdl = imdl.fwd_model;
0806    xyz = zeros(length(fmdl.electrode),size(fmdl.nodes,2))*NaN;
0807    for i = 1:length(fmdl.electrode)
0808       nn = fmdl.nodes(fmdl.electrode(i).nodes,:);
0809       xyz(i,:) = (max(nn,[],1) + min(nn,[],1))/2;
0810       xyz(i,:) = xyz(i,:)*R + X;
0811    end