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

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 
0064 
0065 % (C) 2012-2013 Bartlomiej Grychtol, (C) 2013 Alistair Boyle, License: GPL version 2 or version 3
0066 % $Id: ng_mk_2d_model.m 5560 2017-06-18 18:19:20Z bgrychtol $
0067 
0068 
0069 if ischar(varargin{1}) && strcmp(varargin{1}, 'UNIT_TEST'), mdl = do_unit_test; return, end 
0070 
0071 [shape, elec_pos, elec_shape] = process_input(varargin{:});
0072 
0073 mdl = eidors_cache(@ng_mk_2d_model_do,{shape, elec_pos, elec_shape},'ng_mk_2d_model');
0074 
0075 
0076 
0077 function [shape, elec_pos, elec_shape] = process_input(shape, elec_pos, elec_shape)
0078 
0079 if ~iscell(shape)
0080    shape = {shape};
0081 end
0082 
0083 if nargin < 2
0084     elec_pos = [];
0085 end
0086 if ~iscell(elec_pos)
0087     elec_pos = {elec_pos};
0088 end
0089 
0090 if nargin < 3
0091     elec_shape = [0 10]; % point electrode
0092 end
0093 if size(elec_shape,2) == 1
0094    warning('Refinement factor not specified, using 10');
0095    elec_shape(:,2) = 10;
0096 end
0097 if ~iscell(elec_shape)
0098     elec_shape = {elec_shape};
0099 end
0100 if numel(elec_shape) == 1 && numel(elec_pos) > 1
0101     elec_shape(2:numel(elec_pos)) = elec_shape(1);
0102 end
0103 
0104 
0105 function mdl = ng_mk_2d_model_do(shape, elec_pos, elec_shape)
0106 
0107 [shape,i_wrote_ng_opt] = process_maxsz(shape);
0108 
0109 points = [];
0110 eidx = [];
0111 eref = [];
0112 for i = 1:length(shape)
0113    lp = length(points);
0114    ls = length(shape{i});
0115    if i <= numel(elec_pos) && ~isempty(elec_pos{i})
0116        [pp e_idx elec_pos{i} e_ref] = integrate_elecs(shape{i},elec_pos{i},elec_shape{i});
0117        ls = length(pp);
0118        points  = [points; pp];
0119    else
0120        e_idx = zeros(1,length(shape{i}));
0121        e_ref = [];
0122        points = [points; shape{i}];
0123    end
0124    if ~isempty(eidx)
0125        eidx = [eidx max(double(eidx))*(e_idx>0)+e_idx];
0126    else 
0127        eidx = e_idx;
0128    end
0129    eref = [eref; e_ref];
0130    seg{i} = repmat([0 1],ls,1) + lp + repmat((1:ls)',1,2);
0131    seg{i}(end,2) = lp + 1;
0132 end
0133 
0134 fnamebase = tempname;
0135 fnamein2d = [fnamebase, '.in2d'];
0136 fnamevol =  [fnamebase, '.vol'];
0137 
0138 write_in2d_file(fnamein2d, points, seg, eidx, eref);
0139 
0140 call_netgen( fnamein2d, fnamevol);
0141 if i_wrote_ng_opt; delete('ng.opt'); end
0142 
0143 mdl = ng_mk_fwd_model( fnamevol, [], 'ng', []);
0144 
0145 delete(fnamein2d); 
0146 delete(fnamevol); 
0147 
0148 mdl.nodes(:,3) = [];
0149 if ~all(cellfun(@isempty,elec_pos))
0150     mdl = find_electrodes(mdl, points(find(eidx),:), nonzeros(eidx));
0151 end
0152 mdl.boundary = find_boundary(mdl);
0153 if isfield(mdl, 'electrode')
0154     for i = 1:length(mdl.electrode)
0155         mdl.electrode(i).z_contact = 0.01;
0156     end
0157 end
0158 
0159 function [shape,i_wrote_ng_opt] = process_maxsz(shape)
0160 maxsz = [];
0161 if numel(shape{end})==1
0162     maxsz = shape{end};
0163     shape(end)=[];
0164 end
0165 if ~isempty(maxsz)
0166     ng_write_opt('meshoptions.fineness',6,'options.meshsize',maxsz);
0167     i_wrote_ng_opt = true;
0168 else
0169     i_wrote_ng_opt = false;
0170 end
0171 
0172 function mdl = find_electrodes(mdl, elec_pts, e_idx)
0173 
0174 opt.boundary_face = 1;
0175 mdl = fix_model(mdl, opt); % in case there are multi-point electrodes
0176 
0177 nel = max(e_idx);
0178 npts = size(elec_pts,1);
0179 nn  = length(mdl.nodes);
0180 e = elec_pts';
0181 d = repmat(e(:)',nn,1) - repmat(mdl.nodes,1,npts);
0182 d = sqrt(d(:,1:2:end).^2 + d(:,2:2:end).^2);
0183 for j = 1:nel
0184     epts = find(e_idx==j);
0185     for k = 1:length(epts)
0186         [val mdl.electrode(j).nodes(k)] = min(d(:,epts(k)));
0187     end
0188     if numel(mdl.electrode(j).nodes) > 1
0189         mdl.electrode(j).nodes = fill_in_elec_nodes(mdl, mdl.electrode(j).nodes);
0190     end
0191 end
0192 
0193 function nds = fill_in_elec_nodes(mdl,enodes)
0194 fcs = mdl.faces(mdl.boundary_face,:);
0195 % fcs are ordered such that all(fcs(:,1) < fcs(:,2))
0196 % we assume that enodes are also sorted
0197 nds(1) = enodes(1);
0198 for i = 1:length(enodes)-1
0199     startnode  = enodes(i);
0200     targetnode = enodes(i+1);
0201     nextnode   = startnode;
0202     while nextnode ~= targetnode
0203         % find the two faces the startnode is on
0204         % consider which of the two nodes at their other ends is closer to
0205         % targetnode
0206         idx = find(fcs(:,1) == nextnode);
0207         switch numel(idx)
0208             case 2
0209                 c1 = fcs(idx(1),2);
0210                 c2 = fcs(idx(2),2);
0211             case 1
0212                 c1 = fcs(idx(1),2);
0213                 idx(2) = find(fcs(:,2) == nextnode);
0214                 c2 = fcs(idx(2),1);
0215             case 0
0216                 idx = find(fcs(:,2) == nextnode);
0217                 c1 = fcs(idx(1),1);
0218                 c2 = fcs(idx(2),1);
0219             otherwise
0220                 error('huh?');
0221         end
0222         d1 = sqrt(sum((mdl.nodes(c1,:) - mdl.nodes(targetnode,:)).^2,2));
0223         d2 = sqrt(sum((mdl.nodes(c2,:) - mdl.nodes(targetnode,:)).^2,2));
0224         if d1 < d2
0225             nextnode = c1;
0226         else
0227             nextnode = c2;
0228         end
0229         nds(end+1) = nextnode;
0230     end
0231 end
0232     
0233 
0234     
0235 
0236 
0237 
0238 function [newpoints eidx elec_pos e_ref] = integrate_elecs(points, elec_pos, elec_shape)
0239 
0240 
0241 n_elecs = size(elec_pos,1);
0242 if n_elecs == 1
0243     % the number of electrodes was specified, positions need to be found
0244     n_elecs = elec_pos(1);
0245     start = 0;
0246     try start = elec_pos(2); end
0247     elec_pos = equidistant_elec_pos(points, n_elecs, start);
0248     n_elecs = size(elec_pos,1);
0249 end
0250 
0251 if size(elec_shape,1) == 1;
0252     elec_shape = repmat(elec_shape,n_elecs,1);
0253 end
0254 
0255 newpoints = points;
0256 eidx = zeros(1, length(points));
0257 eref = zeros(1, length(points));
0258 
0259 for i = 1:n_elecs
0260     A = newpoints;
0261     B = circshift(newpoints,-1);
0262     AB = B-A;    L = sqrt(sum((AB).^2,2));
0263 
0264     % 1. find the closest edge
0265     % 2. add between the endpoints
0266     E = elec_pos(i,:);
0267     AE = repmat(E,size(A,1),1) - A;
0268     r = sum(AE .* AB,2)./L.^2;
0269     P = A + r*[1 1].*AB; % E projected on each edge
0270     D = sqrt(sum((repmat(E, size(A,1),1)-P).^2,2));
0271     D(r>1) = Inf; D(r<0) = Inf;
0272     [jnk e] = min(D); % index of closest edge
0273     
0274     if elec_shape(i,1) == 0 % point electrode
0275         if r(e) == 0
0276             eidx(e) = i;
0277             eref(e) = elec_shape(i,2);
0278         elseif r(e) == 1
0279             if e==length(A);
0280                 eidx(1) = i;
0281                 eref(1) = elec_shape(i,2);
0282             else
0283                 eidx(e+1) = i;
0284                 eref(e+1) = elec_shape(i,2);
0285             end
0286         else
0287             newpoints = [newpoints(1:e,:); P(e,:); newpoints(e+1:end,:)];
0288             eidx      = [eidx(1:e) i eidx(e+1:end)];
0289             eref      = [eref(1:e) elec_shape(i,2) eref(e+1:end)];
0290         end
0291     else % multi-point electrode
0292         % e is the first node of the edge the centre lies on
0293         
0294         % 1. Need the perimeter coordinate of the centre
0295         p = sqrt(sum((circshift(newpoints,-1) - newpoints).^2,2));
0296         vec = [0; cumsum(p)];
0297         L = vec(end); % total length
0298         ctr = vec(e) + r(e)*(vec(e+1) - vec(e));
0299         e_fr = linspace(ctr-elec_shape(i,1)/2 , ctr+elec_shape(i,1)/2,2);
0300         e_fr = rem(e_fr, L); % wrap around
0301         e_fr(e_fr<0) = L + e_fr(e_fr<0); % wrap around
0302         for j = 1:numel(e_fr)
0303             k = find(vec <= e_fr(j), 1, 'last');
0304             if k == length(vec)
0305                 % handle the case where electrode falls on the last point
0306                 eidx(1) = i;
0307                 eref(1) = elec_shape(i,2);
0308             else
0309                 r = (e_fr(j) - vec(k)) / (vec(k+1) - vec(k));
0310                 if r == 0
0311                     eidx(k) = i;
0312                     eref(k) = elec_shape(i,2);
0313                 end
0314                 jnkpts = newpoints; jnkpts(end+1,:) = jnkpts(1,:);
0315                 p = newpoints(k,:) + r * (jnkpts(k+1,:) - newpoints(k,:));
0316                 if k < length(eidx)
0317                     newpoints = [newpoints(1:k,:); p; newpoints(k+1:end,:)];
0318                     eidx = [eidx(1:k) i eidx(k+1:end)];
0319                     eref = [eref(1:k) elec_shape(i,2) eref(k+1:end)];
0320                 else
0321                     eidx = [eidx i ];
0322                     eref = [eref elec_shape(i,2)];
0323                     newpoints = [newpoints; p];
0324                 end
0325                 vec = [vec(1:k); e_fr(j); vec(k+1:end)];
0326             end
0327         end
0328         
0329     end
0330         
0331 end
0332 e_ref = nonzeros(eref);
0333 
0334 
0335 function elec_pos = equidistant_elec_pos(points, n_elecs, start)
0336 % 1. Calculate the perimeter
0337 p = sqrt(sum((circshift(points,-1) - points).^2,2));
0338 vec = [0; cumsum(p)];
0339 L = vec(end); % total length
0340 
0341 if n_elecs > 0
0342     e_fr = linspace(start, L+start, n_elecs+1); e_fr(end) = [];
0343 else
0344     e_fr = linspace(L+start, start, -n_elecs+1); e_fr(end) = [];
0345     n_elecs = -n_elecs;
0346 end
0347 e_fr = rem(e_fr, L); % wrap around
0348 e_fr(e_fr<0) = L + e_fr(e_fr<0); % wrap around
0349 elec_pos = NaN(n_elecs,2);
0350 points(end+1,:) = points(1,:);
0351 for i = 1:n_elecs
0352     j = find(vec <= e_fr(i), 1, 'last');
0353     if j == length(vec)
0354         % handle the case where electrode falls on the last point
0355         elec_pos(i,:) = points(1);
0356     else
0357         r = (e_fr(i) - vec(j)) / (vec(j+1) - vec(j));
0358         elec_pos(i,:) = points(j,:) + r * (points(j+1,:) - points(j,:));
0359     end
0360 end
0361 
0362 
0363 
0364 function write_in2d_file(fname,points, seg, e_idx, e_ref)
0365 
0366 if length(e_idx) < length(points);
0367     e_idx(length(points)) = 0;
0368 end
0369 
0370 refine = ones(length(points),1);
0371 if any(e_idx)
0372     refine(logical(e_idx)) = e_ref;  % refinement factor (somehow)
0373 end
0374 fid = fopen(fname,'w');
0375 fprintf(fid, '%s\n','splinecurves2dv2');
0376 fprintf(fid, '%d\n',6); % global grading factor, 6 should force use of ng.opt
0377 fprintf(fid, '%s\n','points');
0378 for i = 1:length(points)
0379    fprintf(fid, '%d   %f   %f   %f\n',i,points(i,:), refine(i));
0380 end
0381 fprintf(fid,'%s\n','segments');
0382 % here we assume the first loop is the boundary, all the others are holes
0383 domains = [ 1 0];
0384 for i = 1:length(seg)
0385    if i > 1
0386       domains = [0 1];
0387    end
0388    for j = 1:length(seg{i})
0389       fprintf(fid,'%d   %d   %d   %d   %d -bc=%d\n',domains, 2, seg{i}(j,:),i);
0390    end
0391 end
0392 fclose(fid);
0393 
0394 function mdl = do_unit_test
0395 xy = [0 0;  1 0; 1 1; 0 1];
0396 for i = 16:20
0397     switch i
0398         case 1
0399             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy});
0400         case 2
0401             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [0.5 1; 0.5 0; 0 0.5]);
0402         case 3
0403             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.3]);
0404         case 4
0405             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25]);
0406         case 5
0407             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [-5, 0.25]);
0408         case 6 
0409             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, -0.25]);
0410         case 6
0411             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0412         case 7
0413             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.1}, {[5, -0.25], [4 0.1]});
0414         case 8
0415             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0416         case 9
0417             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]});
0418         case 10
0419             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]});
0420         case 11 
0421             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},[0 30]);
0422         case 12
0423             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]);
0424         case 13
0425             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]);
0426         case 14
0427             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},...
0428                 [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
0429         case 15
0430             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.05}, {[5, -0.25], [], [4]},...
0431                 [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
0432         case 16
0433             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0434                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0435             xy = flipud(xy);
0436             mdl = ng_mk_2d_model(xy,9,[0.05 10]);
0437         case 17
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 200]);        
0442         case 18
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 0.1},9,[0.05 10]);
0447         case 19
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);
0452     end
0453     show_fem(mdl,[0 1 0]);
0454     drawnow
0455 end