ng_mk_gen_models

 NG_MK_GEN_MODELS: create generic models using netgen
[fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
[fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code);

 INPUT:
 shape_str = string of netgen *.geo file code, with mainobj as the main object

 ELECTRODE POSITIONS:
  elec_pos = [x,y,z,nx,ny,nz] centres of each electrode (N_elecs x 6)
   [x,y,z] is the position, and [nx,ny,nz] is the surface normal

 ELECTRODE SHAPES::
  elec_shape = [width,height, maxsz]  % Rectangular elecs
     OR
  elec_shape = [radius, 0, maxsz ]    % Circular elecs
     OR 
  elec_shape = [0, sz, maxsz ]         % Point electrodes
    (point elecs does some tricks with netgen, using sz square, so the elecs aren't exactly where you ask)

 Specify either a common electrode shape or for each electrode

 ELECTRODE DEFITIONS:
  elec_obj = 'obj_name' or {'name','for','each','elec'} where
    the name is the primitive netgen object (cylinder, plane, etc) which the electrode intersects

 OUTPUT:
  fmdl    - fwd_model object
  mat_idx - indices of materials

 USAGE EXAMPLES:
shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
             'solid bottom = plane(0,0,0;0,0,-1);\n' ...
             'solid top    = plane(0,0,2;0,0,1);\n' ...
             'solid mainobj= top and bottom and cyl -maxh=0.3;\n'];
elec_pos = [  1,  0,  1,   1,  0,  0;
              0,  1,1.2,   0,  1,  0;
              0.8,  0,  0, 0,  0, -1]; 
elec_shape=[0.1];
elec_obj = {'cyl','cyl','bottom'};
fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);

0001 function [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos,  elec_shape, elec_obj, extra_ng_code, mszpoints);
0002 % NG_MK_GEN_MODELS: create generic models using netgen
0003 %[fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0004 %[fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code);
0005 %
0006 % INPUT:
0007 % shape_str = string of netgen *.geo file code, with mainobj as the main object
0008 %
0009 % ELECTRODE POSITIONS:
0010 %  elec_pos = [x,y,z,nx,ny,nz] centres of each electrode (N_elecs x 6)
0011 %   [x,y,z] is the position, and [nx,ny,nz] is the surface normal
0012 %
0013 % ELECTRODE SHAPES::
0014 %  elec_shape = [width,height, maxsz]  % Rectangular elecs
0015 %     OR
0016 %  elec_shape = [radius, 0, maxsz ]    % Circular elecs
0017 %     OR
0018 %  elec_shape = [0, sz, maxsz ]         % Point electrodes
0019 %    (point elecs does some tricks with netgen, using sz square, so the elecs aren't exactly where you ask)
0020 %
0021 % Specify either a common electrode shape or for each electrode
0022 %
0023 % ELECTRODE DEFITIONS:
0024 %  elec_obj = 'obj_name' or {'name','for','each','elec'} where
0025 %    the name is the primitive netgen object (cylinder, plane, etc) which the electrode intersects
0026 %
0027 % OUTPUT:
0028 %  fmdl    - fwd_model object
0029 %  mat_idx - indices of materials
0030 %
0031 % USAGE EXAMPLES:
0032 %shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0033 %             'solid bottom = plane(0,0,0;0,0,-1);\n' ...
0034 %             'solid top    = plane(0,0,2;0,0,1);\n' ...
0035 %             'solid mainobj= top and bottom and cyl -maxh=0.3;\n'];
0036 %elec_pos = [  1,  0,  1,   1,  0,  0;
0037 %              0,  1,1.2,   0,  1,  0;
0038 %              0.8,  0,  0, 0,  0, -1];
0039 %elec_shape=[0.1];
0040 %elec_obj = {'cyl','cyl','bottom'};
0041 %fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0042 
0043 % (C) Andy Adler, 2010. (C) Alistair Boyle, 2013. Licenced under GPL v2 or v3
0044 % $Id: ng_mk_gen_models.m 6100 2021-03-14 19:10:43Z aadler $
0045 
0046 if ischar(shape_str) && strcmp(shape_str,'UNIT_TEST'); do_unit_test; return; end
0047 
0048 if nargin <= 4; extra_ng_code = ''; end
0049 if nargin <= 5; mszpoints = []; else
0050    warning('Use ng_write_opt for mszpoints.');
0051 end
0052 args = { shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code, mszpoints };
0053 copt.cache_on_ng_opt = true;
0054 % Adding ng.opt makes it cache the mtime of the file
0055 copt.fstr = 'ng_mk_gen_models';
0056 
0057 fmdl = eidors_cache(@mk_gen_model,args, copt);
0058 
0059 mat_idx = fmdl.mat_idx;
0060 
0061 function fmdl = mk_gen_model( shape_str, elec_pos, elec_shape, ...
0062                          elec_obj, extra_ng_code, mszpoints);
0063 
0064    fnstem = tempname;
0065    geofn= [fnstem,'.geo'];
0066    ptsfn= [fnstem,'.msz'];
0067    meshfn= [fnstem,'.vol'];
0068 
0069    is2D = 0;
0070    [elecs, centres] = parse_elecs( elec_pos, elec_shape, is2D, elec_obj );
0071 
0072    n_pts = write_geo_file(geofn, ptsfn, shape_str, elecs, extra_ng_code, mszpoints);
0073    if n_pts == 0 
0074       call_netgen( geofn, meshfn);
0075    else
0076       call_netgen( geofn, meshfn, ptsfn);
0077    end
0078 
0079    fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []);
0080 
0081    delete(geofn); delete(meshfn); delete(ptsfn); % remove temp files
0082    if is2D
0083       fmdl = mdl2d_from3d(fmdl);
0084    end
0085 
0086    % convert CEM to PEM if so configured
0087    % TODO shunt model is unsupported
0088    if isfield(fmdl,'electrode');
0089      fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes);
0090    end
0091 
0092 % for the newest netgen, we can't call msz file unless there are actually points in  it
0093 function n_pts_elecs = write_geo_file(geofn, ptsfn, shape_str, ...
0094                           elecs, extra_ng_code, mszpoints);
0095    fid=fopen(geofn,'w');
0096    write_header(fid, shape_str);
0097 
0098    n_elecs = length(elecs);
0099    %  elecs(i).pos   = [x,y,z]
0100    %  elecs(i).shape = 'C' or 'R' or 'P'
0101    %  elecs(i).dims  = [radius] or [width,height]
0102    %  elecs(i).maxh  = '-maxh=#' or '';
0103    pts_elecs_idx = []; 
0104 
0105    % elec_depth is the depth of the netgen object for the electrode
0106    if n_elecs > 1
0107       elec_depth = min(nonzeros(distmat(vertcat( elecs(:).pos ))))/2;
0108    elseif n_elecs==1
0109       % tank_radius = norm( std( vertcat( elecs(:).pos ), 1), 2);
0110       % NOTE: all functions but the point electrode used to use
0111       % tank_radius/4. Point electrode used just tank_radius.
0112       elec_depth = norm( std( vertcat( elecs(:).pos ), 1), 2)/4;
0113    end
0114    for i=1:n_elecs
0115       name = sprintf('elec%04d',i);
0116       pos = elecs(i).pos;
0117       dirn= elecs(i).dirn;
0118       switch elecs(i).shape
0119        case 'C'
0120          write_circ_elec(fid,name, pos, dirn,  ...
0121                elecs(i).dims, elec_depth, elecs(i).maxh);
0122        case 'R'
0123          write_rect_elec(fid,name, pos, dirn,  ...
0124                elecs(i).dims, elec_depth, elecs(i).maxh);
0125        case 'P'
0126          if 0 % Netgen doesn't put elecs where you ask
0127             pts_elecs_idx = [ pts_elecs_idx, i]; 
0128             continue; % DON'T print solid cyl
0129          end
0130          write_rect_elec(fid,name, pos, dirn,  ...
0131                elecs(i).dims, elec_depth, elecs(i).maxh);
0132 
0133        case 'U'
0134          eidors_msg('user defined electrode %d',i, 4);
0135          continue;
0136  
0137        otherwise; error('huh? shouldnt get here');
0138       end
0139       fprintf(fid,'solid cyl%04d = mainobj    and %s; \n',i,name);
0140    end
0141 
0142    % In some cases, ng can't handle an extra ';'
0143    if ~isempty(extra_ng_code)
0144       fprintf(fid,'%s;\n', extra_ng_code);
0145    end
0146    % SHOULD tank_maxh go here?
0147    fprintf(fid,'tlo mainobj;\n');
0148    for i=1:n_elecs
0149       if any(i == pts_elecs_idx); continue; end
0150       if elecs(i).shape == 'U';   continue; end % USER WILL DEFINE
0151       fprintf(fid,'tlo cyl%04d %s -col=[1,0,0];\n',i, elecs(i).ngobj);
0152    end
0153 
0154 %  if ~isempty(extra_ng_code{1})
0155 %     fprintf(fid,'tlo %s  -col=[0,1,0];\n',extra_ng_code{1});
0156 %  end
0157 
0158    fclose(fid); % geofn
0159 % From Documentation: Syntax is
0160 % np
0161 % x1 y1 z1 h1
0162 % x2 y2 z2 h2
0163    n_pts_elecs= length(pts_elecs_idx) + size(mszpoints,1);
0164    fid=fopen(ptsfn,'w');
0165    fprintf(fid,'%d\n',n_pts_elecs);
0166    for i = pts_elecs_idx;
0167       posxy = elecs(i).pos(1:2);
0168       fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) );
0169    end
0170    for i = 1:size(mszpoints,1);
0171       fprintf(fid,'%10f %10f %10f %10f\n', mszpoints(i,:));
0172    end
0173    fprintf(fid,'0\n'); % lines
0174    fclose(fid); % ptsfn
0175 
0176 
0177 % ELECTRODE POSITIONS:
0178 %  elec_pos = [n_elecs_per_plane,z_planes]
0179 %     OR
0180 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0181 %
0182 % ELECTRODE SHAPES::
0183 %  elec_shape = [width,height, {maxsz}]  % Rectangular elecs
0184 %     OR
0185 %  elec_shape = [radius, {0, maxsz} ]  % Circular elecs
0186 %     maxsz  (OPT)  -> max size of mesh elems (default = courase mesh)
0187 %
0188 % OUTPUT:
0189 %  elecs(i).pos   = [x,y,z]
0190 %  elecs(i).shape = 'C' or 'R'
0191 %  elecs(i).dims  = [radius] or [width,height]
0192 %  elecs(i).maxh  = '-maxh=#' or '';
0193 function [elecs, centres] = parse_elecs( elec_pos, elec_shape, is2D, elec_obj );
0194 
0195    n_elecs= size(elec_pos,1); 
0196    if n_elecs == 0
0197       elecs= struct([]); % empty
0198       centres= [];
0199       return;
0200    end
0201 
0202 
0203    if size(elec_shape,1) == 1
0204       elec_shape = ones(n_elecs,1) * elec_shape;
0205    end
0206    
0207    centres = elec_pos(:,1:3);
0208 
0209    for i= 1:n_elecs
0210      if ischar(elec_obj); 
0211         elecs(i) = elec_spec( elec_shape(i,:), elec_pos(i,:), elec_obj );
0212      else
0213         elecs(i) = elec_spec( elec_shape(i,:), elec_pos(i,:), elec_obj{i}  );
0214      end
0215    end
0216 
0217 function elec = elec_spec( row, posrow, elec_obj );
0218   elec.pos =  posrow(1:3);
0219   elec.dirn=  posrow(4:6);
0220 
0221   if all(isnan(elec.dirn))
0222      elec.shape = 'U'; %user defined
0223      elec.dims = NaN;
0224      elec.ngobj = '';
0225      elec.maxh = '';
0226      elec = orderfields(elec); % UNBELIEVABLY STUPID MATLAB
0227      return
0228   end
0229 
0230   elec.ngobj= elec_obj;
0231 
0232   if row(1) == 0
0233      elec.shape = 'P';
0234      elec.dims  = row(2)*[1,1];
0235   elseif length(row)<2 || row(2) == 0 % Circular electrodes
0236      elec.shape = 'C';
0237      elec.dims  = row(1);
0238   elseif row(2)>0      % Rectangular electrodes
0239      elec.shape = 'R';
0240      elec.dims  = row(1:2);
0241   else
0242      error('negative electrode width');
0243   end
0244 
0245   if length(row)>=3 && row(3) > 0
0246      elec.maxh = sprintf('-maxh=%f', row(3));
0247   else
0248      elec.maxh = '';
0249   end
0250 
0251      elec = orderfields(elec); % UNBELIEVABLY STUPID MATLAB
0252 
0253 function write_header(fid, shape_str);
0254    fprintf(fid,'#Automatically generated by ng_mk_gen_models\n');
0255    fprintf(fid,'algebraic3d\n');
0256    fprintf(fid,shape_str);
0257 
0258 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
0259 % writes the specification for a netgen cuboid on fid, named name, centerd on c,
0260 % in the direction given by vector dirn,
0261 % hw = [height, width]  and depth d
0262 % direction is in the xy plane
0263 %
0264    d= min(d);
0265    w = wh(1); h= wh(2);
0266    dirn = dirn/norm(dirn);
0267 % get perpendicular vector by cross product with non-colinear vector
0268 if norm(abs(dirn)-[0,0,1])>1e-1
0269    dirnp = cross(dirn,[0,0,1]);
0270 else
0271 % For electrodes oriented vertically, we arbitrarily choose
0272 %   width = x axis, height = y axis
0273    dirnp = cross(dirn,[0,1,0]);
0274 end
0275    dirnp = dirnp/norm(dirnp);
0276    dirnk = cross(dirn,dirnp); % the other surface normal
0277    if any(isnan(dirnk)); warning('normal calc weird'); keyboard; end % debugging
0278 
0279    fprintf(fid,'solid %s  = ', name);
0280    fprintf(fid,' plane (%f,%f,%f;%f,%f,%f) and\n',c+d/2*dirn , dirn );
0281    fprintf(fid,' plane (%f,%f,%f;%f,%f,%f) and\n',c-d/2*dirn ,-dirn );
0282    fprintf(fid,' plane (%f,%f,%f;%f,%f,%f) and\n',c+w/2*dirnp, dirnp);
0283    fprintf(fid,' plane (%f,%f,%f;%f,%f,%f) and\n',c-w/2*dirnp,-dirnp);
0284    fprintf(fid,' plane (%f,%f,%f;%f,%f,%f) and\n',c+h/2*dirnk,+dirnk);
0285    fprintf(fid,' plane (%f,%f,%f;%f,%f,%f) %s;\n',c-h/2*dirnk,-dirnk,maxh);
0286 
0287 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh)
0288 % writes the specification for a netgen cylindrical rod on fid,
0289 %  named name, centerd on c,
0290 % in the direction given by vector d, radius rd  lenght ln
0291 % direction is in the xy plane
0292 % the direction vector
0293    dirn = dirn/norm(dirn);
0294 
0295    ln = min(ln);
0296  % I would divide by 2 here (shorted tube in cyl), but ng doesn't like
0297  % That - it fails for 16 (but no 15 or 17) electrodes
0298    inpt = c - dirn.*(ln/1);
0299    outpt =c + dirn.*(ln/1);
0300 
0301    fprintf(fid,'solid %s  = ', name);
0302    fprintf(fid,'  plane(%f,%f,%f;%f,%f,%f) and\n',       inpt, -dirn);
0303    fprintf(fid,'  plane(%f,%f,%f;%f,%f,%f) and\n',       outpt, dirn);
0304    fprintf(fid,'  cylinder(%f,%f,%f;%f,%f,%f;%f) %s;\n', inpt, outpt, rd,maxh);
0305 
0306 
0307 function electrode = pem_from_cem(elecs, electrode, nodes)
0308 % elecs = electrode structure of model, from the parse_elecs function
0309 % electrode = the forward electrode model
0310 % nodes = the coordinates for the nodes
0311 % Can only have one node per electrode so we get a Point Electrode Model.
0312 % Choose the node with the greatest angle, so we atlest pick a consistent
0313 % side of the electrode: NetGen seems to give a random order to the nodes
0314 % in the electrode listing so we can't just pick the first one.
0315 % The nodes aside from those on the edges are not garanteed to be at any
0316 % particular location, so won't be consistent between meshes.
0317 % TODO should probably also adjust contact impedance too: its found later
0318 % by taking the average of the edges around the PEM's node, and those
0319 % will vary for each mesh -- should adjust so all electrodes get a
0320 % consistent effective impedance later.
0321   Ne = length(electrode);
0322   for i = 1:Ne
0323     if elecs(i).shape == 'P'
0324       % find the angles of the nodes for this electrode relative to (0,0)
0325       xy = nodes(electrode(i).nodes,:);
0326       ang = atan2(xy(:,2),xy(:,1));
0327       % if the angles cover more than 180 degrees, must be an angle
0328       % roll-over from -pi to +pi, so take all the negative angles
0329       % and move them up
0330       if (max(ang) - min(ang)) > pi
0331         ang = ang + (ang <0)*2*pi;
0332       end
0333       % choose the counter-clockwise most node only
0334       if size(xy,2) == 3 ; ang = ang - xy(:,3); end
0335       [jnk, ind] = max(ang);
0336       electrode(i).nodes = electrode(i).nodes(ind);
0337     end
0338   end
0339 
0340 function square_elec_test
0341     % problem is how to define sides of vertical electrode
0342     shape_str = ['solid top    = plane(0,0,  0;0,0, 1);\n' ...
0343                  'solid bot    = plane(0,0,-10;0,0,-1);\n' ...
0344                  'solid ob     = orthobrick(-3,-3,-99;3,3, 99);\n' ...
0345                  'solid mainobj= top and bot and ob -maxh=0.5;\n'];
0346     elec_pos = [ 0,0,0,0,0,1; 1,0,0,0,0,1];
0347     elec_shape=[0.2,0.2,0.05]; elec_obj = {'top','top'};
0348     fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0349     show_fem(fmdl);
0350 
0351 
0352 function do_unit_test
0353   square_elec_test
0354 
0355   for tn = 1:16
0356      eidors_msg('ng_mk_gen_models: unit_test %02d',tn,1);
0357      fmdl= do_test_number(tn);
0358      show_fem(fmdl); drawnow
0359   end
0360 
0361 function fmdl= do_test_number(tn)
0362 switch tn
0363    case 1;
0364  shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0365               'solid mainobj= orthobrick(-2,-2,0;2,2,2) and cyl -maxh=0.3;\n'];
0366  elec_pos = []; elec_shape = []; elec_obj = {};
0367  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0368 
0369    case 2;
0370  shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0371               'solid mainobj= plane(0,0,0;0,0,-1)\n' ...
0372                         'and  plane(0,0,2;0,0,1)\n' ...
0373                         'and  cyl -maxh=0.3;\n'];
0374  elec_pos = [  1,  0,  1,   1,  0,  0;
0375                0,  1,1.2,   0,  1,  0]; 
0376  elec_shape=[0.1];
0377  elec_obj = 'cyl';
0378  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0379 
0380    case 3;
0381  shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0382               'solid mainobj= orthobrick(-2,-2,0;2,2,2) and cyl -maxh=0.3;\n'];
0383  th = linspace(0,2*pi,15)'; th(end) = [];
0384  cs = [cos(th), sin(th)];
0385  elec_pos = [  cs, th/2/pi + 0.5, cs, 0*th];
0386  elec_shape=[0.1];
0387  elec_obj = 'cyl';
0388  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0389 
0390    case 4;
0391  shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0392               'solid mainobj= orthobrick(-2,-2,0;2,2,2) and cyl -maxh=0.3;\n'];
0393  th = linspace(0,2*pi,15)'; th(end) = [];
0394  cs = [cos(th), sin(th)];
0395  elec_pos = [  cs, th/2/pi + 0.5, cs, 0*th];
0396  elec_shape=[0.1*th/2/pi + 0.05];
0397  elec_obj = 'cyl';
0398  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0399 
0400    case 5;
0401  shape_str = ['solid cyl    = cylinder (0,0,0; 1,0,0; 1); \n', ...
0402               'solid mainobj= plane(0,0,0;-1,0,0)\n' ...
0403                         'and  plane(2,0,0;1,0,0)\n' ...
0404                         'and  cyl -maxh=0.3;\n'];
0405  elec_pos = [  1,  0,  1,   0,  0,  1;
0406              1.2,  1,  0,   0,  1,  0]; 
0407  elec_shape=[0.1];
0408  elec_obj = 'cyl';
0409  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0410 
0411    case 6;
0412  shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0413               'solid bottom = plane(0,0,0;0,0,-1);\n' ...
0414               'solid top    = plane(0,0,2;0,0,1);\n' ...
0415               'solid mainobj= top and bottom and cyl -maxh=0.3;\n'];
0416  elec_pos = [  1,  0,  1,   1,  0,  0;
0417                0,  1,1.2,   0,  1,  0;
0418                0.8,  0,  0, 0,  0, -1]; 
0419  elec_shape=[0.1];
0420  elec_obj = {'cyl','cyl','bottom'};
0421  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0422 
0423    case 7;
0424  shape_str = ['solid top    = plane(0,0,0;0,0,1);\n' ...
0425               'solid mainobj= top and orthobrick(-2,-2,-2;2,2,0);\n'];
0426  elec_pos = [  1,  0,  0,   0,  0,  1;
0427                0,  0,  0,   0,  0,  1;
0428               -1,  0,  0,   0,  0,  1];
0429  elec_shape=[0.1];
0430  elec_obj = 'top';
0431  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0432 
0433    case 8;
0434  shape_str = ['solid top    = plane(0,0,0;0,0,1);\n' ...
0435               'solid cyl    = ellipticcylinder(0,0,0;2.5,0,0;0,1,0);\n' ...
0436               'solid mainobj= top and cyl and orthobrick(-2,-2,-2;2,2,0);\n'];
0437  elec_pos = [  1,  0,  0,   0,  0,  1;
0438                0,  0,  0,   0,  0,  1;
0439               -1,  0,  0,   0,  0,  1;
0440                1, -1,-1.2,  0, -1,  0;
0441                0, -1,-1.0,  0, -1,  0;
0442               -1, -1,-0.8,  0, -1,  0];
0443  elec_shape=[0.1];
0444  elec_obj = {'top','top','top','cyl','cyl','cyl'};
0445  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0446 
0447    case 9;
0448  shape_str = ['solid top    = plane(0,0,0;0,0,1);\n' ...
0449               'solid ball   = sphere(-1.25,0,-1;0.5); tlo ball;\n' ...
0450               'solid mainobj= top and orthobrick(-2,-1,-2;2,1,0) and not ball -maxh=0.5;\n'];
0451  elec_pos = linspace( -1.5,1.5,5)';
0452  elec_pos = [  elec_pos, elec_pos*[0,0,0,0], elec_pos*0+1];
0453  elec_shape=[0.3];
0454  elec_obj = 'top';
0455  [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0456  img = mk_image( fmdl, 1);
0457  img.elem_data(mat_idx{2}) = 1.1; 
0458  
0459  fmdl = img; % so that the code shows the image
0460 
0461   case 10;
0462  shape_str = ['solid top    = plane(0,0,0;0,0,1);\n' ...
0463               'solid mainobj= top and orthobrick(-3,-3,-2;3,3,0) -maxh=0.5;\n'];
0464  [elec_pos_x,elec_pos_y] = meshgrid(linspace( -1.5,1.5,5),linspace(-2,2,7));
0465  elec_pos = [  elec_pos_x(:), elec_pos_y(:), ones(size(elec_pos_x(:)))*[0,0,0,1] ];
0466  elec_shape=[0.2];
0467  elec_obj = 'top';
0468  [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0469 
0470    case 11;
0471  shape_str = ['solid cyl    = cylinder (0,0,0; 0,0,1; 1.0); \n', ...
0472               'solid tank   = orthobrick(-2,-2,0;2,2,0.4) and cyl; \n', ...
0473               'solid fish   = ellipsoid(0.2,0.2,0.2;0.2,0,0;0,0.1,0;0,0,0.1); tlo fish;\n', ...
0474               'solid mainobj= tank and not fish -maxh=0.3;\n'];
0475  n_elec = 7;
0476  th = linspace(0,2*pi,n_elec+1)'; th(end) = [];
0477  cs = [cos(th), sin(th)];
0478  elec_pos = [  cs, 0.2+0*th, cs, 0*th];
0479  elec_shape=[0.05];
0480  for i=1:n_elec; elec_obj{i} = 'cyl'; end
0481  i=i+1;elec_pos(i,:) = [ 0  ,0.2,0.2,-1,0,0]; elec_obj{i} = 'fish';
0482  i=i+1;elec_pos(i,:) = [ 0.4,0.2,0.2, 1,0,0]; elec_obj{i} = 'fish';
0483  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0484 
0485    case 12;
0486 shape_str = ['solid top     = ellipsoid(0,0,0; 0,0,1; 1,0,0; 0,1,0); \n' ...
0487     'solid mainobj= top and orthobrick(-2,-2,0;2,2,2) -maxh=0.1;\n'];
0488 deg2rad = pi/180;
0489 th = (-70:20:70)'*deg2rad;
0490  elec_pos = [0*th,sin(th),cos(th),0*th,sin(th),cos(th); ...
0491              sin(th),0*th,cos(th),sin(th),0*th,cos(th)];
0492  elec_shape=[0.05];
0493  elec_obj = 'top';
0494 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0495 
0496    case 13;
0497  shape_str = ['solid cyl    = cylinder (0,0,0; 0,1,0; 1); \n', ...
0498               'solid bottom = plane(0, 0,0;0,-1,0);\n' ...
0499               'solid top    = plane(0,10,0;0, 1,0);\n' ...
0500               'solid cut1   = plane(0, 4,0;0,-1,0);\n' ...
0501               'solid cut2   = plane(0, 6,0;0, 1,0);\n' ...
0502               'solid ball   = cyl and cut1 and cut2;  tlo ball;\n' ...
0503               'solid mainobj= ( top and (not cut2) and cyl ) or ' ...
0504                       '(bottom      and (not cut1) and cyl ) -maxh=0.8;\n'];
0505  elec_pos = [ 0, 10,  0, 0,  1,  0; 
0506               0,  0,  0, 0, -1,  0]; 
0507  elec_shape=[1.0];
0508  elec_obj = {'top','bottom'};
0509  [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0510  fmdl = mk_image(fmdl,1); 
0511  fmdl.elem_data(mat_idx{2}) = 1.1;
0512 
0513    case 14;
0514  shape_str = [ ...
0515   'solid cyl    = cylinder (0,0,0; 0,0,1; 1); \n', ...
0516   'solid mainobj= plane(0,0,0;0,0,-1)\n' ...
0517         'and  plane(0,0,2;0,0,1)\n' ...
0518         'and  cyl -maxh=0.3;\n' ...
0519   'solid in_elec = sphere(0,-1,1;0.2)' ...
0520         'and not    sphere(0,-1,1;0.15) -maxh=0.05;\n' ...
0521         'solid in_elec0= in_elec  and mainobj;\n' ...
0522         'tlo in_elec0 cyl;\n' ...
0523   'solid out_elec = sphere(0,-1,1;0.4)' ...
0524         'and not    sphere(0,-1,1;0.35) -maxh=0.05;\n' ...
0525         'solid out_elec0= out_elec  and mainobj;\n' ...
0526         'tlo out_elec0 cyl;\n'];
0527  % Find a background electrode (for all) first. This will stop errors in the next
0528  elec_pos = [  0, -1,   0, NaN,NaN,NaN;
0529                1,  0,   1,   1,  0,  0;
0530                0,  1, 1.2,   0,  1,  0;
0531                0, -1, 1.2, NaN,NaN,NaN;
0532                0, -1, 1.4, NaN,NaN,NaN];
0533  elec_shape=[0.1];
0534  elec_obj = 'cyl';
0535  fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0536  % Throw away the first electrode (background). We don't need it
0537  fmdl.electrode = fmdl.electrode(2:end);
0538 
0539    case 15;
0540 d = [-10,-20,-10,10,20,10];
0541 body_str= sprintf([ ...
0542    'solid left   = plane( %f,  0,  0; -1, 0, 0);\n', ...
0543    'solid back   = plane(  0, %f,  0;  0,-1, 0);\n', ...
0544    'solid bot    = plane(  0,  0, %f;  0, 0,-1);\n', ...
0545    'solid right  = plane( %f,  0,  0;  1, 0, 0);\n', ...
0546    'solid front  = plane(  0, %f,  0;  0, 1, 0);\n', ...
0547    'solid top    = plane(  0,  0, %f;  0, 0, 1);\n', ...
0548    'solid mainobj= top and bot and left and right and back and front -maxh=2;\n', ...
0549     ], d);
0550 elec_pos = [ ...
0551    d(1), 0,    0,     1,  0,  0;
0552    0,    d(2), 0,     0,  1,  0;
0553    0,    0,    d(3),  0,  0,  1;
0554    d(4), 0,    0,    -1,  0,  0;
0555    0,    d(5), 0,     0, -1,  0;
0556    0,    0,    d(6),  0,  0, -1;];
0557 elec_shape = ones(6,1)*[3,1,.5];
0558 elec_shape(4,1:2) = [2,0];
0559 elec_shape(5,1:2) = [1,8];
0560 elec_shape(6,1:2) = [1,3];
0561 elec_obj = {'left','back','bot','right','front','top'};
0562 fmdl = ng_mk_gen_models(body_str, elec_pos, elec_shape, elec_obj);
0563 
0564    case 16;
0565 d = [-10,-20,-10,10,20,10];
0566 body_str= sprintf([ ...
0567    'solid left   = plane( %f,  0,  0; -3,-1, 0);\n', ...
0568    'solid back   = plane(  0, %f,  0;  0,-3,-1);\n', ...
0569    'solid bot    = plane(  0,  0, %f; -1, 0,-3);\n', ...
0570    'solid right  = plane( %f,  0,  0;  3, 1, 0);\n', ...
0571    'solid front  = plane(  0, %f,  0;  0, 3, 1);\n', ...
0572    'solid top    = plane(  0,  0, %f;  1, 0, 3);\n', ...
0573    'solid mainobj= top and bot and left and right and back and front -maxh=2;\n', ...
0574     ], d);
0575 elec_pos = [ ...
0576    d(1), 0,    0,     3,  1,  0;
0577    0,    d(2), 0,     0,  3,  1;
0578    0,    0,    d(3),  1,  0,  3;
0579    d(4), 0,    0,    -3, -1,  0;
0580    0,    d(5), 0,     0, -3, -1;
0581    0,    0,    d(6), -1,  0, -3;];
0582 elec_shape = ones(6,1)*[3,1,.5];
0583 elec_shape(4,1:2) = [2,0];
0584 elec_shape(5,1:2) = [1,8];
0585 elec_shape(6,1:2) = [1,3];
0586 elec_obj = {'left','back','bot','right','front','top'};
0587 fmdl = ng_mk_gen_models(body_str, elec_pos, elec_shape, elec_obj);
0588 
0589 
0590   otherwise;
0591      error('huh?')
0592 end