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