0001 function [fmdl,mat_idx] = ng_mk_extruded_model(shape, elec_pos, elec_shape, ...
0002 extra_ng_code)
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
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067 if ischar(shape) && strcmp(shape,'UNIT_TEST'); fmdl = do_unit_test; return; end
0068
0069 citeme(mfilename);
0070
0071 if nargin < 4; extra_ng_code = {'',''}; end
0072 copt.cache_obj = { shape, elec_pos, elec_shape, extra_ng_code};
0073 copt.fstr = 'ng_mk_extruded_models';
0074 args = { shape, elec_pos, elec_shape, extra_ng_code};
0075
0076 fmdl = eidors_cache(@mk_extruded_model, args, copt);
0077
0078 mat_idx = fmdl.mat_idx;
0079
0080 function fmdl = mk_extruded_model(shape, elec_pos, elec_shape, ...
0081 extra_ng_code)
0082
0083 fnstem = tempname;
0084 geofn= [fnstem,'.geo'];
0085 meshfn= [fnstem,'.vol'];
0086
0087 [tank_height, tank_shape, tank_maxh, is2D] = parse_shape(shape);
0088 [elecs, centres] = parse_elecs(elec_pos, elec_shape, tank_shape, tank_height, is2D );
0089 write_geo_file(geofn, tank_height, tank_shape, ...
0090 tank_maxh, elecs, extra_ng_code);
0091
0092 if 0
0093 plot(tank_shape.vertices(:,1),tank_shape.vertices(:,2));
0094 hold on
0095 plot(centres(:,1),centres(:,2),'sk')
0096 for i = 1:size(elecs,2)
0097 dirn = elecs(i).normal;
0098 quiver(centres(i,1),centres(i,2),dirn(1),dirn(2),'k');
0099 end
0100 hold off
0101 axis equal
0102 end
0103
0104 call_netgen( geofn, meshfn );
0105
0106 fmdl = ng_mk_fwd_model( meshfn, centres, 'ng', [],0.01,...
0107 @ng_remove_electrodes);
0108
0109 delete(geofn); delete(meshfn);
0110
0111 if is2D
0112 fmdl = mdl2d_from3d(fmdl);
0113 end
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131 function [tank_height, tank_shape, tank_maxh, is2D] = parse_shape(shape)
0132
0133
0134
0135 is2D = false;
0136 tank_maxh = 0;
0137 tank_shape = [];
0138 tank_shape.curve_type = 1;
0139 curve_info = [];
0140
0141 if iscell(shape) && length(shape)>2
0142 tank_height = shape{1};
0143 if ~iscell(shape{2})
0144 points = shape{2};
0145 else
0146 c = shape{2};
0147 points = c{1};
0148 if numel(shape{2}) > 1
0149 tank_shape.additional_shapes = c(2:end);
0150 end
0151 end
0152
0153 if ~iscell(shape{3})
0154 tank_shape.curve_type = shape{3};
0155 if iscell(tank_shape.curve_type)
0156 tank_shape.curve_type = tank_shape.curve_type{1};
0157 end
0158 else
0159 c = shape{3};
0160 tank_shape.curve_type = c{1};
0161 if numel(shape{3}) > 1
0162 tank_shape.additional_curve_type = c(2:end);
0163 end
0164 end
0165
0166 if max(size(tank_shape.curve_type)) > 1
0167 curve_info = tank_shape.curve_type;
0168 tank_shape.curve_type = curve_info(1);
0169 end
0170
0171
0172
0173 if length(shape) > 3
0174 tank_maxh = shape{4};
0175 end
0176 else
0177 points = shape;
0178 end
0179
0180 spln_sgmnts = zeros(size(points));
0181 if tank_shape.curve_type == 2
0182 [points, spln_sgmnts] = remove_linear_control_points(points);
0183 end
0184
0185 if ~isempty(curve_info)
0186 N = curve_info(2);
0187 else
0188 N = 50;
0189 end
0190 points = interpolate(points,N, tank_shape.curve_type);
0191 spln_sgmnts = zeros(size(points));
0192
0193 if isfield(tank_shape, 'additional_curve_type')
0194 for i = 1:numel(tank_shape.additional_curve_type)
0195 if numel(tank_shape.additional_curve_type{i}) == 1
0196 N = 50;
0197 else
0198 N = tank_shape.additional_curve_type{i}(2);
0199 end
0200 tank_shape.additional_shapes{i} = interpolate(...
0201 tank_shape.additional_shapes{i},N, tank_shape.additional_curve_type{i}(1));
0202 end
0203 end
0204
0205
0206 if tank_shape.curve_type == 3
0207 if ~isempty(curve_info)
0208 n_samples = curve_info(2);
0209 else
0210 n_samples = 50;
0211 end
0212 points = interpolate_shape(points, n_samples);
0213 spln_sgmnts = zeros(size(points));
0214 end
0215
0216
0217 if tank_shape.curve_type == 4
0218 if ~isempty(curve_info)
0219 n_samples = curve_info(2);
0220 else
0221 n_samples = 50;
0222 end
0223 points = fourier_interpolate_shape(points, n_samples);
0224 spln_sgmnts = zeros(size(points));
0225 end
0226
0227 tank_shape.centroid = calc_centroid(points);
0228 tank_shape.spln_sgmnts = spln_sgmnts;
0229
0230 tank_shape.vertices = points;
0231
0232 range_points = max(points) - min(points);
0233 tank_shape.size = 0.5 * sqrt(sum(range_points.^2));
0234
0235 if tank_height==0
0236 is2D = 1;
0237
0238
0239 tank_height = tank_shape.size/5;
0240 if tank_maxh>0
0241 tank_height = min(tank_height,2*tank_maxh);
0242 end
0243 end
0244
0245
0246 tank_shape.edge_normals = [];
0247 tank_shape.vertex_dir = [];
0248
0249 tmp = points;
0250 tmp(end+1,:) = tmp(1,:);
0251
0252 edges = diff(tmp,1);
0253 tmp = [];
0254
0255
0256 tmp = circshift(edges, [0 1]);
0257
0258 lngth = sqrt(sum(tmp.^2, 2));
0259 tmp(:,1) = -tmp(:,1) ./ lngth;
0260 tmp(:,2) = tmp(:,2) ./ lngth;
0261 tank_shape.edge_normals = tmp;
0262
0263 tank_shape.vertex_dir = calc_vertex_dir(points, edges, ...
0264 tank_shape.edge_normals);
0265
0266
0267 tmp = [];
0268 polar = zeros(size(points));
0269 for i = 1:length(points)
0270 tmp = points(i,:) - tank_shape.centroid;
0271 [polar(i,1) polar(i,2)] = cart2pol(tmp(1),tmp(2));
0272 end
0273 tank_shape.vertices_polar = polar;
0274
0275 tank_shape.convex = calc_convex(tank_shape.vertices);
0276
0277
0278 if 0
0279 pts = edges./2 + points;
0280 plot(tank_shape.vertices(:,1),tank_shape.vertices(:,2),'-o'); hold on;
0281 plot(tank_shape.centroid(:,1),tank_shape.centroid(:,2),'+');
0282 plot(tank_shape.vertices(:,1)+0.05*tank_shape.vertex_dir(:,1),...
0283 tank_shape.vertices(:,2)+0.05*tank_shape.vertex_dir(:,2),'ro-')
0284 quiver(pts(:,1),pts(:,2),tank_shape.edge_normals(:,1),tank_shape.edge_normals(:,2));
0285 hold off
0286 end
0287
0288
0289 function new_points = interpolate(points, N, curve_type)
0290 switch curve_type
0291 case 3
0292
0293 new_points = interpolate_shape(points, N);
0294 case 4
0295
0296 new_points = fourier_interpolate_shape(points, N);
0297 otherwise
0298
0299 new_points = points;
0300 end
0301
0302
0303
0304
0305
0306
0307
0308
0309 function [points, spln_sgmnts] = remove_linear_control_points(points)
0310 n_points = length(points);
0311 points(end+1,:) = points(1,:);
0312 spln_sgmnts(1:(n_points/2)) = 1;
0313 for i = 1:2:n_points
0314 a = (points(i+1,:) - points(i,:));
0315 a = a/norm(a);
0316 b = (points(i+2,:) - points(i,:));
0317 b = b/norm(b);
0318 if a(1) == b(1) && a(2) == b(2)
0319 spln_sgmnts(i/2 + 0.5) = 0;
0320 end
0321 end
0322 idx = find(spln_sgmnts==0) * 2;
0323 points(idx,:) = [];
0324 points(end,:) = [];
0325
0326
0327
0328
0329
0330
0331
0332 function out = interpolate_shape(points, n_points)
0333
0334
0335
0336 [pp m] = piece_poly_fit(points);
0337 p = linspace(0,1,n_points+1)'; p(end) = [];
0338 [th xy] = piece_poly_fit(pp,0,p);
0339 tmp = [th xy];
0340 tmp = sortrows(tmp,-1);
0341 xy = tmp(:,2:3);
0342
0343 out = xy + repmat(m, [n_points,1]);
0344
0345
0346
0347
0348
0349
0350 function out = fourier_interpolate_shape(points, n_points)
0351
0352
0353
0354 pp = fourier_fit(points, size(points,1)-1);
0355 p = linspace(0,1,n_points+1)'; p(end) = [];
0356 xy = fourier_fit(pp,p);
0357
0358
0359
0360
0361
0362 out = xy;
0363
0364
0365 function out = calc_vertex_dir(points, edges, edgnrm)
0366
0367
0368
0369
0370 edg = [edges(end,:) ; edges];
0371 edgnrm = [edgnrm(end,:) ; edgnrm];
0372
0373 out = zeros(size(points));
0374 for i = 1:length(points)
0375 p1 = points(i,:) + edgnrm(i,:);
0376 p2 = points(i,:) + edgnrm(i+1,:);
0377
0378 dir1(1) = edgnrm(i,2); dir1(2) = -edgnrm(i,1);
0379 dir2(1) = edgnrm(i+1,2); dir2(2) = -edgnrm(i+1,1);
0380
0381
0382 if isempty(find(abs(dir1 - dir2) > 1e-14))
0383 out(i,:) = edgnrm(i,:);
0384 else
0385 A = [dir1' , -dir2'];
0386 u = (p2 - p1)';
0387 x = A\u;
0388 out(i,:) = x(1) * dir1 + p1 - points(i,:);
0389 end
0390 end
0391
0392 function out = calc_centroid(points)
0393
0394
0395
0396
0397
0398
0399
0400
0401 n_points = size(points,1);
0402 if n_points == 3
0403 out = mean(points);
0404 return
0405 end
0406
0407 out = 0;
0408 pts = [points ; points(1,:)];
0409
0410
0411 m = mean(points);
0412
0413 if ~inpolygon(m(1),m(2),points(:,1),points(:,2))
0414 f1 = figure;
0415 set(f1,'Name', 'Select a point within the shape');
0416 plot(pts(:,1),pts(:,2));
0417 m = ginput(1);
0418 close(f1)
0419 end
0420
0421 tmp = 0;
0422 tot_area = 0;
0423 for i = 1:n_points
0424 a = pts(i,:);
0425 b = pts(i+1,:);
0426 cntrd = (m + a + b)/3;
0427 area = 0.5 * abs(det([m 1; a 1; b 1]));
0428 tmp = tmp + cntrd*area;
0429 tot_area = tot_area + area;
0430 end
0431
0432 out = tmp./tot_area;
0433
0434 function out = calc_convex(verts)
0435
0436
0437
0438
0439
0440 n_verts = size(verts,1);
0441 tmp = [verts(end,:); verts; verts(1,:)];
0442 verts = tmp;
0443
0444 for i = 2:n_verts+1
0445 v1 = [verts(i-1,:) - verts(i,:), 0];
0446 v2 = [verts(i+1,:) - verts(i,:), 0];
0447 cp = cross(v1,v2);
0448 out(i-1) = cp(3) >= 0;
0449 end
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468 function [elecs, centres] = parse_elecs(elec_pos, elec_shape, tank_shape, hig, is2D )
0469
0470 if isempty(elec_pos)
0471 elecs = [];
0472 centres = [];
0473 return;
0474 end
0475
0476 if is2D
0477 elec_pos(:,3) = hig/2;
0478 end
0479
0480
0481 rad = tank_shape.size;
0482
0483
0484
0485 if size(elec_pos,1) == 1
0486
0487 n_elecs= elec_pos(1);
0488 offset = elec_pos(2) - floor(elec_pos(2));
0489 switch floor(elec_pos(2))
0490 case 0
0491 th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0492 th = th + offset*2*pi;
0493 ind = th >= 2*pi;
0494 th(ind) = th(ind) - 2*pi;
0495 case 1
0496
0497 if 1
0498 pp = fourier_fit(tank_shape.vertices,...
0499 size(tank_shape.vertices,1) - 1,tank_shape.vertices(1,:));
0500 p = linspace(0,1,n_elecs+1)'; p(end) = [];
0501 p = p + offset;
0502 xy = fourier_fit(pp,p);
0503
0504
0505 th = atan2(xy(:,2) - tank_shape.centroid(2), ...
0506 xy(:,1) - tank_shape.centroid(1));
0507
0508 elseif any( tank_shape.curve_type == [1,2,3] )
0509
0510
0511 pp= piece_poly_fit(tank_shape.vertices);
0512 p = linspace(1,0,n_elecs+1)'; p(end) = [];
0513 off = offset*2*pi + tank_shape.vertices_polar(1,1);
0514 th = piece_poly_fit(pp,off,p);
0515 else
0516 error('curve_type unrecognized');
0517 end
0518 end
0519
0520 on_elecs = ones(n_elecs, 1);
0521 el_th = [];
0522 el_z = [];
0523 for i=3:length(elec_pos)
0524 el_th = [el_th; th];
0525 el_z = [el_z ; on_elecs*elec_pos(i)];
0526 end
0527 else
0528 el_th = elec_pos(:,1)*2*pi/360;
0529 el_z = elec_pos(:,2);
0530 end
0531
0532 n_elecs= size(el_z,1);
0533
0534 if size(elec_shape,1) == 1
0535 elec_shape = ones(n_elecs,1) * elec_shape;
0536 end
0537
0538 for i= 1:n_elecs
0539 row = elec_shape(i,:);
0540 elecs(i) = elec_spec( row, is2D, hig, rad );
0541 end
0542
0543
0544
0545 for i= 1:n_elecs;
0546
0547
0548 [centres(i,1:2), normal] = calc_elec_centre(tank_shape, el_th(i));
0549
0550
0551
0552
0553
0554 centres(i,3) = el_z(i);
0555 elecs(i).pos = centres(i,:);
0556 if elecs(i).discretize > 0
0557
0558
0559 th = cart2pol(normal(1),normal(2));
0560 frac = 2*pi /elecs(i).discretize ;
0561 th = frac * round( th / frac);
0562 [normal(1) normal(2)] = pol2cart(th,1);
0563 end
0564 elecs(i).normal = normal;
0565
0566 end
0567
0568 if n_elecs == 0
0569 elecs= struct([]);
0570 centres= [];
0571 end
0572
0573
0574
0575 function [pos, normal] = calc_elec_centre(tank_shape, th)
0576
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586 if th > pi; th = th - 2*pi; end
0587
0588
0589 vert_pol = tank_shape.vertices_polar;
0590
0591
0592 n_vert = size(vert_pol,1);
0593 vert_pol = [vert_pol , (1:n_vert)'];
0594
0595 vert_pol = sortrows(vert_pol,1);
0596
0597
0598 idx = find(vert_pol(:,1) > th, 1, 'first');
0599 if isempty(idx); idx = 1; end
0600 edg_no = vert_pol(idx,3);
0601
0602
0603 normal = tank_shape.edge_normals(edg_no,:);
0604
0605 v1 = edg_no;
0606 if edg_no == n_vert
0607 v2 = 1;
0608 else
0609 v2 = v1+1;
0610 end
0611 vert_pol = [];
0612
0613
0614 vert_pol = tank_shape.vertices_polar;
0615 vert = tank_shape.vertices;
0616 cntr = tank_shape.centroid;
0617
0618 AB = sqrt(sum( (vert(v1,:) - cntr).^2 ));
0619 AC = sqrt(sum( (vert(v2,:) - cntr).^2 ));
0620 DAB = abs(vert_pol(v1,1)-th);
0621 if DAB > pi, DAB = abs( DAB - 2*pi); end;
0622 DAC = abs(vert_pol(v2,1)-th);
0623 if DAC > pi, DAC = abs( DAC - 2*pi); end;
0624 if DAC ~= 0
0625 ratio = AB * sin(DAB) / (AC * sin(DAC));
0626 pos = vert(v1,:) + ( ratio / (1 + ratio) ) * (vert(v2,:) - vert(v1,:));
0627 else
0628 pos = vert(v2,:);
0629 end
0630
0631
0632
0633 function [pos, normal] = calc_elec_centre_spline(tank_shape, th)
0634
0635
0636
0637
0638
0639 if th > pi; th = th - 2*pi; end
0640
0641 vert_pol = tank_shape.vertices_polar;
0642
0643
0644 if mod(size(vert_pol,1),2)
0645 error(['The number of points must be even. '...
0646 'One de Boor control point for every vertex']);
0647 end
0648
0649
0650
0651 if tank_shape.curve_type == 2 || tank_shape.curve_type == 3
0652 vert_pol(2:2:end,:) = [];
0653 end
0654
0655 n_vert = size(vert_pol,1);
0656
0657 vert_pol = [vert_pol , (1:n_vert)'];
0658
0659 vert_pol = sortrows(vert_pol,1);
0660
0661
0662 idx = find(vert_pol(:,1) > th, 1, 'first');
0663 if isempty(idx); idx = 1; end
0664 edg_no = vert_pol(idx,3);
0665
0666 v1 = edg_no;
0667 if edg_no == n_vert
0668 v2 = 1;
0669 else
0670 v2 = v1+1;
0671 end
0672 vert_pol = [];
0673
0674
0675 v1 = 2 * v1 - 1;
0676 v2 = 2 * v2 - 1;
0677
0678
0679
0680 C = tank_shape.centroid;
0681 P0 = tank_shape.vertices(v1,:) - C;
0682 P1 = tank_shape.vertices(v1+1,:) - C;
0683 P2 = tank_shape.vertices(v2,:) - C;
0684
0685
0686
0687 [x, y] = pol2cart(th, 1);
0688
0689 g = y/x;
0690
0691
0692
0693
0694
0695 f = @(t) (P2 - 2*P1 + P0)*t^2 + 2*(P1 - P0)*t + P0;
0696
0697 df = @(t) 2*(P2 - 2*P1 + P0)*t + 2*(P1 - P0);
0698
0699
0700 p0 = P0(2) - g * P0(1);
0701 p1 = P1(2) - g * P1(1);
0702 p2 = P2(2) - g * P2(1);
0703
0704
0705 a = (p2 - 2*p1 + p0);
0706 b = 2* (p1 - p0);
0707 c = p0;
0708
0709 if abs(a) < 1e-10
0710 t = -c/b;
0711 pos = f(t) + C;
0712 tmp = df(t);
0713 normal = [-tmp(2), tmp(1)] / sqrt(sum(tmp.^2));
0714 return;
0715 end
0716
0717
0718 D = b^2 - 4*a*c;
0719
0720
0721 if D == 0
0722 t = -b / (2 * a);
0723
0724 elseif D > 0
0725 t1 = (-b - sqrt(D) ) / (2 * a);
0726 t2 = (-b + sqrt(D) ) / (2 * a);
0727 if t1 >= 0 && t1 <= 1
0728 t = t1;
0729 else
0730 t = t2;
0731 end
0732 else
0733 error('Something went wrong, cannot place electrode on spline');
0734 end
0735
0736 pos = f(t) + C;
0737 tmp = df(t);
0738 normal = [-tmp(2), tmp(1)]/ sqrt(sum(tmp.^2));
0739
0740
0741
0742
0743 function elec = elec_spec( row, is2D, hig, rad )
0744 if is2D
0745 if length(row)>=2 && row(2) == -1
0746
0747 elec.shape = 'P' ;
0748 if length(row)>=3 && row(3) > 0
0749 elec.dims = row(3);
0750 else
0751 elec.dims = rad;
0752 end
0753 else
0754
0755
0756
0757 elec.shape = 'R';
0758 elec.dims = [row(1),hig];
0759 end
0760 else
0761 if length(row)<2 || row(2) == 0
0762 elec.shape = 'C';
0763 elec.dims = row(1);
0764 elseif row(2) == -1
0765
0766 elec.shape = 'P';
0767 if length(row)>=3 && row(3) > 0
0768 elec.dims = row(3);
0769 else
0770 elec.dims = rad;
0771 end
0772 elseif row(2)>0
0773 elec.shape = 'R';
0774 elec.dims = row(1:2);
0775 else
0776 error('negative electrode width');
0777 end
0778 end
0779
0780 if length(row)>=3 && row(3) > 0
0781 elec.maxh = sprintf('-maxh=%f', row(3));
0782 else
0783 elec.maxh = '';
0784 end
0785
0786 if length(row)<4 || row(4) == 0
0787 elec.model = 'cem';
0788 else
0789 elec.model = 'pem';
0790 end
0791
0792
0793 if length(row) < 5 || row(5) == 0
0794 elec.discretize = 0;
0795 else
0796 elec.discretize = row(5);
0797 end
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811 function write_geo_file(geofn, tank_height, tank_shape, ...
0812 tank_maxh, elecs, extra_ng_code)
0813 fid=fopen(geofn,'w');
0814 write_header(fid,tank_height,tank_shape,tank_maxh,extra_ng_code);
0815
0816 n_verts = size(tank_shape.vertices,1);
0817 n_elecs = length(elecs);
0818
0819
0820
0821
0822
0823 pts_elecs_idx = [];
0824
0825 for i=1:n_elecs
0826 name = sprintf('elec%04d',i);
0827 pos = elecs(i).pos;
0828 dirn = elecs(i).normal;
0829 switch elecs(i).shape
0830 case 'C'
0831 write_circ_elec(fid,name, pos, dirn, ...
0832 elecs(i).dims, tank_shape.centroid, elecs(i).maxh);
0833 case 'R'
0834 write_rect_elec(fid,name, pos, dirn, ...
0835 elecs(i).dims, tank_shape.size/10, elecs(i).maxh);
0836
0837
0838
0839
0840 otherwise; error('unknown electrode shape');
0841 end
0842
0843 end
0844 fprintf(fid,'solid trunk = bound');
0845 if isfield(tank_shape,'additional_shapes')
0846 for i = 1:length(tank_shape.additional_shapes)
0847 fprintf(fid,' and not add_obj%04d',i);
0848 end
0849 end
0850 fprintf(fid,';\n');
0851
0852 if isfield(tank_shape,'additional_shapes')
0853 for i = 1:length(tank_shape.additional_shapes)
0854 fprintf(fid,'solid add_obj%04dc = add_obj%04d',i,i);
0855 for j = (i+1):length(tank_shape.additional_shapes)
0856 fprintf(fid,' and not add_obj%04d',j);
0857 end
0858
0859
0860
0861
0862
0863
0864
0865 fprintf(fid,[' and plane(0,0,0;0,0,-1)\n' ...
0866 ' and plane(0,0,%6.2f;0,0,1)'],tank_height);
0867 fprintf(fid,';\n');
0868 end
0869 end
0870
0871 if tank_maxh ~= 0
0872 fprintf(fid,'tlo trunk -transparent -maxh=%f;\n',tank_maxh);
0873 else
0874 fprintf(fid,'tlo trunk -transparent;\n');
0875 end
0876 if ~isempty(extra_ng_code{1})
0877 fprintf(fid,'tlo %s -col=[0,1,0];\n',extra_ng_code{1});
0878 end
0879
0880 if isfield(tank_shape,'additional_shapes')
0881 for i = 1:length(tank_shape.additional_shapes)
0882 fprintf(fid,'tlo add_obj%04dc -col=[0,1,0];\n',i);
0883 end
0884 end
0885
0886 for i=1:n_elecs
0887 if any(i == pts_elecs_idx); continue; end
0888 fprintf(fid,'tlo elec%04d -col=[1,0,0] %s;\n',i,elecs(i).maxh);
0889 end
0890
0891
0892 fclose(fid);
0893
0894
0895
0896 function write_header(fid,tank_height,tank_shape,maxsz,extra)
0897 if maxsz==0;
0898 maxsz = '';
0899 else
0900 maxsz = sprintf('-maxh=%f',maxsz);
0901 end
0902
0903 if ~isempty( extra{1} )
0904 extra{1} = [' and not ',extra{1}];
0905 end
0906
0907
0908 fprintf(fid,'#Automatically generated by ng_mk_extruded_model\n');
0909 fprintf(fid,'algebraic3d\n');
0910 fprintf(fid,'%s\n',extra{2});
0911
0912 fprintf(fid,'curve3d extrsncurve=(2; 0,0,0; 0,0,%6.2f; 1; 2,1,2);\n', ...
0913 tank_height+1);
0914
0915
0916 write_curve(fid,tank_shape,'outer', 1.15);
0917 write_curve(fid,tank_shape,'inner', 0.99);
0918 write_curve(fid,tank_shape,'surf', 1);
0919
0920 fprintf(fid,['solid bound= plane(0,0,0;0,0,-1)\n' ...
0921 ' and plane(0,0,%6.2f;0,0,1)\n' ...
0922 ' and extrusion(extrsncurve;surf;0,1,0)'...
0923 '%s %s;\n'],tank_height,extra{1},maxsz);
0924
0925 fprintf(fid,['solid inner_bound= plane(0,0,0;0,0,-1)\n' ...
0926 ' and plane(0,0,%6.2f;0,0,1)\n' ...
0927 ' and extrusion(extrsncurve;inner;0,1,0)'...
0928 '%s %s;\n'],tank_height,extra{1},maxsz);
0929
0930 fprintf(fid,['solid outer_bound= plane(0,0,0;0,0,-1)\n' ...
0931 ' and plane(0,0,%6.2f;0,0,1)\n' ...
0932 ' and extrusion(extrsncurve;outer;0,1,0)'...
0933 '%s %s;\n'],tank_height,extra{1},maxsz);
0934
0935
0936 if ~isfield(tank_shape, 'additional_shapes'), return, end
0937
0938 for i = 1:length(tank_shape.additional_shapes)
0939 name_curve = sprintf('add_curve%04d',i);
0940 write_curve(fid,tank_shape.additional_shapes{i},name_curve);
0941 name_obj = sprintf('add_obj%04d',i);
0942 fprintf(fid,['solid %s= plane(0,0,%6.2f;0,0,-1)\n' ...
0943 ' and plane(0,0,%6.2f;0,0,1)\n' ...
0944 ' and extrusion(extrsncurve;%s;0,1,0)'...
0945 '%s %s;\n'],name_obj,-i,tank_height+i,name_curve,extra{1},maxsz);
0946 end
0947
0948
0949 function write_curve(fid, tank_shape, name, scale)
0950 if nargin <4
0951 scale = 1;
0952 end
0953
0954 is_struct = isstruct(tank_shape);
0955 if ~is_struct
0956 vertices = tank_shape;
0957 STRUCT = false;
0958 if scale ~= 1
0959 warning('Scale is ignored when second input is an array');
0960 scale = 1;
0961 end
0962 elseif scale ~= 1
0963 vertices = tank_shape.vertices + ...
0964 (scale-1)*tank_shape.vertex_dir*tank_shape.size;
0965 else
0966 vertices = tank_shape.vertices;
0967 end
0968 n_vert = size(vertices,1);
0969
0970 fprintf(fid,'curve2d %s=(%d; \n', name, n_vert);
0971
0972 for i = 1:n_vert
0973
0974
0975
0976
0977
0978 fprintf(fid,' %6.4f, %6.4f;\n',[-1 1].*vertices(n_vert-i+1,:));
0979
0980 end
0981 if is_struct
0982 spln_sgmnts = tank_shape.spln_sgmnts;
0983 else
0984 spln_sgmnts = zeros(max(size(vertices)));
0985 end
0986 n_sgmnts = length(spln_sgmnts);
0987 fprintf(fid,' %d;\n',n_sgmnts);
0988 cv = 1;
0989 for i = 1:n_sgmnts
0990 if spln_sgmnts(i)
0991 if i == n_sgmnts
0992 fprintf(fid,' %d, %d, %d, %d );\n\n\n', 3, cv,cv+1, 1);
0993 else
0994 fprintf(fid,' %d, %d, %d, %d; \n', 3, cv, cv+1, cv+2);
0995 end
0996 cv = cv + 2;
0997 else
0998 if i == n_sgmnts
0999 fprintf(fid,' %d, %d, %d );\n\n\n', 2, cv, 1);
1000 else
1001 fprintf(fid,' %d, %d, %d; \n', 2, cv, cv+1);
1002 end
1003 cv = cv + 1;
1004 end
1005 end
1006
1007
1008 function write_circ_elec(fid,name,c, dirn, rd, centroid, maxh)
1009
1010
1011
1012
1013
1014
1015 dirn(3) = 0; dirn = dirn/norm(dirn);
1016
1017 fprintf(fid,'solid %s = ', name);
1018 fprintf(fid,[' outer_bound and not inner_bound and '...
1019 'cylinder(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f;%6.3f) '...
1020 'and plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) '...
1021 'and not bound;\n'], ...
1022 c(1)-dirn(1),c(2)-dirn(2),c(3)-dirn(3),...
1023 c(1)+dirn(1),c(2)+dirn(2),c(3)+dirn(3), rd, ...
1024 centroid(1), centroid(2), 0, -dirn(1), -dirn(2), dirn(3));
1025
1026 function write_rect_elec(fid,name,c, dirn,wh,d,maxh)
1027
1028
1029
1030
1031 w = wh(1); h= wh(2);
1032 dirn(3) = 0; dirn = dirn/norm(dirn);
1033 dirnp = [-dirn(2),dirn(1),0];
1034 dirnp = dirnp/norm(dirnp);
1035
1036 bl = c - (d/2)* dirn + (w/2)*dirnp - [0,0,h/2];
1037 tr = c + (d/2)* dirn - (w/2)*dirnp + [0,0,h/2];
1038 fprintf(fid,'solid %s = outer_bound and not inner_bound and', name);
1039 fprintf(fid,' plane (%6.3f,%6.3f,%6.3f;0, 0, -1) and\n', ...
1040 bl(1),bl(2),bl(3));
1041 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1042 bl(1),bl(2),bl(3),-dirn(1),-dirn(2),0);
1043 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1044 bl(1),bl(2),bl(3),dirnp(1),dirnp(2),0);
1045 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;0, 0, 1) and\n', ...
1046 tr(1),tr(2),tr(3));
1047 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f) and\n', ...
1048 tr(1),tr(2),tr(3),dirn(1),dirn(2),0);
1049 fprintf(fid,' plane(%6.3f,%6.3f,%6.3f;%6.3f,%6.3f,%6.3f )\n', ...
1050 tr(1),tr(2),tr(3),-dirnp(1),-dirnp(2),0);
1051 fprintf(fid,' and not bound;\n');
1052
1053
1054
1055
1056
1057
1058
1059
1060 function [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_remove_electrodes...
1061 (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1062
1063 fc = [];
1064
1065
1066 N_obj = max(mat_ind);
1067
1068
1069 elec_ind = mat_ind > (N_obj - N_elec);
1070
1071 in = unique(simp(~elec_ind,:));
1072 out = unique(simp(elec_ind,:));
1073 boundary = intersect(in,out);
1074 out = setdiff(out,boundary);
1075
1076
1077 remove_simp = any( ismember(simp,out), 2);
1078 simp0 = simp;
1079 simp( remove_simp,:) = [];
1080
1081
1082 vtx_renum = logical( zeros(size(vtx,1),1) );
1083 vtx_renum( in ) = logical(1);
1084 vtx_renum = cumsum(vtx_renum);
1085
1086 vtx(out,:) = [];
1087 simp = reshape( vtx_renum(simp), size(simp));
1088
1089
1090
1091 srf= double( find_boundary(simp) );
1092 bc = ones(size(srf,1),1);
1093
1094
1095 for i=1:N_elec;
1096 eleci_obj = mat_ind == (N_obj - N_elec + i);
1097 this_elec = unique( simp0( eleci_obj, : ));
1098 eleci_nodes = vtx_renum( intersect( this_elec, in ));
1099
1100
1101
1102
1103
1104 eleci_srf = all( ismember(srf, eleci_nodes), 2);
1105 bc( eleci_srf ) = i+1;
1106 end
1107
1108 mat_ind( remove_simp) = [];
1109
1110
1111
1112
1113
1114
1115 function [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_remove_electrodes_old...
1116 (srf,vtx,fc,bc,simp,edg,mat_ind, N_elec)
1117
1118
1119 N_obj = max(mat_ind);
1120
1121
1122 e_simp_ind = mat_ind > (N_obj - N_elec);
1123
1124 in = unique(simp(~e_simp_ind,:));
1125 out = unique(simp(e_simp_ind,:));
1126 boundary = intersect(in,out);
1127 out = setdiff(out,boundary);
1128
1129 ext_srf_ind = ismember(srf,out);
1130 ext_srf_ind = ext_srf_ind(:,1) | ext_srf_ind(:,2) | ext_srf_ind(:,3);
1131
1132 srf(ext_srf_ind,:) = [];
1133 bc(ext_srf_ind,:) = [];
1134 fc(ext_srf_ind,:) = [];
1135 simp = simp(~e_simp_ind,:);
1136 mat_ind = mat_ind(~e_simp_ind);
1137
1138
1139 n_unique = numel(unique(bc));
1140 missing = setdiff(1:n_unique, unique(bc));
1141 spare = setdiff(unique(bc), 1:n_unique);
1142
1143 for i = 1:length(missing)
1144 bc( bc==spare(i) ) = missing(i);
1145 end
1146
1147
1148 v = 1:size(vtx,1);
1149 unused_v = setdiff(v, union(unique(simp),unique(srf)));
1150 v(unused_v) = [];
1151 for i = 1:size(vtx,1);
1152
1153
1154 new_v_ind = find(v == i);
1155 simp( simp == i ) = new_v_ind;
1156 srf( srf == i ) = new_v_ind;
1157 end
1158 vtx(unused_v,:) = [];
1159
1160
1161
1162 function [fmdl, mat_idx] = do_unit_test
1163 fmdl = [];
1164 mat_idx = [];
1165 a = [
1166 -0.8981 -0.7492 -0.2146 0.3162 0.7935 0.9615 0.6751 0.0565 -0.3635 -0.9745
1167 0.1404 0.5146 0.3504 0.5069 0.2702 -0.2339 -0.8677 -0.6997 -0.8563 -0.4668 ]';
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191 xx=[
1192 -88.5777 -11.4887 4.6893 49.8963 122.7033 150.3033 195.5103 249.7573 ...
1193 258.8013 279.7393 304.9623 309.2443 322.0923 337.7963 340.6503 348.2633 ...
1194 357.3043 358.7333 361.5873 364.9183 365.3943 366.3453 366.3453 365.3943 ...
1195 362.5393 351.5943 343.5053 326.8513 299.2503 288.3073 264.9923 224.0703 ...
1196 206.4633 162.6833 106.5313 92.2543 57.5153 7.0733 -8.6297 -42.4167 ...
1197 -90.9547 -105.7057 -134.2577 -178.0367 -193.2647 -222.7687 -265.5957 -278.9197 ...
1198 -313.1817 -355.5337 -363.6237 -379.3267 -397.8857 -400.7407 -401.6927 -398.8377 ...
1199 -395.0307 -384.0867 -368.3837 -363.6247 -351.7277 -334.1217 -328.4117 -314.1357 ...
1200 -291.2947 -282.7297 -267.0257 -236.5707 -221.8187 -196.5977 -159.4807 -147.5837];
1201
1202 yy=[
1203 -385.8513 -386.8033 -386.3273 -384.8993 -368.7193 -353.9673 -323.0363 -283.5403 ...
1204 -274.9743 -254.0363 -225.4843 -217.8703 -187.4153 -140.7813 -124.6013 -86.0573 ...
1205 -38.4703 -29.4273 -9.9173 21.0137 32.4347 53.3727 83.8257 93.3437 ...
1206 114.7587 149.0237 161.8717 187.5677 222.3037 231.3447 247.5237 267.5087 ...
1207 271.3177 277.0297 281.3127 279.4097 274.6507 273.2227 276.5547 284.6447 ...
1208 295.1127 297.4927 301.7757 304.1557 302.2537 297.4947 287.5017 282.2667 ...
1209 259.9017 225.6387 213.7427 185.6677 141.4127 125.2337 88.5917 34.8187 ...
1210 17.6897 -22.2803 -73.6723 -85.0923 -117.9263 -163.6083 -176.4573 -205.9613 ...
1211 -245.9343 -256.4023 -275.4373 -304.9403 -315.4083 -332.0623 -352.0473 -355.3783];
1212
1213 a = [xx; yy]';
1214 a = flipud(a);
1215
1216
1217
1218
1219 fmdl = ng_mk_extruded_model({300,a,[4,25]},[16,1.11,150],[1]);
1220 show_fem(fmdl);