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