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