ng_mk_gen_models

PURPOSE ^

NG_MK_GEN_MODELS: create generic models using netgen

SYNOPSIS ^

function [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code);

DESCRIPTION ^

 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);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005