ng_mk_2d_model

NG_MG_2D_MODELS create a 2D mesh with Netgen via the in2d interface
 mdl = ng_mk_2d_model(shape)
 mdl = ng_mk_2d_model(shape, elec_pos)
 mdl = ng_mk_2d_model(shape, elec_pos, elec_shape)

 SHAPE can be:
  - xy (Nx2)             : a counter- clockwise list of points in 2D 
                           defining the outer contour
  - {xy, xy1, xy2, ...}  : allows specifying additional counter-clockwise 
                           loops  xy1, xy2, etc, which represent holes in  
                           the bigger contour xy contour
  - {..., maxsz}         : specifies maximum element size of the mesh.
                           If absent, mesh paremeters are controlled by
                           the ng.opt file in the current directory.

 WARNING: Specifying maxsz overwrites the ng.opt in the current directory.

 ELEC_POS (optional) defines electrodes:
  - ep (Nx2)              : a list of points in 2D (will be projected on
                            closest edge of the first contour specified in
                            SHAPE
  - ep (1x1) = N          : the number of equidistant electrodes to create
                            with first electrode on the first point in XY
                            and counter-clockwise ordering. Specify a
                            negative number for clockwise ordering.
  - ep (1x2) = [N offset] : specify offset of the first electrode with
                            respect to the first point of XY
                            (clockwise if negative, counter-clockwise
                            otherwise)
  - {ep1, ep2, ...}       : allows specifying electrodes on the internal 
                            contours specified in SHAPE. Use an empty 
                            array [] if a contour has no electrodes

 ELEC_SHAPE (optional) defines the electrode shape
  - es (1x2) = [wd rfnum] : WD defines width of the electrode (default: 0 
                            i.e. point electrode
                            RFNUM controls amount of local refinement
                            around the electrode.
  - es (NEx2)             : specifies the above for each electrode
                            individually
 
 NOTE: smaller MAXSZ generally requires a lower RFNUM than a coarser mesh
 would.

 Examples:

 xy = [0 0;  1 0; 1 1; 0 1];
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy});
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [0.5 1; 0.5 0; 0 0.5]);
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25]);
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [-5, 0.25]);
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, -0.25]);
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
 mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
 mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1], [4]});
 mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]});
 mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},[0 30]);
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25],[0.2,10;0 20; 0 30; 0 20; 0 10]);
 mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0],[0.2,10;0 20; 0 20; 0 20; 0 20]);
 mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},...
     [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
 th = linspace(2*pi,0,101)';th(1) = [];
 xy = [sin(th),cos(th)];
 mdl = ng_mk_2d_model({xy,0.2});

0001 function mdl = ng_mk_2d_model(varargin)
0002 %NG_MG_2D_MODELS create a 2D mesh with Netgen via the in2d interface
0003 % mdl = ng_mk_2d_model(shape)
0004 % mdl = ng_mk_2d_model(shape, elec_pos)
0005 % mdl = ng_mk_2d_model(shape, elec_pos, elec_shape)
0006 %
0007 % SHAPE can be:
0008 %  - xy (Nx2)             : a counter- clockwise list of points in 2D
0009 %                           defining the outer contour
0010 %  - {xy, xy1, xy2, ...}  : allows specifying additional counter-clockwise
0011 %                           loops  xy1, xy2, etc, which represent holes in
0012 %                           the bigger contour xy contour
0013 %  - {..., maxsz}         : specifies maximum element size of the mesh.
0014 %                           If absent, mesh paremeters are controlled by
0015 %                           the ng.opt file in the current directory.
0016 %
0017 % WARNING: Specifying maxsz overwrites the ng.opt in the current directory.
0018 %
0019 % ELEC_POS (optional) defines electrodes:
0020 %  - ep (Nx2)              : a list of points in 2D (will be projected on
0021 %                            closest edge of the first contour specified in
0022 %                            SHAPE
0023 %  - ep (1x1) = N          : the number of equidistant electrodes to create
0024 %                            with first electrode on the first point in XY
0025 %                            and counter-clockwise ordering. Specify a
0026 %                            negative number for clockwise ordering.
0027 %  - ep (1x2) = [N offset] : specify offset of the first electrode with
0028 %                            respect to the first point of XY
0029 %                            (clockwise if negative, counter-clockwise
0030 %                            otherwise)
0031 %  - {ep1, ep2, ...}       : allows specifying electrodes on the internal
0032 %                            contours specified in SHAPE. Use an empty
0033 %                            array [] if a contour has no electrodes
0034 %
0035 % ELEC_SHAPE (optional) defines the electrode shape
0036 %  - es (1x2) = [wd rfnum] : WD defines width of the electrode (default: 0
0037 %                            i.e. point electrode
0038 %                            RFNUM controls amount of local refinement
0039 %                            around the electrode.
0040 %  - es (NEx2)             : specifies the above for each electrode
0041 %                            individually
0042 %
0043 % NOTE: smaller MAXSZ generally requires a lower RFNUM than a coarser mesh
0044 % would.
0045 %
0046 % Examples:
0047 %
0048 % xy = [0 0;  1 0; 1 1; 0 1];
0049 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy});
0050 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [0.5 1; 0.5 0; 0 0.5]);
0051 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25]);
0052 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [-5, 0.25]);
0053 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, -0.25]);
0054 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0055 % mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0056 % mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1], [4]});
0057 % mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]});
0058 % mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},[0 30]);
0059 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25],[0.2,10;0 20; 0 30; 0 20; 0 10]);
0060 % mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0],[0.2,10;0 20; 0 20; 0 20; 0 20]);
0061 % mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},...
0062 %     [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
0063 % th = linspace(2*pi,0,101)';th(1) = [];
0064 % xy = [sin(th),cos(th)];
0065 % mdl = ng_mk_2d_model({xy,0.2});
0066 
0067 
0068 % (C) 2012-2013 Bartlomiej Grychtol, (C) 2013 Alistair Boyle, License: GPL version 2 or version 3
0069 % $Id: ng_mk_2d_model.m 6926 2024-05-31 15:34:13Z bgrychtol $
0070 
0071 
0072 if ischar(varargin{1}) && strcmp(varargin{1}, 'UNIT_TEST'), mdl = do_unit_test; return, end 
0073 
0074 [shape, elec_pos, elec_shape] = process_input(varargin{:});
0075 
0076 mdl = eidors_cache(@ng_mk_2d_model_do,{shape, elec_pos, elec_shape},'ng_mk_2d_model');
0077 
0078 
0079 
0080 function [shape, elec_pos, elec_shape] = process_input(shape, elec_pos, elec_shape)
0081 
0082 if ~iscell(shape)
0083    shape = {shape};
0084 end
0085 
0086 if nargin < 2
0087     elec_pos = [];
0088 end
0089 if ~iscell(elec_pos)
0090     elec_pos = {elec_pos};
0091 end
0092 
0093 if nargin < 3
0094     elec_shape = [0 10]; % point electrode
0095 end
0096 if size(elec_shape,2) == 1
0097    warning('Refinement factor not specified, using 10');
0098    elec_shape(:,2) = 10;
0099 end
0100 if ~iscell(elec_shape)
0101     elec_shape = {elec_shape};
0102 end
0103 if numel(elec_shape) == 1 && numel(elec_pos) > 1
0104     elec_shape(2:numel(elec_pos)) = elec_shape(1);
0105 end
0106 
0107 
0108 function mdl = ng_mk_2d_model_do(shape, elec_pos, elec_shape)
0109 
0110 [shape,i_wrote_ng_opt] = process_maxsz(shape);
0111 
0112 points = [];
0113 eidx = [];
0114 eref = [];
0115 for i = 1:length(shape)
0116    lp = length(points);
0117    ls = length(shape{i});
0118    if i <= numel(elec_pos) && ~isempty(elec_pos{i})
0119        [pp e_idx elec_pos{i} e_ref] = integrate_elecs(shape{i},elec_pos{i},elec_shape{i});
0120        ls = length(pp);
0121        points  = [points; pp];
0122    else
0123        e_idx = zeros(1,length(shape{i}));
0124        e_ref = [];
0125        points = [points; shape{i}];
0126    end
0127    if ~isempty(eidx)
0128        eidx = [eidx max(double(eidx))*(e_idx>0)+e_idx];
0129    else 
0130        eidx = e_idx;
0131    end
0132    eref = [eref; e_ref];
0133    seg{i} = repmat([0 1],ls,1) + lp + repmat((1:ls)',1,2);
0134    seg{i}(end,2) = lp + 1;
0135 end
0136 
0137 fnamebase = tempname;
0138 fnamein2d = [fnamebase, '.in2d'];
0139 fnamevol =  [fnamebase, '.vol'];
0140 
0141 write_in2d_file(fnamein2d, points, seg, eidx, eref);
0142 
0143 call_netgen( fnamein2d, fnamevol);
0144 if i_wrote_ng_opt; delete('ng.opt'); end
0145 
0146 mdl = ng_mk_fwd_model( fnamevol, [], 'ng', []);
0147 
0148 delete(fnamein2d); 
0149 delete(fnamevol); 
0150 
0151 mdl.nodes(:,3) = [];
0152 if ~all(cellfun(@isempty,elec_pos))
0153     mdl = find_electrodes(mdl, points(find(eidx),:), nonzeros(eidx));
0154 end
0155 mdl.boundary = find_boundary(mdl);
0156 if isfield(mdl, 'electrode')
0157     for i = 1:length(mdl.electrode)
0158         mdl.electrode(i).z_contact = 0.01;
0159     end
0160 end
0161 
0162 mdl = eidors_obj('fwd_model',mdl,'name','ng_mk_2d_model');
0163 
0164 function [shape,i_wrote_ng_opt] = process_maxsz(shape)
0165 maxsz = [];
0166 if numel(shape{end})==1
0167     maxsz = shape{end};
0168     shape(end)=[];
0169 end
0170 if ~isempty(maxsz)
0171     ng_write_opt('meshoptions.fineness',6,'options.meshsize',maxsz);
0172     i_wrote_ng_opt = true;
0173 else
0174     i_wrote_ng_opt = false;
0175 end
0176 
0177 function mdl = find_electrodes(mdl, elec_pts, e_idx)
0178 
0179 opt.boundary_face = 1;
0180 mdl = fix_model(mdl, opt); % in case there are multi-point electrodes
0181 
0182 nel = max(e_idx);
0183 npts = size(elec_pts,1);
0184 nn  = length(mdl.nodes);
0185 e = elec_pts';
0186 d = repmat(e(:)',nn,1) - repmat(mdl.nodes,1,npts);
0187 d = sqrt(d(:,1:2:end).^2 + d(:,2:2:end).^2);
0188 for j = 1:nel
0189     epts = find(e_idx==j);
0190     for k = 1:length(epts)
0191         [val mdl.electrode(j).nodes(k)] = min(d(:,epts(k)));
0192     end
0193     if numel(mdl.electrode(j).nodes) > 1
0194         mdl.electrode(j).nodes = fill_in_elec_nodes(mdl, mdl.electrode(j).nodes);
0195     end
0196 end
0197 
0198 function nds = fill_in_elec_nodes(mdl,enodes)
0199 fcs = mdl.faces(mdl.boundary_face,:);
0200 % fcs are ordered such that all(fcs(:,1) < fcs(:,2))
0201 % we assume that enodes are also sorted
0202 nds(1) = enodes(1);
0203 for i = 1:length(enodes)-1
0204     startnode  = enodes(i);
0205     targetnode = enodes(i+1);
0206     nextnode   = startnode;
0207     while nextnode ~= targetnode
0208         % find the two faces the startnode is on
0209         % consider which of the two nodes at their other ends is closer to
0210         % targetnode
0211         idx = find(fcs(:,1) == nextnode);
0212         switch numel(idx)
0213             case 2
0214                 c1 = fcs(idx(1),2);
0215                 c2 = fcs(idx(2),2);
0216             case 1
0217                 c1 = fcs(idx(1),2);
0218                 idx(2) = find(fcs(:,2) == nextnode);
0219                 c2 = fcs(idx(2),1);
0220             case 0
0221                 idx = find(fcs(:,2) == nextnode);
0222                 c1 = fcs(idx(1),1);
0223                 c2 = fcs(idx(2),1);
0224             otherwise
0225                 error('huh?');
0226         end
0227         d1 = sqrt(sum((mdl.nodes(c1,:) - mdl.nodes(targetnode,:)).^2,2));
0228         d2 = sqrt(sum((mdl.nodes(c2,:) - mdl.nodes(targetnode,:)).^2,2));
0229         if d1 < d2
0230             nextnode = c1;
0231         else
0232             nextnode = c2;
0233         end
0234         nds(end+1) = nextnode;
0235     end
0236 end
0237     
0238 
0239     
0240 
0241 
0242 
0243 function [newpoints eidx elec_pos e_ref] = integrate_elecs(points, elec_pos, elec_shape)
0244 
0245 
0246 n_elecs = size(elec_pos,1);
0247 if n_elecs == 1
0248     % the number of electrodes was specified, positions need to be found
0249     n_elecs = elec_pos(1);
0250     start = 0;
0251     try start = elec_pos(2); end
0252     elec_pos = equidistant_elec_pos(points, n_elecs, start);
0253     n_elecs = size(elec_pos,1);
0254 end
0255 
0256 if size(elec_shape,1) == 1;
0257     elec_shape = repmat(elec_shape,n_elecs,1);
0258 end
0259 
0260 newpoints = points;
0261 eidx = zeros(1, length(points));
0262 eref = zeros(1, length(points));
0263 
0264 for i = 1:n_elecs
0265     A = newpoints;
0266     B = circshift(newpoints,-1);
0267     AB = B-A;    L = sqrt(sum((AB).^2,2));
0268 
0269     % 1. find the closest edge
0270     % 2. add between the endpoints
0271     E = elec_pos(i,:);
0272     AE = repmat(E,size(A,1),1) - A;
0273     r = sum(AE .* AB,2)./L.^2;
0274     P = A + r*[1 1].*AB; % E projected on each edge
0275     D = sqrt(sum((repmat(E, size(A,1),1)-P).^2,2));
0276     D(r>1) = Inf; D(r<0) = Inf;
0277     [jnk e] = min(D); % index of closest edge
0278     
0279     if elec_shape(i,1) == 0 % point electrode
0280         if r(e) == 0
0281             eidx(e) = i;
0282             eref(e) = elec_shape(i,2);
0283         elseif r(e) == 1
0284             if e==length(A);
0285                 eidx(1) = i;
0286                 eref(1) = elec_shape(i,2);
0287             else
0288                 eidx(e+1) = i;
0289                 eref(e+1) = elec_shape(i,2);
0290             end
0291         else
0292             newpoints = [newpoints(1:e,:); P(e,:); newpoints(e+1:end,:)];
0293             eidx      = [eidx(1:e) i eidx(e+1:end)];
0294             eref      = [eref(1:e) elec_shape(i,2) eref(e+1:end)];
0295         end
0296     else % multi-point electrode
0297         % e is the first node of the edge the centre lies on
0298         
0299         % 1. Need the perimeter coordinate of the centre
0300         p = sqrt(sum((circshift(newpoints,-1) - newpoints).^2,2));
0301         vec = [0; cumsum(p)];
0302         L = vec(end); % total length
0303         ctr = vec(e) + r(e)*(vec(e+1) - vec(e));
0304         e_fr = linspace(ctr-elec_shape(i,1)/2 , ctr+elec_shape(i,1)/2,2);
0305         e_fr = rem(e_fr, L); % wrap around
0306         e_fr(e_fr<0) = L + e_fr(e_fr<0); % wrap around
0307         for j = 1:numel(e_fr)
0308             k = find(vec <= e_fr(j), 1, 'last');
0309             if k == length(vec)
0310                 % handle the case where electrode falls on the last point
0311                 eidx(1) = i;
0312                 eref(1) = elec_shape(i,2);
0313             else
0314                 r = (e_fr(j) - vec(k)) / (vec(k+1) - vec(k));
0315                 if r == 0
0316                     eidx(k) = i;
0317                     eref(k) = elec_shape(i,2);
0318                 end
0319                 jnkpts = newpoints; jnkpts(end+1,:) = jnkpts(1,:);
0320                 p = newpoints(k,:) + r * (jnkpts(k+1,:) - newpoints(k,:));
0321                 if k < length(eidx)
0322                     newpoints = [newpoints(1:k,:); p; newpoints(k+1:end,:)];
0323                     eidx = [eidx(1:k) i eidx(k+1:end)];
0324                     eref = [eref(1:k) elec_shape(i,2) eref(k+1:end)];
0325                 else
0326                     eidx = [eidx i ];
0327                     eref = [eref elec_shape(i,2)];
0328                     newpoints = [newpoints; p];
0329                 end
0330                 vec = [vec(1:k); e_fr(j); vec(k+1:end)];
0331             end
0332         end
0333         
0334     end
0335         
0336 end
0337 e_ref = nonzeros(eref);
0338 
0339 
0340 function elec_pos = equidistant_elec_pos(points, n_elecs, start)
0341 % 1. Calculate the perimeter
0342 p = sqrt(sum((circshift(points,-1) - points).^2,2));
0343 vec = [0; cumsum(p)];
0344 L = vec(end); % total length
0345 
0346 if n_elecs > 0
0347     e_fr = linspace(start, L+start, n_elecs+1); e_fr(end) = [];
0348 else
0349     e_fr = linspace(L+start, start, -n_elecs+1); e_fr(end) = [];
0350     n_elecs = -n_elecs;
0351 end
0352 e_fr = rem(e_fr, L); % wrap around
0353 e_fr(e_fr<0) = L + e_fr(e_fr<0); % wrap around
0354 elec_pos = NaN(n_elecs,2);
0355 points(end+1,:) = points(1,:);
0356 for i = 1:n_elecs
0357     j = find(vec <= e_fr(i), 1, 'last');
0358     if j == length(vec)
0359         % handle the case where electrode falls on the last point
0360         elec_pos(i,:) = points(1);
0361     else
0362         r = (e_fr(i) - vec(j)) / (vec(j+1) - vec(j));
0363         elec_pos(i,:) = points(j,:) + r * (points(j+1,:) - points(j,:));
0364     end
0365 end
0366 
0367 
0368 
0369 function write_in2d_file(fname,points, seg, e_idx, e_ref)
0370 
0371 if length(e_idx) < length(points);
0372     e_idx(length(points)) = 0;
0373 end
0374 
0375 refine = ones(length(points),1);
0376 if any(e_idx)
0377     refine(logical(e_idx)) = e_ref;  % refinement factor (somehow)
0378 end
0379 fid = fopen(fname,'w');
0380 fprintf(fid, '%s\n','splinecurves2dv2');
0381 fprintf(fid, '%d\n',6); % global grading factor, 6 should force use of ng.opt
0382 fprintf(fid, '%s\n','points');
0383 for i = 1:length(points)
0384    fprintf(fid, '%d   %f   %f   %f\n',i,points(i,:), refine(i));
0385 end
0386 fprintf(fid,'%s\n','segments');
0387 % here we assume the first loop is the boundary, all the others are holes
0388 domains = [ 1 0];
0389 for i = 1:length(seg)
0390    if i > 1
0391       domains = [0 1];
0392    end
0393    for j = 1:length(seg{i})
0394       fprintf(fid,'%d   %d   %d   %d   %d -bc=%d\n',domains, 2, seg{i}(j,:),i);
0395    end
0396 end
0397 fclose(fid);
0398 
0399 function mdl = do_unit_test
0400 xy = [0 0;  1 0; 1 1; 0 1];
0401 for i = 1:20
0402     switch i
0403         case 1
0404             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy});
0405         case 2
0406             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [0.5 1; 0.5 0; 0 0.5]);
0407         case 3
0408             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.3]);
0409         case 4
0410             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25]);
0411         case 5
0412             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [-5, 0.25]);
0413         case 6 
0414             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, -0.25]);
0415         case 6
0416             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0417         case 7
0418             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.1}, {[5, -0.25], [4 0.1]});
0419         case 8
0420             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0421         case 9
0422             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1], [4]});
0423         case 10
0424             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]});
0425         case 11 
0426             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},[0 30]);
0427         case 12
0428             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25],[0.2,10;0 20; 0 30; 0 20; 0 10]);
0429         case 13
0430             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0],[0.2,10;0 20; 0 20; 0 20; 0 20]);
0431         case 14
0432             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},...
0433                 [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
0434         case 15
0435             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.05}, {[5, -0.25], [], [4]},...
0436                 [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
0437         case 16
0438             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0439                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0440             xy = flipud(xy);
0441             mdl = ng_mk_2d_model(xy,9,[0.05 10]);
0442         case 17
0443             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0444                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0445             xy = flipud(xy);
0446             mdl = ng_mk_2d_model(xy,9,[0.05 200]);        
0447         case 18
0448             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0449                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0450                 xy = flipud(xy);
0451             mdl = ng_mk_2d_model({xy 0.1},9,[0.05 10]);
0452         case 19
0453             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0454                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0455                 xy = flipud(xy);
0456             mdl = ng_mk_2d_model({xy 0.1},9,0.05);
0457     end
0458     show_fem(mdl,[0 1 0]);
0459     drawnow
0460 end