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