mk_circ_tank

MK_CIRC_TANK: make a cylindrical tank FEM geometry in 2D or 3D
 param= mk_circ_tank(rings, levels, elec_spec );
 
 rings:  number of horizontal plane rings (divisible by 4)
 levels: vector of vertical placement of levels
     for 2D mesh, levels = []
 
 elec_spec: parameter to specify number of electrodes
        specified as { 'opt1', val11, val12 , 'opt2', val21, val22 }

 elec_spec = scalar (divisible by 4)
      - puts a single plane of electrodes in centre of cylinder
 eg. elec_spec  = 16

 elec_spec = { 'planes', n_elecs, elec_planes }
      - puts plane each of n_elecs at planes specified by elec_planes
 eg. elec_spec  =  {'planes', 16, [2,6,8]}

 elec_spec = { 'zigzag', n_elecs, elec_planes }
      - puts plane of n_elecs 'zigzagged' electrodes onto planes specified
        1st elec on plane 2, 2nd elec on plane 6, 3rd on plane 2, etc 
 eg. elec_spec  =  {'zigzag', 16, [2,6]}
      - Note, based on the restults of Graham et al (2006), zigzag
        electrode placement is not recommended
      - In order to implement the 'planar3d' pattern from this paper,
        puts 2d electrodes onto rings ie [ ...  7  8  1  2  ...
                                           ... 15 16  9 10  ... ]
      ->use  elec_spec = { 'planes', n_elecs/2, elec_planes }

 mk_circ_tank creates simple, point electrodes. Improved models
  may be created with ng_mk_cyl_models

 output:
  param.name        Model name (if known) 
  param.nodes       position of FEM nodes (Nodes x Dims) 
  param.elems       definition of FEM elements (Elems x Dims+1) 
  param.boundary    nodes of element faces on the medium surface 
  param.gnd_node    Number of node connected to ground 
  param.electrode   Vector (Num_elecs x 1) of electrode models (elec_model)

0001 function param= mk_circ_tank(rings, levels, elec_spec );
0002 %MK_CIRC_TANK: make a cylindrical tank FEM geometry in 2D or 3D
0003 % param= mk_circ_tank(rings, levels, elec_spec );
0004 %
0005 % rings:  number of horizontal plane rings (divisible by 4)
0006 % levels: vector of vertical placement of levels
0007 %     for 2D mesh, levels = []
0008 %
0009 % elec_spec: parameter to specify number of electrodes
0010 %        specified as { 'opt1', val11, val12 , 'opt2', val21, val22 }
0011 %
0012 % elec_spec = scalar (divisible by 4)
0013 %      - puts a single plane of electrodes in centre of cylinder
0014 % eg. elec_spec  = 16
0015 %
0016 % elec_spec = { 'planes', n_elecs, elec_planes }
0017 %      - puts plane each of n_elecs at planes specified by elec_planes
0018 % eg. elec_spec  =  {'planes', 16, [2,6,8]}
0019 %
0020 % elec_spec = { 'zigzag', n_elecs, elec_planes }
0021 %      - puts plane of n_elecs 'zigzagged' electrodes onto planes specified
0022 %        1st elec on plane 2, 2nd elec on plane 6, 3rd on plane 2, etc
0023 % eg. elec_spec  =  {'zigzag', 16, [2,6]}
0024 %      - Note, based on the restults of Graham et al (2006), zigzag
0025 %        electrode placement is not recommended
0026 %      - In order to implement the 'planar3d' pattern from this paper,
0027 %        puts 2d electrodes onto rings ie [ ...  7  8  1  2  ...
0028 %                                           ... 15 16  9 10  ... ]
0029 %      ->use  elec_spec = { 'planes', n_elecs/2, elec_planes }
0030 %
0031 % mk_circ_tank creates simple, point electrodes. Improved models
0032 %  may be created with ng_mk_cyl_models
0033 %
0034 % output:
0035 %  param.name        Model name (if known)
0036 %  param.nodes       position of FEM nodes (Nodes x Dims)
0037 %  param.elems       definition of FEM elements (Elems x Dims+1)
0038 %  param.boundary    nodes of element faces on the medium surface
0039 %  param.gnd_node    Number of node connected to ground
0040 %  param.electrode   Vector (Num_elecs x 1) of electrode models (elec_model)
0041 
0042 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0043 % $Id: mk_circ_tank.m 5112 2015-06-14 13:00:41Z aadler $
0044 
0045 if ischar(rings) && strcmp(rings,'UNIT_TEST'); do_unit_test; return; end
0046 
0047 
0048 % parse easy case of electrode specifications
0049 n_elec= [];
0050 if size(elec_spec) == [1,1] if isnumeric(elec_spec)
0051    n_elec= elec_spec;
0052 end; end
0053 
0054 [elem, node, bdy, point_elec_nodes, node_order] = mk_2D_model( rings );
0055 
0056 if isempty( levels ) % 2D
0057    
0058    if n_elec==0
0059       elec_nodes= [];
0060    elseif ~isempty( n_elec )
0061       idx= (0:n_elec-1)*length(point_elec_nodes)/n_elec + 1;
0062       if any(rem(idx,1) ~= 0);
0063          error('The requested number of electrodes (%d) is not compatible with this FEM mesh', n_elec)
0064       end
0065       elec_nodes= point_elec_nodes( idx );
0066    else
0067       error('2D models only support scalar electrode patterns');
0068    end
0069 else  %3D
0070    [elem, node, bdy, point_elec_nodes] = mk_3D_model( elem, node, ...
0071                   levels, bdy, point_elec_nodes, node_order );
0072 
0073    if n_elec==0
0074       elec_nodes= [];
0075    elseif ~isempty( n_elec )
0076       idx= (0:n_elec-1)*length(point_elec_nodes)/n_elec + 1;
0077       half_lev = ceil( length(levels)/2 );
0078       elec_nodes= point_elec_nodes( half_lev, idx );
0079    else
0080       elec_nodes= electrode_pattern( point_elec_nodes, elec_spec );
0081    end
0082 
0083 end
0084 
0085 param.name= sprintf('EIT FEM by mk_circ_tank with N=%d levs=%d', ...
0086                     rings, length(levels) );
0087 param.nodes = node';
0088 param.elems = elem';
0089 param.boundary = bdy';
0090 param.gnd_node = 1; % node at bottom and center of the tank
0091 if ~isempty( elec_nodes)
0092    param.electrode =  mk_electrodes( elec_nodes );
0093 end
0094 param.type = 'fwd_model';
0095 
0096 param.normalize_measurements = 0; % default
0097 
0098 return;
0099 
0100 % parse the elec_spec parameter
0101 % elec_spec = { 'planes', n_elecs, elec_planes }
0102 %      - puts plane each of n_elecs at planes specified by elec_planes
0103 % eg. elec_spec  =  {'planes', 16, [2,6,8]}
0104 %
0105 % elec_spec = { 'zigzag', n_elecs, elec_planes }
0106 %      - puts plane of n_elecs 'zigzagged' electrodes onto planes specified
0107 %        1st elec on plane 2, 2nd elec on plane 6, 3rd on plane 2, etc
0108 % eg. elec_spec  =  {'zigzag', 16, [2,6]}
0109 function elec_nodes= electrode_pattern( point_elec_nodes, elec_spec )
0110    elec_nodes= [];
0111    lpe = size(point_elec_nodes,2);
0112    nlev= size(point_elec_nodes,1);
0113    for i=1:3:length(elec_spec)-2
0114       spec = elec_spec{i};
0115       if      strcmp( spec,'planes' )
0116           n_elec= elec_spec{i+1};
0117           levs =  elec_spec{i+2};
0118 
0119           eidx= (0:n_elec-1);
0120           idx= round(eidx*lpe/n_elec) + 1;
0121           nodes= point_elec_nodes( levs, idx )';
0122           elec_nodes= [ elec_nodes; nodes(:) ];
0123       elseif  strcmp( spec,'zigzag' )
0124           n_elec= elec_spec{i+1};
0125           levs =  elec_spec{i+2};
0126           if any(levs > size(point_elec_nodes,1))
0127              error('requested electrode plane larger than FEModel');
0128           end
0129 
0130           eidx= (0:n_elec-1);
0131           idx= round(eidx*lpe/n_elec)*nlev + ...
0132                levs( rem( eidx, length(levs))+1);
0133           nodes= point_elec_nodes( idx );
0134           elec_nodes= [ elec_nodes; nodes(:) ];
0135       else
0136         error('elec_spec parameter not understood');
0137       end
0138    end
0139 
0140 % Create a simple 2D regular mesh, based on N circular rings
0141 %   and n_elec electrodes
0142 function [ELEM, NODE, bdy_nodes, point_elec_nodes, NODE_order] =  ...
0143           mk_2D_model( N );
0144   ELEM=[];
0145   NODE= [0;0];
0146   NODE_order= [1];
0147   int=1;
0148   for k=1:N
0149     phi= (0:4*k-1)*pi/2/k;
0150     NODE= [NODE k/N*[sin(phi);cos(phi)]];
0151 
0152 % NODE_order for extruded 3D model      3 1 2 3 1
0153 %                                     1 2 3 1 2 3
0154     NOq= rem(k+(0:k),3)+1;
0155     NODE_order= [NODE_order, NOq([1:k, k+1:-1:2, 1:k, k+1:-1:2])];
0156 
0157     ext= 2*(k*k-k+1);
0158     idxe=[0:k-1; 1:k];
0159     idxi=[0:k-1]; 
0160     elem= [ ext+idxe, ext+2*k+[-idxe idxe], ...
0161                      ext+rem(4*k-idxe,4*k), ...
0162             ext+idxe, ext+2*k+[-idxe idxe], ...
0163                      ext+rem(4*k-idxe,4*k);
0164             int+idxi, int+2*(k-1)+[-idxi, idxi], ... 
0165             int+rem(4*(k-1)-idxi, 4*(k-1)+(k==1) ) ...
0166             ext+4*k+1+idxi, ...
0167             ext+6*k+ [1-idxi 3+idxi], ...
0168             ext+8*k+3-idxi ];
0169     for j=1:k
0170       r1= rem(j+k-1,3)+1;
0171       r2= rem(j+k,3)+1;
0172       r3= 6-r1-r2;
0173       elem([r1 r2 r3],j+k*(0:7) )= elem(:,j+k*(0:7));
0174     end
0175 
0176     ELEM=[ ELEM elem(:,1:(8-4*(k==N))*k) ];
0177     int=ext;
0178   end %for k=1:N
0179 
0180   bdy_nodes= [ (ext  :ext+N*4-1) ; ...
0181                (ext+1:ext+N*4-1), ext ];
0182   point_elec_nodes= (ext):(ext+N*4-1) ;
0183  
0184 
0185 % 'extrude' a 2D model defined by ELEM and NODE into a 3D model
0186 % levels are defined by 'niveaux',
0187 % 2D parameters are ELEM, NODE, and bdy
0188 %
0189 % FIXME: The boundary calculated in 3D is no good. Instead
0190 %   it needs to be fixed using find_boundary, later
0191 function [ELEM, NODE, BDY, elec_nodes] = mk_3D_model( ...
0192      elem0, node0, niveaux, bdy, elec_nodes0, node_order );
0193 
0194   elem0= node_reorder( elem0, node_order);
0195 
0196   d= size(elem0,1);       %dimentions+1
0197   n= size(node0,2);       %NODEs
0198   e= size(elem0,2);       %ELEMents
0199 
0200 %                   D     U
0201   elem_odd= [elem0([3,2,1,1],:), ... % 1 up 1 2 3 down
0202              elem0([3,2,2,1],:), ... % 1 2 up 2 3 down 
0203              elem0([3,3,2,1],:)];    % 1 2 3 up 3 down
0204   elem_even=[elem0([1,2,3,3],:), ... % 3 up 1 2 3 down
0205              elem0([1,2,2,3],:), ... % 3 2 up 2 1 down 
0206              elem0([1,1,2,3],:)];    % 3 2 1 up 1 down
0207 
0208   NODE= [node0; niveaux(1)*ones(1,n) ];
0209   ELEM= [];
0210   bl= size(bdy,2);
0211 % Interlaced bdy idx
0212 
0213   bdy_order =node_order(bdy);
0214   bdy_up= find(bdy_order>[1;1]*min(bdy_order));
0215   bdy_dn= find(bdy_order<[1;1]*max(bdy_order));
0216   
0217   bdy_odd = [bdy; bdy(bdy_up')];
0218   bdy_even= [bdy; bdy(bdy_dn')];
0219   BDY = [];
0220  
0221   ln= length(niveaux);
0222   for k=2:ln
0223     NODE=[NODE  [node0; niveaux(k)*ones(1,n)] ];
0224     if rem(k,2)==1
0225         elem= elem_odd;
0226         bdy_e0= bdy_even;
0227         bdy_e1= bdy_odd;
0228     else
0229         elem= elem_even;
0230         bdy_e1= bdy_even;
0231         bdy_e0= bdy_odd;
0232     end
0233     el_add = (k-2)*n+[[zeros(3,e);n*ones(1,e)], ...
0234                       [zeros(2,e);n*ones(2,e)], ...
0235                       [zeros(1,e);n*ones(3,e)]];
0236     ELEM= [ELEM,elem + el_add];
0237     BDY= [BDY, bdy_e0+(k-2)*n+[zeros(2,bl);n*ones(1,bl)], ...
0238                bdy_e1+(k-2)*n+[n*ones(2,bl);zeros(1,bl)] ];
0239   end %for k
0240 
0241   % Now add top and bottom boundary
0242   BDY= [elem0, BDY, elem0+n*(ln-1) ];
0243 
0244   % elec_nodes is all nodes for all layers
0245   elec_nodes= ones(ln,1) * elec_nodes0 + ...
0246               (0:ln-1)'  * n*ones(1, length(elec_nodes0) );
0247 
0248 
0249 %param.electrode = mk_electrodes( elec_nodes );
0250 % Create the electrode structure from elec_nodes
0251 % Currently implements point electrodes with
0252 %   contact impedance of near zero.
0253 function elec_struct = mk_electrodes( elec_nodes)
0254    for i= 1:length( elec_nodes )
0255       elec_struct(i).nodes     = elec_nodes(i);
0256       elec_struct(i).z_contact = 0.001; % corresponds to 1 ohm
0257    end
0258    % Need to do it this way to be compatible accross versions
0259    if ~exist('elec_struct');
0260        elec_struct= [];
0261    end
0262 
0263 function elem=  node_reorder( elem0, node_order);
0264   e= size(elem0,2);       %ELEMents
0265 
0266   no_test=  node_order(elem0);
0267   no_test=  (0:e-1)'*[3,3,3]+no_test';
0268   elem=     elem0(no_test');
0269 
0270   no_test = node_order(elem);
0271   ok= ~norm(no_test - [1;2;3]*ones(1,e));
0272 
0273   if ~ok; error('test_node_order fails - cant do 3D meshes'); end
0274 
0275 function do_unit_test
0276   subplot(3,3,1)
0277   mdl= mk_circ_tank(2, [], 2 );
0278   show_fem(mdl, [0,1,1]);
0279   unit_test_cmp('2D mdl', length(mdl.elems), 16);
0280 
0281   subplot(3,3,2)
0282   mdl= mk_circ_tank(4, [], 16 );
0283   show_fem(mdl);
0284   unit_test_cmp('2D mdl', length(mdl.nodes), 41);
0285 
0286   subplot(3,3,3)
0287   mdl= mk_circ_tank(4, linspace(-1,1,3), {'planes', 8, 2} );
0288   show_fem(mdl);
0289   unit_test_cmp('2D mdl', length(mdl.elems), 384);
0290 
0291   subplot(3,3,4)
0292   mdl= mk_circ_tank(4, linspace(-1,1,5), {'zigzag', 4, [2,4]} );
0293   show_fem(mdl);
0294   unit_test_cmp('2D mdl', length(mdl.elems), 768);
0295 
0296   try 
0297      mdl=  mk_circ_tank(2, [], 18 );  error
0298      unit_test_cmp('test for error', 1,0);
0299   catch
0300      unit_test_cmp('test for error', 1,1);
0301   end
0302 
0303   subplot(3,3,5)
0304   mdl= mk_circ_tank(4, [], 0);
0305   show_fem(mdl);
0306   title 'no electodes'
0307 
0308   subplot(3,3,6)
0309   mdl= mk_circ_tank(4, linspace(-1,1,5), 0);
0310   show_fem(mdl);
0311   title 'no electodes'