ng_mk_cyl_models

PURPOSE ^

NG_MAKE_CYL_MODELS: create cylindrical models using netgen

SYNOPSIS ^

function [fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos,elec_shape, extra_ng_code);

DESCRIPTION ^

 NG_MAKE_CYL_MODELS: create cylindrical models using netgen
[fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos, ...
                 elec_shape, extra_ng_code);
 INPUT:
 cyl_shape = {height, [radius, [maxsz]]}
    if height = 0 -> calculate a 2D shape
    radius (OPT)  -> (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 elecs
    (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 cylinder. Radius = 1. No electrodes
   fmdl= ng_mk_cyl_models(3,[0],[]); 
 Simple 2D cylinder. Radius = 2. Set minsize to refine
   fmdl= ng_mk_cyl_models([0,2,.2],[0],[]); 
 3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1
   fmdl= ng_mk_cyl_models(3,[8,1,2],[0.1]); 
 3D cylinder. Radius = 1. 6 circ elecs with elec refinement
   fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0,0.05]); 
 3D cylinder. Radius = 1. 7 rect elecs with no refinement
   fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0.3]); 
 2D cylinder. Radius = 1. 11 rect elecs with refinement
   fmdl= ng_mk_cyl_models(0,[11],[0.2,0,0.05]); 
 2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement
   fmdl= ng_mk_cyl_models([0,1,0.1],[11],[0.2,0,0.02]); 
 2D cylinder. elecs at 0, 90 and 120 degrees
   fmdl= ng_mk_cyl_models(0,[0;90;120],[0.2,0,0.03]); 
 2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees
   fmdl= ng_mk_cyl_models(0,[0;120],[0.4,0,0.01;0.1,0,0.1]); 
 3D cylinder. elecs at 0, 30, 60, 90 in planes
   fmdl= ng_mk_cyl_models(3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]); 
 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_cyl_models(3,el_pos,el_sz); 
 Simple 3D cylinder with a ball
   extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'}
   [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra); 
   img= eidors_obj('image','ball'); img.fwd_model= fmdl;
   img.elem_data(mat_idx{1}) = 1; img.elem_data(mat_idx{2}) = 2;
 3D cylinder with 8 electrodes and cube
   extra={'cube','solid cube = orthobrick(0.5,0.5,0.5;0,0,1.5);'}
   [fmdl,mat_idx]= ng_mk_cyl_models(2,[8,0.5,1.5],[0.1],extra); 
 3D cylinder with inner cylinder
   extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,1;1,1,2) -maxh=0.05;'}
   [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra); 
 2D cylinder with 8 electrodes and hole
   extra={'ball','solid ball = sphere(0.2,0.2,0;0.2) -maxh=0.05;'}
   fmdl= ng_mk_cyl_models(0,[8],[0.1,0,0.05],extra); 
 2D cylinder with 9 electrodes and inner cylinder
   extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,0;1,1,0.05) -maxh=0.03;'}
   fmdl= ng_mk_cyl_models(0,[9],[0.2,0,0.05],extra); 
   img= eidors_obj('image','ball'); img.fwd_model= fmdl;
   ctr = interp_mesh(fmdl); ctr=(ctr(:,1)-0.2).^2 + (ctr(:,2)-0.2).^2;
   img.elem_data = 1 + 0.1*(ctr<0.2^2);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos, ...
0002                   elec_shape, extra_ng_code);
0003 % NG_MAKE_CYL_MODELS: create cylindrical models using netgen
0004 %[fmdl,mat_idx] = ng_mk_cyl_models(cyl_shape, elec_pos, ...
0005 %                 elec_shape, extra_ng_code);
0006 % INPUT:
0007 % cyl_shape = {height, [radius, [maxsz]]}
0008 %    if height = 0 -> calculate a 2D shape
0009 %    radius (OPT)  -> (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 elecs
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 cylinder. Radius = 1. No electrodes
0039 %   fmdl= ng_mk_cyl_models(3,[0],[]);
0040 % Simple 2D cylinder. Radius = 2. Set minsize to refine
0041 %   fmdl= ng_mk_cyl_models([0,2,.2],[0],[]);
0042 % 3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1
0043 %   fmdl= ng_mk_cyl_models(3,[8,1,2],[0.1]);
0044 % 3D cylinder. Radius = 1. 6 circ elecs with elec refinement
0045 %   fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0,0.05]);
0046 % 3D cylinder. Radius = 1. 7 rect elecs with no refinement
0047 %   fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0.3]);
0048 % 2D cylinder. Radius = 1. 11 rect elecs with refinement
0049 %   fmdl= ng_mk_cyl_models(0,[11],[0.2,0,0.05]);
0050 % 2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement
0051 %   fmdl= ng_mk_cyl_models([0,1,0.1],[11],[0.2,0,0.02]);
0052 % 2D cylinder. elecs at 0, 90 and 120 degrees
0053 %   fmdl= ng_mk_cyl_models(0,[0;90;120],[0.2,0,0.03]);
0054 % 2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees
0055 %   fmdl= ng_mk_cyl_models(0,[0;120],[0.4,0,0.01;0.1,0,0.1]);
0056 % 3D cylinder. elecs at 0, 30, 60, 90 in planes
0057 %   fmdl= ng_mk_cyl_models(3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]);
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_cyl_models(3,el_pos,el_sz);
0062 % Simple 3D cylinder with a ball
0063 %   extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'}
0064 %   [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra);
0065 %   img= eidors_obj('image','ball'); img.fwd_model= fmdl;
0066 %   img.elem_data(mat_idx{1}) = 1; img.elem_data(mat_idx{2}) = 2;
0067 % 3D cylinder with 8 electrodes and cube
0068 %   extra={'cube','solid cube = orthobrick(0.5,0.5,0.5;0,0,1.5);'}
0069 %   [fmdl,mat_idx]= ng_mk_cyl_models(2,[8,0.5,1.5],[0.1],extra);
0070 % 3D cylinder with inner cylinder
0071 %   extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,1;1,1,2) -maxh=0.05;'}
0072 %   [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra);
0073 % 2D cylinder with 8 electrodes and hole
0074 %   extra={'ball','solid ball = sphere(0.2,0.2,0;0.2) -maxh=0.05;'}
0075 %   fmdl= ng_mk_cyl_models(0,[8],[0.1,0,0.05],extra);
0076 % 2D cylinder with 9 electrodes and inner cylinder
0077 %   extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,0;1,1,0.05) -maxh=0.03;'}
0078 %   fmdl= ng_mk_cyl_models(0,[9],[0.2,0,0.05],extra);
0079 %   img= eidors_obj('image','ball'); img.fwd_model= fmdl;
0080 %   ctr = interp_mesh(fmdl); ctr=(ctr(:,1)-0.2).^2 + (ctr(:,2)-0.2).^2;
0081 %   img.elem_data = 1 + 0.1*(ctr<0.2^2);
0082 
0083 % (C) Andy Adler, 2009. (C) Alistair Boyle 2013. Licenced under GPL v2 or v3
0084 % $Id: ng_mk_cyl_models.m 6100 2021-03-14 19:10:43Z aadler $
0085 
0086 if ischar(cyl_shape) && strcmp(cyl_shape,'UNIT_TEST'); do_unit_test; return; end
0087 
0088 if nargin < 4; extra_ng_code = {'',''}; end
0089 copt.cache_obj = { cyl_shape, elec_pos, elec_shape, extra_ng_code};
0090 copt.fstr = 'ng_mk_cyl_models';
0091 copt.cache_on_ng_opt = true;
0092 args = {cyl_shape, elec_pos, elec_shape, extra_ng_code};
0093 
0094 fmdl = eidors_cache(@mk_cyl_model, args, copt);
0095 
0096 mat_idx = fmdl.mat_idx;
0097 
0098 function fmdl = mk_cyl_model( cyl_shape, elec_pos, elec_shape, extra_ng_code );
0099 
0100    fnstem = tempname;
0101    geofn= [fnstem,'.geo'];
0102    ptsfn= [fnstem,'.msz'];
0103    meshfn= [fnstem,'.vol'];
0104 
0105    [tank_height, tank_radius, tank_maxh, is2D] = parse_shape(cyl_shape);
0106    [elecs, centres] = parse_elecs( elec_pos, elec_shape,  ...
0107                           tank_height, tank_radius, is2D );
0108 
0109    n_pts = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0110                   tank_maxh, elecs, extra_ng_code);
0111    if n_pts == 0 
0112       call_netgen( geofn, meshfn);
0113    else
0114       call_netgen( geofn, meshfn, ptsfn);
0115    end
0116 
0117    fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []);
0118 
0119 %  delete(geofn); delete(meshfn); delete(ptsfn); % remove temp files
0120    if is2D
0121       fmdl = mdl2d_from3d(fmdl);
0122    end
0123 
0124    % convert CEM to PEM if so configured
0125    % TODO shunt model is unsupported
0126    if isfield(fmdl,'electrode');
0127       fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes);
0128    end
0129 
0130 % for the newest netgen, we can't call msz file unless there are
0131 % actually points in  it (ie. empty msz files break netgen)
0132 function n_pts_elecs = write_geo_file(geofn, ptsfn, tank_height, tank_radius, ...
0133                         tank_maxh, elecs, extra_ng_code);
0134    fid=fopen(geofn,'w');
0135    write_header(fid,tank_height,tank_radius,tank_maxh,extra_ng_code);
0136 
0137    n_elecs = length(elecs);
0138    %  elecs(i).pos   = [x,y,z]
0139    %  elecs(i).shape = 'C' or 'R'
0140    %  elecs(i).dims  = [radius] or [width,height]
0141    %  elecs(i).maxh  = '-maxh=#' or '';
0142    pts_elecs_idx = []; 
0143 
0144    for i=1:n_elecs
0145       name = sprintf('elec%04d',i);
0146       pos = elecs(i).pos;
0147       switch elecs(i).shape
0148        case 'C'
0149          write_circ_elec(fid,name, pos, pos,  ...
0150                elecs(i).dims, tank_radius, elecs(i).maxh);
0151        case 'R'
0152          write_rect_elec(fid,name, pos, pos,  ...
0153                elecs(i).dims, tank_radius, elecs(i).maxh);
0154        case 'P'
0155          % I had the good idea of trying to specify points for the point electrodes,
0156          % but netgen doesn't really listen to this. So instead, it only puts points
0157          % close to where you ask. Instead, specifc a rectangular elec where you want it.
0158          if 0 % OLD technique - keep in case we can figure out netgen better
0159             pts_elecs_idx = [ pts_elecs_idx, i]; 
0160             continue; % DON'T print solid cyl
0161          else
0162             write_rect_elec(fid,name, pos, pos,  ...
0163                   elecs(i).dims, tank_radius, elecs(i).maxh);
0164          end
0165 
0166        otherwise; error('huh? shouldnt get here');
0167       end
0168       fprintf(fid,'solid cyl%04d = bigcyl    and %s; \n',i,name);
0169    end
0170 
0171    % SHOULD tank_maxh go here? - right now it seems to make the whole model refined
0172    fprintf(fid,'tlo bigcyl;\n');
0173    for i=1:n_elecs
0174       if any(i == pts_elecs_idx); continue; end
0175       fprintf(fid,'tlo cyl%04d cyl -col=[1,0,0];\n ',i);
0176    end
0177 
0178    for i=1:length(extra_ng_code)-1
0179       if ~isempty(extra_ng_code{i})
0180          fprintf(fid,'tlo %s  -col=[0,1,0];\n',extra_ng_code{i});
0181       end
0182    end
0183 
0184    fclose(fid); % geofn
0185 % From Documentation: Syntax is
0186 % np
0187 % x1 y1 z1 h1
0188 % x2 y2 z2 h2
0189    n_pts_elecs= length(pts_elecs_idx);
0190    fid=fopen(ptsfn,'w');
0191    fprintf(fid,'%d\n',n_pts_elecs);
0192    for i = pts_elecs_idx;
0193       posxy = elecs(i).pos(1:2);
0194       fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) );
0195    end
0196    fclose(fid); % ptsfn
0197 
0198 function [tank_height, tank_radius, tank_maxh, is2D] = ...
0199               parse_shape(cyl_shape);
0200    tank_height = cyl_shape(1);
0201    tank_radius = 1;
0202    tank_maxh   = 0;
0203    is2D = 0;
0204 
0205    if length(cyl_shape)>1;
0206       tank_radius=cyl_shape(2);
0207    end
0208    if length(cyl_shape)>2; 
0209       tank_maxh  =cyl_shape(3);
0210    end
0211    if tank_height==0;
0212       is2D = 1;
0213 
0214       %Need some width to let netgen work, but not too much so
0215       % that it meshes the entire region
0216       tank_height = tank_radius/5; % initial extimate
0217       if tank_maxh>0
0218          tank_height = min(tank_height,2*tank_maxh);
0219       end
0220    end
0221 
0222 % ELECTRODE POSITIONS:
0223 %  elec_pos = [n_elecs_per_plane,z_planes]
0224 %     OR
0225 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0226 %
0227 % ELECTRODE SHAPES::
0228 %  elec_shape = [width,height, {maxsz}]  % Rectangular elecs
0229 %     OR
0230 %  elec_shape = [radius, {0, maxsz} ]  % Circular elecs
0231 %     maxsz  (OPT)  -> max size of mesh elems (default = courase mesh)
0232 %
0233 % OUTPUT:
0234 %  elecs(i).pos   = [x,y,z]
0235 %  elecs(i).shape = 'C' or 'R'
0236 %  elecs(i).dims  = [radius] or [width,height]
0237 %  elecs(i).maxh  = '-maxh=#' or '';
0238 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, hig, rad, is2D );
0239     
0240     n_elecs= size(elec_pos,1); 
0241     
0242     if n_elecs == 0
0243       elecs= struct([]); % empty
0244       centres= [];
0245       return;
0246     end
0247    
0248    if is2D
0249       elec_pos(:,2) = hig/2;
0250    end
0251 
0252    % It never makes sense to specify only one elec
0253    % So elec_pos means the number of electrodes in this case
0254    if size(elec_pos,1) == 1
0255        % Parse elec_pos = [n_elecs_per_plane,z_planes]
0256       n_elecs= elec_pos(1); % per plane
0257       th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0258 
0259       on_elecs = ones(n_elecs, 1);
0260       el_th = []; 
0261       el_z  = []; 
0262       for i=2:length(elec_pos)
0263         el_th = [el_th; th];
0264         el_z  = [el_z ; on_elecs*elec_pos(i)];
0265       end
0266    else
0267       el_th = elec_pos(:,1)*2*pi/360;
0268       el_z  = elec_pos(:,2);
0269    end
0270       
0271    n_elecs= size(el_z,1); 
0272 
0273    if size(elec_shape,1) == 1
0274       elec_shape = ones(n_elecs,1) * elec_shape;
0275    end
0276 
0277    for i= 1:n_elecs
0278      row = elec_shape(i,:); 
0279      elecs(i) = elec_spec( row, is2D, hig, rad );
0280    end
0281 
0282 %   %MC FIX 05.07.13 - COORDINATE SYSTEM - THETA = 0 AT 3PM AND THETA INCREASE ANTICLOCK
0283 % NO - we will use clockwise coordinate system
0284    centres = [rad*sin(el_th),rad*cos(el_th),el_z];
0285    
0286    for i= 1:n_elecs; elecs(i).pos  = centres(i,:); end
0287 
0288    if n_elecs == 0
0289       elecs= struct([]); % empty
0290    end
0291 
0292 function elec = elec_spec( row, is2D, hig, rad )
0293   if     is2D
0294      if row(1) == 0;
0295         elec.shape = 'P';
0296 % To create a PEM, we make a square and take the corner. This isn't perfect, since
0297 % the elec isn't quite where we asked for it, but that's as good is I can do. I tried
0298 % asking for two rectangles to touch, but that freaks netgen out.
0299         elec.dims  =  [rad/20, hig]; 
0300      else
0301         elec.shape = 'R';
0302         elec.dims  = [row(1),hig];
0303      end
0304   else
0305      if row(1) == 0
0306         elec.shape = 'P' 
0307         elec.dims  = [rad/20, hig/10];
0308      elseif length(row)<2 || row(2) == 0 % Circular electrodes
0309         elec.shape = 'C';
0310         elec.dims  = row(1);
0311      elseif row(2)>0      % Rectangular electrodes
0312         elec.shape = 'R';
0313         elec.dims  = row(1:2);
0314      else
0315         error('negative electrode width');
0316      end
0317   end
0318 
0319   if length(row)>=3 && row(3) > 0
0320      elec.maxh = sprintf('-maxh=%f', row(3));
0321   else
0322      elec.maxh = '';
0323   end
0324 
0325 function write_header(fid,tank_height,tank_radius,maxsz,extra);
0326    if maxsz==0; 
0327       maxsz = '';
0328    else
0329       maxsz = sprintf('-maxh=%f',maxsz);
0330    end
0331 
0332    extra_ng= '';
0333    for i=1:length(extra)-1
0334       if ~isempty( extra{i} )
0335          extra_ng = sprintf(' %s and (not %s) ', ...
0336             extra_ng,extra{i});
0337       end
0338    end
0339 
0340    fprintf(fid,'#Automatically generated by ng_mk_cyl_models\n');
0341    fprintf(fid,'algebraic3d\n');
0342    fprintf(fid,'%s\n',extra{end}); % Define extra stuff here
0343    fprintf(fid,'solid cyl=cylinder (0,0,0;0,0,%6.2f;%6.2f); \n', ...
0344            tank_height, tank_radius);
0345    fprintf(fid,['solid bigcyl= plane(0,0,0;0,0,-1)\n' ...
0346                 'and  plane(0,0,%6.2f;0,0,1)\n' ...
0347                 'and  cyl %s %s;\n'],tank_height,extra_ng,maxsz);  
0348 
0349 
0350 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
0351 % writes the specification for a netgen cuboid on fid, named name, centerd on c,
0352 % in the direction given by vector dirn,
0353 % hw = [height, width]  and depth d
0354 % direction is in the xy plane
0355    w = wh(1); h= wh(2);
0356    dirn(3) = 0; dirn = dirn/norm(dirn);
0357    dirnp = [-dirn(2),dirn(1),0];
0358    dirnp = dirnp/norm(dirnp);
0359 
0360    bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
0361    tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
0362    fprintf(fid,'solid %s  = ', name);
0363    fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
0364            bl(1),bl(2),bl(3));
0365    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0366            bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
0367    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0368            bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
0369    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
0370            tr(1),tr(2),tr(3));
0371    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0372            tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
0373    fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f  )%s;\n', ...
0374            tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0,maxh);
0375 
0376 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh)
0377 % writes the specification for a netgen cylindrical rod on fid,
0378 %  named name, centerd on c,
0379 % in the direction given by vector d, radius rd  lenght ln
0380 % direction is in the xy plane
0381 % the direction vector
0382    dirn(3) = 0; dirn = dirn/norm(dirn);
0383 
0384  % I would divide by 2 here (shorted tube in cyl), but ng doesn't like
0385  % That - it fails for 16 (but no 15 or 17) electrodes
0386    inpt = c - dirn.*(ln/1);
0387    outpt =c + dirn.*(ln/1);
0388 
0389    fprintf(fid,'solid %s  = ', name);
0390    fprintf(fid,'  plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0391          inpt(1),inpt(2),inpt(3),-dirn(1),-dirn(2),-dirn(3));
0392    fprintf(fid,'  plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0393          outpt(1),outpt(2),outpt(3),dirn(1),dirn(2),dirn(3));
0394    fprintf(fid,'  cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) %s;\n', ...
0395          inpt(1),inpt(2),inpt(3),outpt(1),outpt(2),outpt(3), rd,maxh);
0396 
0397 
0398 function electrode = pem_from_cem(elecs, electrode, nodes)
0399 % elecs = electrode structure of model, from the parse_elecs function
0400 % electrode = the forward electrode model
0401 % nodes = the coordinates for the nodes
0402 % Can only have one node per electrode so we get a Point Electrode Model.
0403 % Choose the node with the greatest angle, so we atlest pick a consistent
0404 % side of the electrode: NetGen seems to give a random order to the nodes
0405 % in the electrode listing so we can't just pick the first one.
0406 % The nodes aside from those on the edges are not garanteed to be at any
0407 % particular location, so won't be consistent between meshes.
0408 % TODO should probably also adjust contact impedance too: its found later
0409 % by taking the average of the edges around the PEM's node, and those
0410 % will vary for each mesh -- should adjust so all electrodes get a
0411 % consistent effective impedance later.
0412   Ne = length(electrode);
0413   for i = 1:Ne
0414     if elecs(i).shape == 'P'
0415       % find the angles of the nodes for this electrode relative to (0,0)
0416       xy = nodes(electrode(i).nodes,:);
0417       ang = atan2(xy(:,2),xy(:,1));
0418       % if the angles cover more than 180 degrees, must be an angle
0419       % roll-over from -pi to +pi, so take all the negative angles
0420       % and move them up
0421       if (max(ang) - min(ang)) > pi
0422         ang = ang + (ang <0)*2*pi;
0423       end
0424       % choose the counter-clockwise most node only
0425       % subtract the height so we get bottom left. This is OK, since we have rect elecs
0426       if size(xy,2) == 3 ; ang = ang - xy(:,3); end
0427       [jnk, ind] = max(ang);
0428       electrode(i).nodes = electrode(i).nodes(ind);
0429     end
0430   end
0431 
0432 function do_unit_test
0433   for tn = 1:do_test_number(0)
0434      fprintf('>>> ng_mk_cyl_models: TEST #%d\n',tn);
0435      fmdl= do_test_number(tn);
0436      show_fem(fmdl);
0437   end
0438 
0439 function fmdl= do_test_number(tn)
0440    eidors_msg('ng_mk_cyl_models: UNIT_TEST #%d',tn,1);
0441    switch tn
0442    case 1;
0443 % Simple 3D cylinder. Radius = 1. No electrodes
0444     fmdl= ng_mk_cyl_models(3,[0],[]); 
0445 
0446    case 2;
0447 % Simple 2D cylinder. Radius = 2. Set minsize to refine
0448     fmdl= ng_mk_cyl_models([0,2,.2],[0],[]); 
0449 
0450    case 3;
0451 % 3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1
0452     fmdl= ng_mk_cyl_models(3,[8,1,2],[0.1]); 
0453 
0454    case 4;
0455 % 3D cylinder. Radius = 1. 6 circ elecs with elec refinement
0456     fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0,0.05]); 
0457 
0458    case 5;
0459 % 3D cylinder. Radius = 1. 7 rect elecs with no refinement
0460     fmdl= ng_mk_cyl_models(3,[7,1],[0.2,0.3]); 
0461 
0462    case 6;
0463 % 2D cylinder. Radius = 1. 11 rect elecs with refinement
0464     fmdl= ng_mk_cyl_models(0,[11],[0.2,0,0.05]); 
0465 
0466    case 7;
0467 % 2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement
0468     fmdl= ng_mk_cyl_models([0,1,0.1],[11],[0.2,0,0.02]); 
0469 
0470    case 8;
0471 % 2D cylinder. elecs at 0, 90 and 120 degrees
0472     fmdl= ng_mk_cyl_models(0,[0;90;120],[0.2,0,0.03]); 
0473 
0474    case 9;
0475 % 2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees
0476     fmdl= ng_mk_cyl_models(0,[0;120],[0.4,0,0.01;0.1,0,0.1]); 
0477 
0478    case 10;
0479 % 3D cylinder. elecs at 0, 30, 60, 90 in planes
0480     fmdl= ng_mk_cyl_models(3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]); 
0481 
0482    case 11;
0483 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
0484     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0485     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0486     fmdl= ng_mk_cyl_models(3,el_pos,el_sz); 
0487 
0488    case 12;
0489 % Simple 3D cylinder with a ball
0490     extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'};
0491     fmdl= ng_mk_cyl_models(3,[0],[],extra); 
0492     img= mk_image(fmdl, 1);
0493     img.elem_data(fmdl.mat_idx{2}) = 2;
0494 
0495    case 13;
0496 % 3D cylinder with 8 electrodes and cube
0497     extra={'cube','solid cube = orthobrick(0.5,0.5,0.5;0,0,1.5);'};
0498     [fmdl,mat_idx]= ng_mk_cyl_models(2,[8,0.5,1.5],[0.1],extra); 
0499 
0500    case 14;
0501 % 3D cylinder with inner cylinder
0502     extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,1;1,1,2) -maxh=0.05;'};
0503     [fmdl,mat_idx]= ng_mk_cyl_models(3,[0],[],extra); 
0504 
0505    case 15;
0506 % 2D cylinder with 8 electrodes and hole
0507     extra={'ball','solid ball = sphere(0.2,0.2,0;0.2) -maxh=0.05;'};
0508     fmdl= ng_mk_cyl_models(0,[8],[0.1,0,0.05],extra); 
0509 
0510    case 16;
0511 % 2D cylinder with 9 electrodes and inner cylinder
0512     extra={'ball','solid ball = cylinder(0.2,0.2,0;0.2,0.2,1;0.2) and orthobrick(-1,-1,0;1,1,0.05) -maxh=0.03;'};
0513     fmdl= ng_mk_cyl_models(0,[9],[0.2,0,0.05],extra); 
0514     img= eidors_obj('image','ball'); img.fwd_model= fmdl;
0515     ctr = interp_mesh(fmdl); ctr=(ctr(:,1)-0.2).^2 + (ctr(:,2)-0.2).^2;
0516     img.elem_data = 1 + 0.1*(ctr<0.2^2);
0517 
0518    case 17;
0519 % Simple 3D cylinder with a ball
0520     extra={'ball','solid ball = sphere(0.5,0.5,2;0.4);'};
0521     ng_write_opt('MSZBRICK',[-1,1,-1,1,2,3,0.1]);  
0522     fmdl= ng_mk_cyl_models(3,[0],[],extra); 
0523     delete('ng.opt');
0524     img= mk_image(fmdl, 1);
0525     img.elem_data(fmdl.mat_idx{2}) = 2;
0526 
0527    case 0; fmdl = 17; %%%% RETURN MAXIMUM
0528    otherwise;
0529      error('huh?')
0530    end
0531   
0532    if exist('img'); fmdl = img; end

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005