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