0001 function [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046 if ischar(shape_str) && strcmp(shape_str,'UNIT_TEST'); do_unit_test; return; end
0047
0048 if nargin <= 4; extra_ng_code = ''; end
0049 copt.cache_obj = { shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code };
0050 copt.fstr = 'ng_mk_gen_models';
0051 args = { shape_str, elec_pos, elec_shape, elec_obj, extra_ng_code };
0052
0053 fmdl = eidors_cache(@mk_gen_model,args, copt);
0054
0055 mat_idx = fmdl.mat_idx;
0056
0057 function fmdl = mk_gen_model( shape_str, elec_pos, elec_shape, ...
0058 elec_obj, extra_ng_code);
0059
0060 fnstem = tempname;
0061 geofn= [fnstem,'.geo'];
0062 ptsfn= [fnstem,'.msz'];
0063 meshfn= [fnstem,'.vol'];
0064
0065 is2D = 0;
0066 [elecs, centres] = parse_elecs( elec_pos, elec_shape, is2D, elec_obj );
0067
0068 n_pts = write_geo_file(geofn, ptsfn, shape_str, elecs, extra_ng_code);
0069 if n_pts == 0
0070 call_netgen( geofn, meshfn);
0071 else
0072 call_netgen( geofn, meshfn, ptsfn);
0073 end
0074
0075 fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', []);
0076
0077 delete(geofn); delete(meshfn); delete(ptsfn);
0078 if is2D
0079 fmdl = mdl2d_from3d(fmdl);
0080 end
0081
0082
0083
0084 if isfield(fmdl,'electrode');
0085 fmdl.electrode = pem_from_cem(elecs, fmdl.electrode, fmdl.nodes);
0086 end
0087
0088
0089 function n_pts_elecs = write_geo_file(geofn, ptsfn, shape_str, ...
0090 elecs, extra_ng_code);
0091 fid=fopen(geofn,'w');
0092 write_header(fid, shape_str);
0093
0094 n_elecs = length(elecs);
0095
0096
0097
0098
0099 pts_elecs_idx = [];
0100
0101 if n_elecs > 1
0102 elec_depth = min(nonzeros(distmat(vertcat( elecs(:).pos ))))/2;
0103
0104
0105
0106 end
0107 for i=1:n_elecs
0108 name = sprintf('elec%04d',i);
0109 pos = elecs(i).pos;
0110 dirn= elecs(i).dirn;
0111 switch elecs(i).shape
0112 case 'C'
0113 write_circ_elec(fid,name, pos, dirn, ...
0114 elecs(i).dims, elec_depth, elecs(i).maxh);
0115 case 'R'
0116 write_rect_elec(fid,name, pos, dirn, ...
0117 elecs(i).dims, elec_depth, elecs(i).maxh);
0118 case 'P'
0119 if 0
0120 pts_elecs_idx = [ pts_elecs_idx, i];
0121 continue;
0122 end
0123 write_rect_elec(fid,name, pos, dirn, ...
0124 elecs(i).dims, elec_depth, elecs(i).maxh);
0125
0126 case 'U'
0127 eidors_msg('user defined electrode %d',i, 4);
0128 continue;
0129
0130 otherwise; error('huh? shouldnt get here');
0131 end
0132 fprintf(fid,'solid cyl%04d = mainobj and %s; \n',i,name);
0133 end
0134
0135
0136 if ~isempty(extra_ng_code)
0137 fprintf(fid,'%s;\n', extra_ng_code);
0138 end
0139
0140 fprintf(fid,'tlo mainobj;\n');
0141 for i=1:n_elecs
0142 if any(i == pts_elecs_idx); continue; end
0143 if elecs(i).shape == 'U'; continue; end
0144 fprintf(fid,'tlo cyl%04d %s -col=[1,0,0];\n',i, elecs(i).ngobj);
0145 end
0146
0147
0148
0149
0150
0151 fclose(fid);
0152
0153
0154
0155
0156 n_pts_elecs= length(pts_elecs_idx);
0157 fid=fopen(ptsfn,'w');
0158 fprintf(fid,'%d\n',n_pts_elecs);
0159 for i = pts_elecs_idx;
0160 posxy = elecs(i).pos(1:2);
0161 fprintf(fid,'%10f %10f 0 %10f\n', posxy, elecs(i).dims(1) );
0162 end
0163 fclose(fid);
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182 function [elecs, centres] = parse_elecs( elec_pos, elec_shape, is2D, elec_obj );
0183
0184 n_elecs= size(elec_pos,1);
0185 if n_elecs == 0
0186 elecs= struct([]);
0187 centres= [];
0188 return;
0189 end
0190
0191
0192 if size(elec_shape,1) == 1
0193 elec_shape = ones(n_elecs,1) * elec_shape;
0194 end
0195
0196 centres = elec_pos(:,1:3);
0197
0198 for i= 1:n_elecs
0199 if ischar(elec_obj);
0200 elecs(i) = elec_spec( elec_shape(i,:), elec_pos(i,:), elec_obj );
0201 else
0202 elecs(i) = elec_spec( elec_shape(i,:), elec_pos(i,:), elec_obj{i} );
0203 end
0204 end
0205
0206 function elec = elec_spec( row, posrow, elec_obj );
0207 elec.pos = posrow(1:3);
0208 elec.dirn= posrow(4:6);
0209
0210 if all(isnan(elec.dirn))
0211 elec.shape = 'U';
0212 elec.dims = NaN;
0213 elec.ngobj = '';
0214 elec.maxh = '';
0215 elec = orderfields(elec);
0216 return
0217 end
0218
0219 elec.ngobj= elec_obj;
0220
0221 if row(1) == 0
0222 elec.shape = 'P';
0223 elec.dims = row(2)*[1,1];
0224 elseif length(row)<2 || row(2) == 0
0225 elec.shape = 'C';
0226 elec.dims = row(1);
0227 elseif row(2)>0
0228 elec.shape = 'R';
0229 elec.dims = row(1:2);
0230 else
0231 error('negative electrode width');
0232 end
0233
0234 if length(row)>=3 && row(3) > 0
0235 elec.maxh = sprintf('-maxh=%f', row(3));
0236 else
0237 elec.maxh = '';
0238 end
0239
0240 elec = orderfields(elec);
0241
0242 function write_header(fid, shape_str);
0243 fprintf(fid,'#Automatically generated by ng_mk_gen_models\n');
0244 fprintf(fid,'algebraic3d\n');
0245 fprintf(fid,shape_str);
0246
0247 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
0248
0249
0250
0251
0252
0253
0254
0255 d= min(d);
0256 w = wh(1); h= wh(2);
0257 dirn = dirn/norm(dirn);
0258 dirnp = [-dirn(2),dirn(1),0];
0259
0260 if norm(dirnp)<1e-6
0261 dirnp = [0,1,0];
0262 dirn = [1,0,0];
0263 else
0264 dirnp = dirnp/norm(dirnp);
0265 end
0266
0267 bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
0268 tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
0269 fprintf(fid,'solid %s = ', name);
0270 fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
0271 bl(1),bl(2),bl(3));
0272 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0273 bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
0274 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0275 bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
0276 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
0277 tr(1),tr(2),tr(3));
0278 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
0279 tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
0280 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f )%s;\n', ...
0281 tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0,maxh);
0282
0283 function write_circ_elec(fid,name,c, dirn,rd,ln,maxh)
0284
0285
0286
0287
0288
0289 dirn = dirn/norm(dirn);
0290
0291 ln = min(ln);
0292
0293
0294 inpt = c - dirn.*(ln/1);
0295 outpt =c + dirn.*(ln/1);
0296
0297 fprintf(fid,'solid %s = ', name);
0298 fprintf(fid,' plane(%f,%f,%f;%f,%f,%f) and\n', inpt, -dirn);
0299 fprintf(fid,' plane(%f,%f,%f;%f,%f,%f) and\n', outpt, dirn);
0300 fprintf(fid,' cylinder(%f,%f,%f;%f,%f,%f;%f) %s;\n', inpt, outpt, rd,maxh);
0301
0302
0303 function electrode = pem_from_cem(elecs, electrode, nodes)
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317 Ne = length(electrode);
0318 for i = 1:Ne
0319 if elecs(i).shape == 'P'
0320
0321 xy = nodes(electrode(i).nodes,:);
0322 ang = atan2(xy(:,2),xy(:,1));
0323
0324
0325
0326 if (max(ang) - min(ang)) > pi
0327 ang = ang + (ang <0)*2*pi;
0328 end
0329
0330 if size(xy,2) == 3 ; ang = ang - xy(:,3); end
0331 [jnk, ind] = max(ang);
0332 electrode(i).nodes = electrode(i).nodes(ind);
0333 end
0334 end
0335
0336 function square_elec_test
0337
0338 shape_str = ['solid top = plane(0,0, 0;0,0, 1);\n' ...
0339 'solid bot = plane(0,0,-10;0,0,-1);\n' ...
0340 'solid ob = orthobrick(-3,-3,-99;3,3, 99);\n' ...
0341 'solid mainobj= top and bot and ob -maxh=0.5;\n'];
0342 elec_pos = [ 0,0,0,0,0,1; 1,0,0,0,0,1];
0343 elec_shape=[0.2,0.2,0.05]; elec_obj = {'top','top'};
0344 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0345 show_fem(fmdl);
0346
0347
0348 function do_unit_test
0349 square_elec_test
0350
0351 for tn = 1:14
0352 eidors_msg('ng_mk_gen_models: unit_test %02d',tn,1);
0353 fmdl= do_test_number(tn);
0354 show_fem(fmdl); drawnow
0355 end
0356
0357 function fmdl= do_test_number(tn)
0358 switch tn
0359 case 1;
0360 shape_str = ['solid cyl = cylinder (0,0,0; 0,0,1; 1); \n', ...
0361 'solid mainobj= orthobrick(-2,-2,0;2,2,2) and cyl -maxh=0.3;\n'];
0362 elec_pos = []; elec_shape = []; elec_obj = {};
0363 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0364
0365 case 2;
0366 shape_str = ['solid cyl = cylinder (0,0,0; 0,0,1; 1); \n', ...
0367 'solid mainobj= plane(0,0,0;0,0,-1)\n' ...
0368 'and plane(0,0,2;0,0,1)\n' ...
0369 'and cyl -maxh=0.3;\n'];
0370 elec_pos = [ 1, 0, 1, 1, 0, 0;
0371 0, 1,1.2, 0, 1, 0];
0372 elec_shape=[0.1];
0373 elec_obj = 'cyl';
0374 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0375
0376 case 3;
0377 shape_str = ['solid cyl = cylinder (0,0,0; 0,0,1; 1); \n', ...
0378 'solid mainobj= orthobrick(-2,-2,0;2,2,2) and cyl -maxh=0.3;\n'];
0379 th = linspace(0,2*pi,15)'; th(end) = [];
0380 cs = [cos(th), sin(th)];
0381 elec_pos = [ cs, th/2/pi + 0.5, cs, 0*th];
0382 elec_shape=[0.1];
0383 elec_obj = 'cyl';
0384 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0385
0386 case 4;
0387 shape_str = ['solid cyl = cylinder (0,0,0; 0,0,1; 1); \n', ...
0388 'solid mainobj= orthobrick(-2,-2,0;2,2,2) and cyl -maxh=0.3;\n'];
0389 th = linspace(0,2*pi,15)'; th(end) = [];
0390 cs = [cos(th), sin(th)];
0391 elec_pos = [ cs, th/2/pi + 0.5, cs, 0*th];
0392 elec_shape=[0.1*th/2/pi + 0.05];
0393 elec_obj = 'cyl';
0394 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0395
0396 case 5;
0397 shape_str = ['solid cyl = cylinder (0,0,0; 1,0,0; 1); \n', ...
0398 'solid mainobj= plane(0,0,0;-1,0,0)\n' ...
0399 'and plane(2,0,0;1,0,0)\n' ...
0400 'and cyl -maxh=0.3;\n'];
0401 elec_pos = [ 1, 0, 1, 0, 0, 1;
0402 1.2, 1, 0, 0, 1, 0];
0403 elec_shape=[0.1];
0404 elec_obj = 'cyl';
0405 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0406
0407 case 6;
0408 shape_str = ['solid cyl = cylinder (0,0,0; 0,0,1; 1); \n', ...
0409 'solid bottom = plane(0,0,0;0,0,-1);\n' ...
0410 'solid top = plane(0,0,2;0,0,1);\n' ...
0411 'solid mainobj= top and bottom and cyl -maxh=0.3;\n'];
0412 elec_pos = [ 1, 0, 1, 1, 0, 0;
0413 0, 1,1.2, 0, 1, 0;
0414 0.8, 0, 0, 0, 0, -1];
0415 elec_shape=[0.1];
0416 elec_obj = {'cyl','cyl','bottom'};
0417 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0418
0419 case 7;
0420 shape_str = ['solid top = plane(0,0,0;0,0,1);\n' ...
0421 'solid mainobj= top and orthobrick(-2,-2,-2;2,2,0);\n'];
0422 elec_pos = [ 1, 0, 0, 0, 0, 1;
0423 0, 0, 0, 0, 0, 1;
0424 -1, 0, 0, 0, 0, 1];
0425 elec_shape=[0.1];
0426 elec_obj = 'top';
0427 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0428
0429 case 8;
0430 shape_str = ['solid top = plane(0,0,0;0,0,1);\n' ...
0431 'solid cyl = ellipticcylinder(0,0,0;2.5,0,0;0,1,0);\n' ...
0432 'solid mainobj= top and cyl and orthobrick(-2,-2,-2;2,2,0);\n'];
0433 elec_pos = [ 1, 0, 0, 0, 0, 1;
0434 0, 0, 0, 0, 0, 1;
0435 -1, 0, 0, 0, 0, 1;
0436 1, -1,-1.2, 0, -1, 0;
0437 0, -1,-1.0, 0, -1, 0;
0438 -1, -1,-0.8, 0, -1, 0];
0439 elec_shape=[0.1];
0440 elec_obj = {'top','top','top','cyl','cyl','cyl'};
0441 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0442
0443 case 9;
0444 shape_str = ['solid top = plane(0,0,0;0,0,1);\n' ...
0445 'solid ball = sphere(-1.25,0,-1;0.5); tlo ball;\n' ...
0446 'solid mainobj= top and orthobrick(-2,-1,-2;2,1,0) and not ball -maxh=0.5;\n'];
0447 elec_pos = linspace( -1.5,1.5,5)';
0448 elec_pos = [ elec_pos, elec_pos*[0,0,0,0], elec_pos*0+1];
0449 elec_shape=[0.3];
0450 elec_obj = 'top';
0451 [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0452 img = mk_image( fmdl, 1);
0453 img.elem_data(mat_idx{2}) = 1.1;
0454
0455 fmdl = img;
0456
0457 case 10;
0458 shape_str = ['solid top = plane(0,0,0;0,0,1);\n' ...
0459 'solid mainobj= top and orthobrick(-3,-3,-2;3,3,0) -maxh=0.5;\n'];
0460 [elec_pos_x,elec_pos_y] = meshgrid(linspace( -1.5,1.5,5),linspace(-2,2,7));
0461 elec_pos = [ elec_pos_x(:), elec_pos_y(:), ones(size(elec_pos_x(:)))*[0,0,0,1] ];
0462 elec_shape=[0.2];
0463 elec_obj = 'top';
0464 [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0465
0466 case 11;
0467 shape_str = ['solid cyl = cylinder (0,0,0; 0,0,1; 1.0); \n', ...
0468 'solid tank = orthobrick(-2,-2,0;2,2,0.4) and cyl; \n', ...
0469 'solid fish = ellipsoid(0.2,0.2,0.2;0.2,0,0;0,0.1,0;0,0,0.1); tlo fish;\n', ...
0470 'solid mainobj= tank and not fish -maxh=0.3;\n'];
0471 n_elec = 7;
0472 th = linspace(0,2*pi,n_elec+1)'; th(end) = [];
0473 cs = [cos(th), sin(th)];
0474 elec_pos = [ cs, 0.2+0*th, cs, 0*th];
0475 elec_shape=[0.05];
0476 for i=1:n_elec; elec_obj{i} = 'cyl'; end
0477 i=i+1;elec_pos(i,:) = [ 0 ,0.2,0.2,-1,0,0]; elec_obj{i} = 'fish';
0478 i=i+1;elec_pos(i,:) = [ 0.4,0.2,0.2, 1,0,0]; elec_obj{i} = 'fish';
0479 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0480
0481 case 12;
0482 shape_str = ['solid top = ellipsoid(0,0,0; 0,0,1; 1,0,0; 0,1,0); \n' ...
0483 'solid mainobj= top and orthobrick(-2,-2,0;2,2,2) -maxh=0.1;\n'];
0484 deg2rad = pi/180;
0485 th = (-70:20:70)'*deg2rad;
0486 elec_pos = [0*th,sin(th),cos(th),0*th,sin(th),cos(th); ...
0487 sin(th),0*th,cos(th),sin(th),0*th,cos(th)];
0488 elec_shape=[0.05];
0489 elec_obj = 'top';
0490 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0491
0492 case 13;
0493 shape_str = ['solid cyl = cylinder (0,0,0; 0,1,0; 1); \n', ...
0494 'solid bottom = plane(0, 0,0;0,-1,0);\n' ...
0495 'solid top = plane(0,10,0;0, 1,0);\n' ...
0496 'solid cut1 = plane(0, 4,0;0,-1,0);\n' ...
0497 'solid cut2 = plane(0, 6,0;0, 1,0);\n' ...
0498 'solid ball = cyl and cut1 and cut2; tlo ball;\n' ...
0499 'solid mainobj= ( top and (not cut2) and cyl ) or ' ...
0500 '(bottom and (not cut1) and cyl ) -maxh=0.8;\n'];
0501 elec_pos = [ 0, 10, 0, 0, 1, 0;
0502 0, 0, 0, 0, -1, 0];
0503 elec_shape=[1.0];
0504 elec_obj = {'top','bottom'};
0505 [fmdl,mat_idx] = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0506 fmdl = mk_image(fmdl,1);
0507 fmdl.elem_data(mat_idx{2}) = 1.1;
0508
0509 case 14;
0510 shape_str = [ ...
0511 'solid cyl = cylinder (0,0,0; 0,0,1; 1); \n', ...
0512 'solid mainobj= plane(0,0,0;0,0,-1)\n' ...
0513 'and plane(0,0,2;0,0,1)\n' ...
0514 'and cyl -maxh=0.3;\n' ...
0515 'solid in_elec = sphere(0,-1,1;0.2)' ...
0516 'and not sphere(0,-1,1;0.15) -maxh=0.05;\n' ...
0517 'solid in_elec0= in_elec and mainobj;\n' ...
0518 'tlo in_elec0 cyl;\n' ...
0519 'solid out_elec = sphere(0,-1,1;0.4)' ...
0520 'and not sphere(0,-1,1;0.35) -maxh=0.05;\n' ...
0521 'solid out_elec0= out_elec and mainobj;\n' ...
0522 'tlo out_elec0 cyl;\n'];
0523
0524 elec_pos = [ 0, -1, 0, NaN,NaN,NaN;
0525 1, 0, 1, 1, 0, 0;
0526 0, 1, 1.2, 0, 1, 0;
0527 0, -1, 1.2, NaN,NaN,NaN;
0528 0, -1, 1.4, NaN,NaN,NaN];
0529 elec_shape=[0.1];
0530 elec_obj = 'cyl';
0531 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0532
0533 fmdl.electrode = fmdl.electrode(2:end);
0534 otherwise;
0535 error('huh?')
0536 end