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 6244 2022-03-23 13:33:53Z aadler $
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 function [shape,i_wrote_ng_opt] = process_maxsz(shape)
0163 maxsz = [];
0164 if numel(shape{end})==1
0165     maxsz = shape{end};
0166     shape(end)=[];
0167 end
0168 if ~isempty(maxsz)
0169     ng_write_opt('meshoptions.fineness',6,'options.meshsize',maxsz);
0170     i_wrote_ng_opt = true;
0171 else
0172     i_wrote_ng_opt = false;
0173 end
0174 
0175 function mdl = find_electrodes(mdl, elec_pts, e_idx)
0176 
0177 opt.boundary_face = 1;
0178 mdl = fix_model(mdl, opt); % in case there are multi-point electrodes
0179 
0180 nel = max(e_idx);
0181 npts = size(elec_pts,1);
0182 nn  = length(mdl.nodes);
0183 e = elec_pts';
0184 d = repmat(e(:)',nn,1) - repmat(mdl.nodes,1,npts);
0185 d = sqrt(d(:,1:2:end).^2 + d(:,2:2:end).^2);
0186 for j = 1:nel
0187     epts = find(e_idx==j);
0188     for k = 1:length(epts)
0189         [val mdl.electrode(j).nodes(k)] = min(d(:,epts(k)));
0190     end
0191     if numel(mdl.electrode(j).nodes) > 1
0192         mdl.electrode(j).nodes = fill_in_elec_nodes(mdl, mdl.electrode(j).nodes);
0193     end
0194 end
0195 
0196 function nds = fill_in_elec_nodes(mdl,enodes)
0197 fcs = mdl.faces(mdl.boundary_face,:);
0198 % fcs are ordered such that all(fcs(:,1) < fcs(:,2))
0199 % we assume that enodes are also sorted
0200 nds(1) = enodes(1);
0201 for i = 1:length(enodes)-1
0202     startnode  = enodes(i);
0203     targetnode = enodes(i+1);
0204     nextnode   = startnode;
0205     while nextnode ~= targetnode
0206         % find the two faces the startnode is on
0207         % consider which of the two nodes at their other ends is closer to
0208         % targetnode
0209         idx = find(fcs(:,1) == nextnode);
0210         switch numel(idx)
0211             case 2
0212                 c1 = fcs(idx(1),2);
0213                 c2 = fcs(idx(2),2);
0214             case 1
0215                 c1 = fcs(idx(1),2);
0216                 idx(2) = find(fcs(:,2) == nextnode);
0217                 c2 = fcs(idx(2),1);
0218             case 0
0219                 idx = find(fcs(:,2) == nextnode);
0220                 c1 = fcs(idx(1),1);
0221                 c2 = fcs(idx(2),1);
0222             otherwise
0223                 error('huh?');
0224         end
0225         d1 = sqrt(sum((mdl.nodes(c1,:) - mdl.nodes(targetnode,:)).^2,2));
0226         d2 = sqrt(sum((mdl.nodes(c2,:) - mdl.nodes(targetnode,:)).^2,2));
0227         if d1 < d2
0228             nextnode = c1;
0229         else
0230             nextnode = c2;
0231         end
0232         nds(end+1) = nextnode;
0233     end
0234 end
0235     
0236 
0237     
0238 
0239 
0240 
0241 function [newpoints eidx elec_pos e_ref] = integrate_elecs(points, elec_pos, elec_shape)
0242 
0243 
0244 n_elecs = size(elec_pos,1);
0245 if n_elecs == 1
0246     % the number of electrodes was specified, positions need to be found
0247     n_elecs = elec_pos(1);
0248     start = 0;
0249     try start = elec_pos(2); end
0250     elec_pos = equidistant_elec_pos(points, n_elecs, start);
0251     n_elecs = size(elec_pos,1);
0252 end
0253 
0254 if size(elec_shape,1) == 1;
0255     elec_shape = repmat(elec_shape,n_elecs,1);
0256 end
0257 
0258 newpoints = points;
0259 eidx = zeros(1, length(points));
0260 eref = zeros(1, length(points));
0261 
0262 for i = 1:n_elecs
0263     A = newpoints;
0264     B = circshift(newpoints,-1);
0265     AB = B-A;    L = sqrt(sum((AB).^2,2));
0266 
0267     % 1. find the closest edge
0268     % 2. add between the endpoints
0269     E = elec_pos(i,:);
0270     AE = repmat(E,size(A,1),1) - A;
0271     r = sum(AE .* AB,2)./L.^2;
0272     P = A + r*[1 1].*AB; % E projected on each edge
0273     D = sqrt(sum((repmat(E, size(A,1),1)-P).^2,2));
0274     D(r>1) = Inf; D(r<0) = Inf;
0275     [jnk e] = min(D); % index of closest edge
0276     
0277     if elec_shape(i,1) == 0 % point electrode
0278         if r(e) == 0
0279             eidx(e) = i;
0280             eref(e) = elec_shape(i,2);
0281         elseif r(e) == 1
0282             if e==length(A);
0283                 eidx(1) = i;
0284                 eref(1) = elec_shape(i,2);
0285             else
0286                 eidx(e+1) = i;
0287                 eref(e+1) = elec_shape(i,2);
0288             end
0289         else
0290             newpoints = [newpoints(1:e,:); P(e,:); newpoints(e+1:end,:)];
0291             eidx      = [eidx(1:e) i eidx(e+1:end)];
0292             eref      = [eref(1:e) elec_shape(i,2) eref(e+1:end)];
0293         end
0294     else % multi-point electrode
0295         % e is the first node of the edge the centre lies on
0296         
0297         % 1. Need the perimeter coordinate of the centre
0298         p = sqrt(sum((circshift(newpoints,-1) - newpoints).^2,2));
0299         vec = [0; cumsum(p)];
0300         L = vec(end); % total length
0301         ctr = vec(e) + r(e)*(vec(e+1) - vec(e));
0302         e_fr = linspace(ctr-elec_shape(i,1)/2 , ctr+elec_shape(i,1)/2,2);
0303         e_fr = rem(e_fr, L); % wrap around
0304         e_fr(e_fr<0) = L + e_fr(e_fr<0); % wrap around
0305         for j = 1:numel(e_fr)
0306             k = find(vec <= e_fr(j), 1, 'last');
0307             if k == length(vec)
0308                 % handle the case where electrode falls on the last point
0309                 eidx(1) = i;
0310                 eref(1) = elec_shape(i,2);
0311             else
0312                 r = (e_fr(j) - vec(k)) / (vec(k+1) - vec(k));
0313                 if r == 0
0314                     eidx(k) = i;
0315                     eref(k) = elec_shape(i,2);
0316                 end
0317                 jnkpts = newpoints; jnkpts(end+1,:) = jnkpts(1,:);
0318                 p = newpoints(k,:) + r * (jnkpts(k+1,:) - newpoints(k,:));
0319                 if k < length(eidx)
0320                     newpoints = [newpoints(1:k,:); p; newpoints(k+1:end,:)];
0321                     eidx = [eidx(1:k) i eidx(k+1:end)];
0322                     eref = [eref(1:k) elec_shape(i,2) eref(k+1:end)];
0323                 else
0324                     eidx = [eidx i ];
0325                     eref = [eref elec_shape(i,2)];
0326                     newpoints = [newpoints; p];
0327                 end
0328                 vec = [vec(1:k); e_fr(j); vec(k+1:end)];
0329             end
0330         end
0331         
0332     end
0333         
0334 end
0335 e_ref = nonzeros(eref);
0336 
0337 
0338 function elec_pos = equidistant_elec_pos(points, n_elecs, start)
0339 % 1. Calculate the perimeter
0340 p = sqrt(sum((circshift(points,-1) - points).^2,2));
0341 vec = [0; cumsum(p)];
0342 L = vec(end); % total length
0343 
0344 if n_elecs > 0
0345     e_fr = linspace(start, L+start, n_elecs+1); e_fr(end) = [];
0346 else
0347     e_fr = linspace(L+start, start, -n_elecs+1); e_fr(end) = [];
0348     n_elecs = -n_elecs;
0349 end
0350 e_fr = rem(e_fr, L); % wrap around
0351 e_fr(e_fr<0) = L + e_fr(e_fr<0); % wrap around
0352 elec_pos = NaN(n_elecs,2);
0353 points(end+1,:) = points(1,:);
0354 for i = 1:n_elecs
0355     j = find(vec <= e_fr(i), 1, 'last');
0356     if j == length(vec)
0357         % handle the case where electrode falls on the last point
0358         elec_pos(i,:) = points(1);
0359     else
0360         r = (e_fr(i) - vec(j)) / (vec(j+1) - vec(j));
0361         elec_pos(i,:) = points(j,:) + r * (points(j+1,:) - points(j,:));
0362     end
0363 end
0364 
0365 
0366 
0367 function write_in2d_file(fname,points, seg, e_idx, e_ref)
0368 
0369 if length(e_idx) < length(points);
0370     e_idx(length(points)) = 0;
0371 end
0372 
0373 refine = ones(length(points),1);
0374 if any(e_idx)
0375     refine(logical(e_idx)) = e_ref;  % refinement factor (somehow)
0376 end
0377 fid = fopen(fname,'w');
0378 fprintf(fid, '%s\n','splinecurves2dv2');
0379 fprintf(fid, '%d\n',6); % global grading factor, 6 should force use of ng.opt
0380 fprintf(fid, '%s\n','points');
0381 for i = 1:length(points)
0382    fprintf(fid, '%d   %f   %f   %f\n',i,points(i,:), refine(i));
0383 end
0384 fprintf(fid,'%s\n','segments');
0385 % here we assume the first loop is the boundary, all the others are holes
0386 domains = [ 1 0];
0387 for i = 1:length(seg)
0388    if i > 1
0389       domains = [0 1];
0390    end
0391    for j = 1:length(seg{i})
0392       fprintf(fid,'%d   %d   %d   %d   %d -bc=%d\n',domains, 2, seg{i}(j,:),i);
0393    end
0394 end
0395 fclose(fid);
0396 
0397 function mdl = do_unit_test
0398 xy = [0 0;  1 0; 1 1; 0 1];
0399 for i = 16:20
0400     switch i
0401         case 1
0402             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy});
0403         case 2
0404             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [0.5 1; 0.5 0; 0 0.5]);
0405         case 3
0406             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.3]);
0407         case 4
0408             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [5, 0.25]);
0409         case 5
0410             mdl = ng_mk_2d_model({xy, 0.25 + 0.5*xy, 0.1}, [-5, 0.25]);
0411         case 6 
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], [4 0.1]});
0415         case 7
0416             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.1}, {[5, -0.25], [4 0.1]});
0417         case 8
0418             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [4 0.1]});
0419         case 9
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], [4]});
0421         case 10
0422             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]});
0423         case 11 
0424             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[5, -0.25], [], [4]},[0 30]);
0425         case 12
0426             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]);
0427         case 13
0428             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]);
0429         case 14
0430             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.1}, {[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 15
0433             mdl = ng_mk_2d_model({xy, 0.1 + 0.25*xy, 0.4 + 0.5*xy, 0.05}, {[5, -0.25], [], [4]},...
0434                 [0.2,10;0 20; 0 20; 0 20; 0 20; 0 20; 0 20; 0.2 20; 0 20]);
0435         case 16
0436             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0437                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0438             xy = flipud(xy);
0439             mdl = ng_mk_2d_model(xy,9,[0.05 10]);
0440         case 17
0441             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0442                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0443             xy = flipud(xy);
0444             mdl = ng_mk_2d_model(xy,9,[0.05 200]);        
0445         case 18
0446             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0447                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0448                 xy = flipud(xy);
0449             mdl = ng_mk_2d_model({xy 0.1},9,[0.05 10]);
0450         case 19
0451             xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0452                    0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0453                 xy = flipud(xy);
0454             mdl = ng_mk_2d_model({xy 0.1},9,0.05);
0455     end
0456     show_fem(mdl,[0 1 0]);
0457     drawnow
0458 end