ng_mk_ellip_models

 NG_MAKE_ELLIP_MODELS: create elliptical models using netgen
[fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
                 elec_shape, extra_ng_code);
 INPUT:
 ellip_shape = {height, [x_radius, y_radius, [maxsz]]}
    if height = 0 -> calculate a 2D shape
    x_radius, y_radius (OPT)  -> elliptical eccentricity in x,y directions(default = 1)
    maxsz  (OPT)  -> max size of mesh elems (default = course mesh)

 ELECTRODE POSITIONS:
  elec_pos = [n_elecs_per_plane,z_planes] 
     OR
  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)

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

 Specify either a common electrode shape or for each electrode

 EXTRA_NG_CODE
   string of extra code to put into netgen geo file. Normally this
   would be to insert extra materials into the space

 OUTPUT:
  fmdl    - fwd_model object
  mat_idx - indices of materials (if extra_ng_code is used)
    Note mat_idx does not work in 2D. Netgen does not provide it.


 USAGE EXAMPLES:
 Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
 
 Simple 2D cylinder. Axes = [1.5,0.8]. Refined to 0.1
     fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
 
 3D cylinder. Axes = [1.5,0.8]. 2 planes of 8 elecs with radius 0.1
     fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
 
 3D cylinder. Axes= [1.3,1] = 1. 7 rect elecs with no refinement
     fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
 
 2D cylinder. Axes = [1.2,0.8]. 11 rect elecs with refinement. Rotated 45 degrees
     fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]); 
     th = 45* [2*pi/360];
     fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
 
 2D cylinder. elecs at 0, 90 and 120 degrees
     fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
 
 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
     fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
 
 Simple 3D cylinder with a ball
     extra={'ball','solid ball = sphere(0.5,0.5,1;0.4);'};
     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra); 
     img= mk_image(fmdl, 1);
     img.elem_data(mat_idx{2}) = 2;   show_fem(img);

0001 function [fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0002                   elec_shape, extra_ng_code);
0003 % NG_MAKE_ELLIP_MODELS: create elliptical models using netgen
0004 %[fmdl,mat_idx] = ng_mk_ellip_models(ellip_shape, elec_pos, ...
0005 %                 elec_shape, extra_ng_code);
0006 % INPUT:
0007 % ellip_shape = {height, [x_radius, y_radius, [maxsz]]}
0008 %    if height = 0 -> calculate a 2D shape
0009 %    x_radius, y_radius (OPT)  -> elliptical eccentricity in x,y directions(default = 1)
0010 %    maxsz  (OPT)  -> max size of mesh elems (default = course mesh)
0011 %
0012 % ELECTRODE POSITIONS:
0013 %  elec_pos = [n_elecs_per_plane,z_planes]
0014 %     OR
0015 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0016 %
0017 % ELECTRODE SHAPES::
0018 %  elec_shape = [width,height, maxsz]  % Rectangular elecs
0019 %     OR
0020 %  elec_shape = [radius, 0, maxsz ]    % Circular elecs
0021 %     OR
0022 %  elec_shape = [0, 0, maxsz ]         % Point electrodes
0023 %    (point elecs does some tricks with netgen, so the elecs aren't exactly where you ask)
0024 %
0025 % Specify either a common electrode shape or for each electrode
0026 %
0027 % EXTRA_NG_CODE
0028 %   string of extra code to put into netgen geo file. Normally this
0029 %   would be to insert extra materials into the space
0030 %
0031 % OUTPUT:
0032 %  fmdl    - fwd_model object
0033 %  mat_idx - indices of materials (if extra_ng_code is used)
0034 %    Note mat_idx does not work in 2D. Netgen does not provide it.
0035 %
0036 %
0037 % USAGE EXAMPLES:
0038 % Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
0039 %     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
0040 %
0041 % Simple 2D cylinder. Axes = [1.5,0.8]. Refined to 0.1
0042 %     fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
0043 %
0044 % 3D cylinder. Axes = [1.5,0.8]. 2 planes of 8 elecs with radius 0.1
0045 %     fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
0046 %
0047 % 3D cylinder. Axes= [1.3,1] = 1. 7 rect elecs with no refinement
0048 %     fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
0049 %
0050 % 2D cylinder. Axes = [1.2,0.8]. 11 rect elecs with refinement. Rotated 45 degrees
0051 %     fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]);
0052 %     th = 45* [2*pi/360];
0053 %     fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
0054 %
0055 % 2D cylinder. elecs at 0, 90 and 120 degrees
0056 %     fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
0057 %
0058 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
0059 %     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0060 %     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0061 %     fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
0062 %
0063 % Simple 3D cylinder with a ball
0064 %     extra={'ball','solid ball = sphere(0.5,0.5,1;0.4);'};
0065 %     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra);
0066 %     img= mk_image(fmdl, 1);
0067 %     img.elem_data(mat_idx{2}) = 2;   show_fem(img);
0068 
0069 
0070 
0071 % (C) Andy Adler, 2010. (C) Alistair Boyle 2013. Licenced under GPL v2 or v3
0072 % $Id: ng_mk_ellip_models.m 6100 2021-03-14 19:10:43Z aadler $
0073 
0074 if ischar(ellip_shape) && strcmp(ellip_shape,'UNIT_TEST'), do_unit_test, return, end
0075 if nargin < 4; extra_ng_code = {'',''}; end
0076 
0077 copt.cache_obj = { ellip_shape, elec_pos, elec_shape, extra_ng_code};
0078 copt.fstr = 'ng_mk_ellip_models';
0079 args = {ellip_shape, elec_pos, elec_shape, extra_ng_code};
0080 copt.cache_on_ng_opt = true;
0081 
0082 fmdl = eidors_cache(@mk_ellip_model,args,copt);
0083 
0084 mat_idx = fmdl.mat_idx;
0085 
0086 function fmdl = mk_ellip_model( ellip_shape, elec_pos, elec_shape, extra_ng_code );
0087 
0088    fnstem = tempname;
0089    geofn= [fnstem,'.geo'];
0090    ptsfn= [fnstem,'.msz'];
0091    meshfn= [fnstem,'.vol'];
0092 
0093    [tank_height, tank_radius, tank_maxh, is2D] = parse_shape(ellip_shape);
0094    [elecs, centres] = parse_elecs( elec_pos, elec_shape,  ...
0095                           tank_height, tank_radius, is2D );
0096 
0097    n_pts = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0098                   tank_maxh, elecs, extra_ng_code);
0099    if n_pts == 0 
0100       call_netgen( geofn, meshfn);
0101    else
0102       call_netgen( geofn, meshfn, ptsfn);
0103    end
0104 
0105    fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []);
0106 
0107    delete(geofn); delete(meshfn); delete(ptsfn); % remove temp files
0108    if is2D
0109       fmdl = mdl2d_from3d(fmdl);
0110    end
0111 
0112    % convert CEM to PEM if so configured
0113    % TODO shunt model is unsupported
0114    if isfield(fmdl,'electrode');
0115      fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes);
0116    end
0117 
0118 % for the newest netgen, we can't call msz file unless there are actually points in  it
0119 function n_pts_elecs = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0120                         tank_maxh, elecs, extra_ng_code);
0121    fid=fopen(geofn,'w');
0122    write_header(fid,tank_height,tank_radius,tank_maxh,extra_ng_code);
0123 
0124    n_elecs = length(elecs);
0125    %  elecs(i).pos   = [x,y,z]
0126    %  elecs(i).shape = 'C' or 'R'
0127    %  elecs(i).dims  = [radius] or [width,height]
0128    %  elecs(i).maxh  = '-maxh=#' or '';
0129    pts_elecs_idx = []; 
0130 
0131    for i=1:n_elecs
0132       name = sprintf('elec%04d',i);
0133       pos = elecs(i).pos;
0134       % calculate the normal vector to the shape
0135       ab = tank_radius(1)/tank_radius(2);
0136       dirn= pos.*[inv(ab), ab, 0 ];
0137       switch elecs(i).shape
0138        case 'C'
0139          write_circ_elec(fid,name, pos, dirn,  ...
0140                elecs(i).dims, tank_radius, elecs(i).maxh);
0141        case 'R'
0142          write_rect_elec(fid,name, pos, dirn,  ...
0143                elecs(i).dims, tank_radius, elecs(i).maxh);
0144        case 'P'
0145          if 0 % Netgen doesn't put elecs where you ask
0146             pts_elecs_idx = [ pts_elecs_idx, i]; 
0147             continue; % DON'T print solid cyl
0148          end
0149          write_rect_elec(fid,name, pos, dirn,  ...
0150                elecs(i).dims, tank_radius, elecs(i).maxh);
0151 
0152        otherwise; error('huh? shouldnt get here');
0153       end
0154       fprintf(fid,'solid cyl%04d = %s and not bigcyl; \n',i,name);
0155    end
0156 
0157    % SHOULD tank_maxh go here?
0158    fprintf(fid,'tlo bigcyl;\n');
0159    for i=1:n_elecs
0160       if any(i == pts_elecs_idx); continue; end
0161       fprintf(fid,'tlo cyl%04d cyl -col=[1,0,0];\n ',i);
0162    end
0163 
0164    for i=1:length(extra_ng_code)-1
0165       if ~isempty(extra_ng_code{i})
0166          fprintf(fid,'tlo %s  -col=[0,1,0];\n',extra_ng_code{i});
0167       end
0168    end
0169 
0170    fclose(fid); % geofn
0171 % From Documentation: Syntax is
0172 % np
0173 % x1 y1 z1 h1
0174 % x2 y2 z2 h2
0175    n_pts_elecs= length(pts_elecs_idx);
0176    fid=fopen(ptsfn,'w');
0177    fprintf(fid,'%d\n',n_pts_elecs);
0178    for i = pts_elecs_idx;
0179       posxy = elecs(i).pos(1:2);
0180       fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) );
0181    end
0182    fclose(fid); % ptsfn
0183 
0184 function [tank_height, tank_radius, tank_maxh, is2D] = ...
0185               parse_shape(cyl_shape);
0186    tank_height = cyl_shape(1);
0187    tank_radius = [1,1];
0188    tank_maxh   = 0;
0189    is2D = 0;
0190    lcs = length(cyl_shape);
0191 
0192    if lcs == 2
0193       tank_radius(1)=cyl_shape(2);
0194    elseif lcs >= 3
0195       tank_radius=cyl_shape(2:3);
0196       if diff(tank_radius) == 0;
0197          warning(['Using ng_mk_ellip_models to create cylinder. This may fail for '...
0198                   'even electrode numbers. Recommend use ng_mk_cyl_models']);
0199       end
0200    end
0201    if length(cyl_shape)>=4; 
0202       tank_maxh  =cyl_shape(4);
0203    end
0204    if tank_height==0;
0205       is2D = 1;
0206 
0207       %Need some width to let netgen work, but not too much so
0208       % that it meshes the entire region
0209       tank_height = min(tank_radius)/5; % initial extimate
0210       if tank_maxh>0
0211          tank_height = min(tank_height,2*tank_maxh);
0212       end
0213    end
0214 
0215 % ELECTRODE POSITIONS:
0216 %  elec_pos = [n_elecs_per_plane,z_planes]
0217 %     OR
0218 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0219 %
0220 % ELECTRODE SHAPES::
0221 %  elec_shape = [width,height, {maxsz}]  % Rectangular elecs
0222 %     OR
0223 %  elec_shape = [radius, {0, maxsz} ]  % Circular elecs
0224 %     maxsz  (OPT)  -> max size of mesh elems (default = courase mesh)
0225 %
0226 % OUTPUT:
0227 %  elecs(i).pos   = [x,y,z]
0228 %  elecs(i).shape = 'C' or 'R'
0229 %  elecs(i).dims  = [radius] or [width,height]
0230 %  elecs(i).maxh  = '-maxh=#' or '';
0231 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, hig, rad, is2D );
0232 
0233    if is2D
0234       elec_pos(:,2) = hig/2;
0235    end
0236 
0237    % It never makes sense to specify only one elec
0238    % So elec_pos means the number of electrodes in this case
0239    if size(elec_pos,1) == 1
0240        % Parse elec_pos = [n_elecs_per_plane,z_planes]
0241       n_elecs= elec_pos(1); % per plane
0242       th = ellip_space_elecs( n_elecs, rad );
0243 
0244       on_elecs = ones(n_elecs, 1);
0245       el_th = []; 
0246       el_z  = []; 
0247       for i=2:length(elec_pos)
0248         el_th = [el_th; th];
0249         el_z  = [el_z ; on_elecs*elec_pos(i)];
0250       end
0251    else
0252       el_th = elec_pos(:,1)*2*pi/360;
0253       el_z  = elec_pos(:,2);
0254    end
0255       
0256    n_elecs= size(el_z,1); 
0257 
0258    if size(elec_shape,1) == 1
0259       elec_shape = ones(n_elecs,1) * elec_shape;
0260    end
0261 
0262    for i= 1:n_elecs
0263      row = elec_shape(i,:); 
0264      elecs(i) = elec_spec( row, is2D, hig, rad );
0265    end
0266    
0267 % Electrodes are numbered clockwise from top.
0268    centres = [rad(1)*sin(el_th),rad(2)*cos(el_th),el_z];   
0269    for i= 1:n_elecs; elecs(i).pos  = centres(i,:); end
0270 
0271    if n_elecs == 0
0272       elecs= struct([]); % empty
0273    end
0274 
0275 % equally space n_elecs around an ellipse of outer radius rad(1),rad(2)
0276 function th = ellip_space_elecs( n_elecs, rad )
0277    % The radius is the integral of sqrt((r1*sin(th))^2 + (r2*cos(th))^2)
0278    %  This integral is the incomplete_elliptic_integral(th, 1-r2/r1)*sqrt(r1)
0279    %  Unfortunately, STUPID MATLAB, doesn't have incomplete elliptic integrals
0280    %  by default. So, rather than install a toolkit for it, we integrate numerically.
0281    if n_elecs==0; th=[]; return; end
0282    
0283    th = linspace(0,2*pi, 100*(n_elecs)); th(1)=[]; % Accuracy to 100x spacing
0284    len = cumsum( sqrt( rad(1)*cos(th).^2 + rad(2)*sin(th).^2 ) );
0285    len = len/max(len);
0286    xi = linspace(0,1,n_elecs+1); xi(1)= []; xi(end)=[];
0287    yi = interp1(len,th,xi);
0288 
0289    th= [0;yi(:)];
0290    for exact = 0:3;
0291       eth = exact/2*pi;
0292       ff = abs(th-eth)<1e-3;
0293       th(ff) = eth;
0294    end
0295 
0296 function elec = elec_spec( row, is2D, hig, rad )
0297   if     is2D
0298      if row(1) == 0;
0299         elec.shape = 'P';
0300 % To create a PEM, we make a square and take the corner. This isn't perfect, since
0301 % the elec isn't quite where we asked for it, but that's as good is I can do. I tried
0302 % asking for two rectangles to touch, but that freaks netgen out.
0303         elec.dims  =  [min(rad)/20, hig]; 
0304      else
0305         elec.shape = 'R';
0306         elec.dims  = [row(1),hig];
0307      end
0308   else
0309      if row(1) == 0
0310         elec.shape = 'P' 
0311         elec.dims  = [min(rad)/20, hig/10];
0312      elseif length(row)<2 || row(2) == 0 % Circular electrodes
0313         elec.shape = 'C';
0314         elec.dims  = row(1);
0315      elseif row(2)>0      % Rectangular electrodes
0316         elec.shape = 'R';
0317         elec.dims  = row(1:2);
0318      else
0319         error('negative electrode width');
0320      end
0321   end
0322 
0323   if length(row)>=3 && row(3) > 0
0324      elec.maxh = sprintf('-maxh=%f', row(3));
0325   else
0326      elec.maxh = '';
0327   end
0328 
0329 
0330 function write_header(fid,tank_height,tank_radius,maxsz,extra);
0331    if maxsz==0; 
0332       maxsz = '';
0333    else
0334       maxsz = sprintf('-maxh=%f',maxsz);
0335    end
0336 
0337    extra_ng= '';
0338    for i=1:length(extra)-1
0339       if ~isempty( extra{i} )
0340          extra_ng = sprintf(' %s and (not %s) ', ...
0341             extra_ng,extra{i});
0342       end
0343    end
0344 
0345    fprintf(fid,'#Automatically generated by ng_mk_ellip_models\n');
0346    fprintf(fid,'algebraic3d\n');
0347    fprintf(fid,['solid mainobj_bot= plane(0,0,0;0,0,-1);\n']);
0348    fprintf(fid,['solid mainobj_top= plane(0,0,%f;0,0,1);\n'], ...
0349                  tank_height);
0350    fprintf(fid,'%s\n',extra{end}); % Define extra stuff here
0351    fprintf(fid,'solid cyl=ellipticcylinder (0,0,0;%.4f,0,0;0,%.2f,0); \n', ...
0352             tank_radius);
0353    fprintf(fid,['solid bigcyl= mainobj_top and mainobj_bot and ' ...
0354                 'cyl %s %s;\n'],extra_ng,maxsz);  
0355 
0356 
0357 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
0358 % writes the specification for a netgen cuboid on fid, named name, centerd on c,
0359 % in the direction given by vector dirn,
0360 % hw = [height, width]  and depth d
0361 % direction is in the xy plane
0362    d= min(d);
0363    w = wh(1); h= wh(2);
0364    dirn(3) = 0; dirn = dirn/norm(dirn);
0365    dirnp = [-dirn(2),dirn(1),0];
0366    dirnp = dirnp/norm(dirnp);
0367 
0368    bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
0369    tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
0370    fprintf(fid,'solid %s  = ', name);
0371    fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
0372            bl(1),bl(2),bl(3));
0373    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0374            bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
0375    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0376            bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
0377    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
0378            tr(1),tr(2),tr(3));
0379    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0380            tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
0381    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f  )%s;\n', ...
0382            tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0,maxh);
0383 
0384 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh)
0385 % writes the specification for a netgen cylindrical rod on fid,
0386 %  named name, centerd on c,
0387 % in the direction given by vector d, radius rd  lenght ln
0388 % direction is in the xy plane
0389 % the direction vector
0390    dirn(3) = 0; dirn = dirn/norm(dirn);
0391 
0392    ln = min(ln);
0393  % This is hard to debug here, why does netgen sometime just fail
0394  % fails for 16 (but no 15 or 17) electrodes
0395  % The 'exact' fix seems to fix this, now. Leave comment above to test
0396    inpt = c - dirn.*(ln/6);
0397    outpt =c + dirn.*(ln/6);
0398 
0399    fprintf(fid,'solid %s  = ', name);
0400    fprintf(fid,'  plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0401          inpt(1),inpt(2),inpt(3),-dirn(1),-dirn(2),-dirn(3));
0402    fprintf(fid,'  plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0403          outpt(1),outpt(2),outpt(3),dirn(1),dirn(2),dirn(3));
0404    fprintf(fid,'  cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) %s;\n', ...
0405          inpt(1),inpt(2),inpt(3),outpt(1),outpt(2),outpt(3), rd,maxh);
0406 
0407 
0408 function electrode = pem_from_cem(elecs, electrode, nodes)
0409 % elecs = electrode structure of model, from the parse_elecs function
0410 % electrode = the forward electrode model
0411 % nodes = the coordinates for the nodes
0412 % Can only have one node per electrode so we get a Point Electrode Model.
0413 % Choose the node with the greatest angle, so we atlest pick a consistent
0414 % side of the electrode: NetGen seems to give a random order to the nodes
0415 % in the electrode listing so we can't just pick the first one.
0416 % The nodes aside from those on the edges are not garanteed to be at any
0417 % particular location, so won't be consistent between meshes.
0418 % TODO should probably also adjust contact impedance too: its found later
0419 % by taking the average of the edges around the PEM's node, and those
0420 % will vary for each mesh -- should adjust so all electrodes get a
0421 % consistent effective impedance later.
0422   Ne = length(electrode);
0423   for i = 1:Ne
0424     if elecs(i).shape == 'P'
0425       % find the angles of the nodes for this electrode relative to (0,0)
0426       xy = nodes(electrode(i).nodes,:);
0427       ang = atan2(xy(:,2),xy(:,1));
0428       % if the angles cover more than 180 degrees, must be an angle
0429       % roll-over from -pi to +pi, so take all the negative angles
0430       % and move them up
0431       if (max(ang) - min(ang)) > pi
0432         ang = ang + (ang <0)*2*pi;
0433       end
0434       % choose the counter-clockwise most node only
0435       if size(xy,2) == 3 ; ang = ang - xy(:,3); end
0436       [jnk, ind] = max(ang);
0437       electrode(i).nodes = electrode(i).nodes(ind);
0438     end
0439   end
0440 
0441 
0442 function do_unit_test
0443    sp=1;     subplot(4,4,sp);
0444 % Simple 3D ellipse. Major, minor axes = [1.5, 0.8]. No electrodes
0445     fmdl= ng_mk_ellip_models([1,1.5,0.8],[0],[]);  show_fem(fmdl);
0446 
0447    sp=sp+1;  subplot(4,4,sp);
0448 % Simple 2D cylinder. Axes = [1.5,0.8]. Refined to 0.1
0449     fmdl= ng_mk_ellip_models([0,1.5,0.8,0.1],[],[]); show_fem(fmdl);
0450 
0451    sp=sp+1;  subplot(4,4,sp);
0452 % 3D cylinder. Axes = [1.5,0.8]. 2 planes of 8 elecs with radius 0.1
0453     fmdl= ng_mk_ellip_models([1,1.2,0.8],[8,0.3,0.7],[0.1]); show_fem(fmdl);
0454 
0455    sp=sp+1;  subplot(4,4,sp);
0456 % 3D cylinder. Axes= [1.3,1] = 1. 7 rect elecs with no refinement
0457     fmdl= ng_mk_ellip_models([3,1.3],[7,1],[0.2,0.3]); show_fem(fmdl);
0458 
0459    sp=sp+1;  subplot(4,4,sp);
0460 % 2D cylinder. Axes = [1.2,0.8]. 11 rect elecs with refinement. Rotated 45 degrees
0461     fmdl= ng_mk_ellip_models([0,1.2,0.8],[11],[0.2,0,0.05]); 
0462     th = 45* [2*pi/360];
0463     fmdl.nodes = fmdl.nodes*[cos(th),sin(th);-sin(th),cos(th)]; show_fem(fmdl);
0464 
0465    sp=sp+1;  subplot(4,4,sp);
0466 % 2D cylinder. elecs at 0, 90 and 120 degrees
0467     fmdl= ng_mk_ellip_models([0,1.2,0.8],[0;90;120],[0.2,0,0.03]); show_fem(fmdl);
0468 
0469    sp=sp+1;  subplot(4,4,sp);
0470 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
0471     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0472     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0473     fmdl= ng_mk_ellip_models([3,0.8,1.2],el_pos,el_sz); show_fem(fmdl);
0474 
0475    sp=sp+1;  subplot(4,4,sp);
0476 % Simple 3D cylinder with a ball
0477     extra={'ball','solid ball = sphere(0.3,0.3,1;0.4);'};
0478     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra); 
0479     img= mk_image(fmdl, 1);
0480     img.elem_data(mat_idx{2}) = 2;   show_fem(img);
0481 
0482    sp=sp+1;  subplot(4,4,sp);
0483 % 3D cylinder with a two balls
0484     b1 = 'solid ball1= sphere(0.5, 0.5,1;0.2);';
0485     b2 = 'solid ball2= sphere(0.5,-0.5,1;0.2);';
0486     extra = {'ball1','ball2',[b1,b2]};
0487     [fmdl,mat_idx]= ng_mk_ellip_models([2,1.2,0.8],[8,1],[0.1],extra); 
0488     img= mk_image(fmdl, 1);
0489     img.elem_data(mat_idx{2}) = 2;
0490     img.elem_data(mat_idx{3}) = 0.5;
0491     show_fem(img);
0492      
0493    sp=sp+1;  subplot(4,4,sp);
0494 % Simple 3D cylinder with a ball
0495     extra={'ball','solid ball = sphere(0.3,0.3,1;0.4);'};
0496     [fmdl,mat_idx]= ng_mk_ellip_models([1.15,1.2,0.8],[8,1],[0.1],extra); 
0497     img= mk_image(fmdl, 1);
0498     img.elem_data(mat_idx{2}) = 2;   show_fem(img); view(-30,3);
0499 
0500    sp=sp+1;  subplot(4,4,sp);
0501 % Simple 3D cylinder with a ball
0502     extra={'ball',[ ...
0503        'solid topcut = plane(0,0,1.15;0,0,1);' ...
0504        'solid ball = sphere(0.3,0.3,1;0.4) and topcut;']};
0505     [fmdl,mat_idx]= ng_mk_ellip_models([1.15,1.2,0.8],[8,1],[0.1],extra); 
0506     img= mk_image(fmdl, 1);
0507     img.elem_data(mat_idx{2}) = 2;   show_fem(img); view(-30,3);