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

Generated on Tue 31-Dec-2019 17:03:26 by m2html © 2005