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.html 2819 2011-09-07 16:43:11Z aadler $
0044 
0045 if rem(rings,4) ~= 0
0046    error('parameter rings and must be divisible by 4');
0047 end
0048 
0049 % parse easy case of electrode specifications
0050 n_elec= [];
0051 if size(elec_spec) == [1,1] if isnumeric(elec_spec)
0052    n_elec= elec_spec;
0053 end; end
0054 
0055 [elem, node, bdy, point_elec_nodes, node_order] = mk_2D_model( rings );
0056 
0057 if isempty( levels ) % 2D
0058    
0059    if ~isempty( n_elec )
0060       idx= (0:n_elec-1)*length(point_elec_nodes)/n_elec + 1;
0061       if any(rem(idx,1) ~= 0);
0062          error('The requested number of electrodes is not compatible with this FEM mesh')
0063       end
0064       elec_nodes= point_elec_nodes( idx );
0065    else
0066       error('2D models only support scalar electrode patterns');
0067    end
0068 else  %3D
0069    [elem, node, bdy, point_elec_nodes] = mk_3D_model( elem, node, ...
0070                   levels, bdy, point_elec_nodes, node_order );
0071     % 3D - fixed - don't need this anymore!
0072 %  bdy= find_boundary(elem')';
0073 
0074    if ~isempty( n_elec )
0075       idx= (0:n_elec-1)*length(point_elec_nodes)/n_elec + 1;
0076       half_lev = ceil( length(levels)/2 );
0077       elec_nodes= point_elec_nodes( half_lev, idx );
0078    else
0079       elec_nodes= electrode_pattern( point_elec_nodes, elec_spec );
0080    end
0081 
0082 end
0083 
0084 param.name= sprintf('EIT FEM by mk_circ_tank with N=%d levs=%d', ...
0085                     rings, length(levels) );
0086 param.nodes = node';
0087 param.elems = elem';
0088 param.boundary = bdy';
0089 param.gnd_node = 1; % node at bottom and center of the tank
0090 param.electrode =  mk_electrodes( elec_nodes );
0091 
0092 return;
0093 
0094 % parse the elec_spec parameter
0095 % elec_spec = { 'planes', n_elecs, elec_planes }
0096 %      - puts plane each of n_elecs at planes specified by elec_planes
0097 % eg. elec_spec  =  {'planes', 16, [2,6,8]}
0098 %
0099 % elec_spec = { 'zigzag', n_elecs, elec_planes }
0100 %      - puts plane of n_elecs 'zigzagged' electrodes onto planes specified
0101 %        1st elec on plane 2, 2nd elec on plane 6, 3rd on plane 2, etc
0102 % eg. elec_spec  =  {'zigzag', 16, [2,6]}
0103 function elec_nodes= electrode_pattern( point_elec_nodes, elec_spec )
0104    elec_nodes= [];
0105    lpe = size(point_elec_nodes,2);
0106    nlev= size(point_elec_nodes,1);
0107    for i=1:3:length(elec_spec)-2
0108       spec = elec_spec{i};
0109       if      strcmp( spec,'planes' )
0110           n_elec= elec_spec{i+1};
0111           levs =  elec_spec{i+2};
0112 
0113           eidx= (0:n_elec-1);
0114           idx= round(eidx*lpe/n_elec) + 1;
0115           nodes= point_elec_nodes( levs, idx )';
0116           elec_nodes= [ elec_nodes; nodes(:) ];
0117       elseif  strcmp( spec,'zigzag' )
0118           n_elec= elec_spec{i+1};
0119           levs =  elec_spec{i+2};
0120           if any(levs > size(point_elec_nodes,1))
0121              error('requested electrode plane larger than FEModel');
0122           end
0123 
0124           eidx= (0:n_elec-1);
0125           idx= round(eidx*lpe/n_elec)*nlev + ...
0126                levs( rem( eidx, length(levs))+1);
0127           nodes= point_elec_nodes( idx );
0128           elec_nodes= [ elec_nodes; nodes(:) ];
0129       else
0130         error('elec_spec parameter not understood');
0131       end
0132    end
0133 
0134 % Create a simple 2D regular mesh, based on N circular rings
0135 %   and n_elec electrodes
0136 function [ELEM, NODE, bdy_nodes, point_elec_nodes, NODE_order] =  ...
0137           mk_2D_model( N );
0138   ELEM=[];
0139   NODE= [0;0];
0140   NODE_order= [1];
0141   int=1;
0142   for k=1:N
0143     phi= (0:4*k-1)*pi/2/k;
0144     NODE= [NODE k/N*[sin(phi);cos(phi)]];
0145 
0146 % NODE_order for extruded 3D model      3 1 2 3 1
0147 %                                     1 2 3 1 2 3
0148     NOq= rem(k+(0:k),3)+1;
0149     NODE_order= [NODE_order, NOq([1:k, k+1:-1:2, 1:k, k+1:-1:2])];
0150 
0151     ext= 2*(k*k-k+1);
0152     idxe=[0:k-1; 1:k];
0153     idxi=[0:k-1]; 
0154     elem= [ ext+idxe, ext+2*k+[-idxe idxe], ...
0155                      ext+rem(4*k-idxe,4*k), ...
0156             ext+idxe, ext+2*k+[-idxe idxe], ...
0157                      ext+rem(4*k-idxe,4*k);
0158             int+idxi, int+2*(k-1)+[-idxi, idxi], ... 
0159             int+rem(4*(k-1)-idxi, 4*(k-1)+(k==1) ) ...
0160             ext+4*k+1+idxi, ...
0161             ext+6*k+ [1-idxi 3+idxi], ...
0162             ext+8*k+3-idxi ];
0163     for j=1:k
0164       r1= rem(j+k-1,3)+1;
0165       r2= rem(j+k,3)+1;
0166       r3= 6-r1-r2;
0167       elem([r1 r2 r3],j+k*(0:7) )= elem(:,j+k*(0:7));
0168     end
0169 
0170     ELEM=[ ELEM elem(:,1:(8-4*(k==N))*k) ];
0171     int=ext;
0172   end %for k=1:N
0173 
0174   bdy_nodes= [ (ext  :ext+N*4-1) ; ...
0175                (ext+1:ext+N*4-1), ext ];
0176   point_elec_nodes= (ext):(ext+N*4-1) ;
0177  
0178 
0179 % 'extrude' a 2D model defined by ELEM and NODE into a 3D model
0180 % levels are defined by 'niveaux',
0181 % 2D parameters are ELEM, NODE, and bdy
0182 %
0183 % FIXME: The boundary calculated in 3D is no good. Instead
0184 %   it needs to be fixed using find_boundary, later
0185 function [ELEM, NODE, BDY, elec_nodes] = mk_3D_model( ...
0186      elem0, node0, niveaux, bdy, elec_nodes0, node_order );
0187 
0188   elem0= node_reorder( elem0, node_order);
0189 
0190   d= size(elem0,1);       %dimentions+1
0191   n= size(node0,2);       %NODEs
0192   e= size(elem0,2);       %ELEMents
0193 
0194 %                   D     U
0195   elem_odd= [elem0([3,2,1,1],:), ... % 1 up 1 2 3 down
0196              elem0([3,2,2,1],:), ... % 1 2 up 2 3 down 
0197              elem0([3,3,2,1],:)];    % 1 2 3 up 3 down
0198   elem_even=[elem0([1,2,3,3],:), ... % 3 up 1 2 3 down
0199              elem0([1,2,2,3],:), ... % 3 2 up 2 1 down 
0200              elem0([1,1,2,3],:)];    % 3 2 1 up 1 down
0201 
0202   NODE= [node0; niveaux(1)*ones(1,n) ];
0203   ELEM= [];
0204   bl= size(bdy,2);
0205 % Interlaced bdy idx
0206 
0207   bdy_order =node_order(bdy);
0208   bdy_up= find(bdy_order>[1;1]*min(bdy_order));
0209   bdy_dn= find(bdy_order<[1;1]*max(bdy_order));
0210   
0211   bdy_odd = [bdy; bdy(bdy_up')];
0212   bdy_even= [bdy; bdy(bdy_dn')];
0213   BDY = [];
0214  
0215   ln= length(niveaux);
0216   for k=2:ln
0217     NODE=[NODE  [node0; niveaux(k)*ones(1,n)] ];
0218     if rem(k,2)==1
0219         elem= elem_odd;
0220         bdy_e0= bdy_even;
0221         bdy_e1= bdy_odd;
0222     else
0223         elem= elem_even;
0224         bdy_e1= bdy_even;
0225         bdy_e0= bdy_odd;
0226     end
0227     el_add = (k-2)*n+[[zeros(3,e);n*ones(1,e)], ...
0228                       [zeros(2,e);n*ones(2,e)], ...
0229                       [zeros(1,e);n*ones(3,e)]];
0230     ELEM= [ELEM,elem + el_add];
0231     BDY= [BDY, bdy_e0+(k-2)*n+[zeros(2,bl);n*ones(1,bl)], ...
0232                bdy_e1+(k-2)*n+[n*ones(2,bl);zeros(1,bl)] ];
0233   end %for k
0234 
0235   % Now add top and bottom boundary
0236   BDY= [elem0, BDY, elem0+n*(ln-1) ];
0237 
0238   % elec_nodes is all nodes for all layers
0239   elec_nodes= ones(ln,1) * elec_nodes0 + ...
0240               (0:ln-1)'  * n*ones(1, length(elec_nodes0) );
0241 
0242 
0243 %param.electrode = mk_electrodes( elec_nodes );
0244 % Create the electrode structure from elec_nodes
0245 % Currently implements point electrodes with
0246 %   contact impedance of near zero.
0247 function elec_struct = mk_electrodes( elec_nodes)
0248    for i= 1:length( elec_nodes )
0249       elec_struct(i).nodes     = elec_nodes(i);
0250       elec_struct(i).z_contact = 0.001; % corresponds to 1 ohm
0251    end
0252    % Need to do it this way to be compatible accross versions
0253    if ~exist('elec_struct');
0254        elec_struct= [];
0255    end
0256 
0257 function elem=  node_reorder( elem0, node_order);
0258   e= size(elem0,2);       %ELEMents
0259 
0260   no_test=  node_order(elem0);
0261   no_test=  (0:e-1)'*[3,3,3]+no_test';
0262   elem=     elem0(no_test');
0263 
0264   no_test = node_order(elem);
0265   ok= ~norm(no_test - [1;2;3]*ones(1,e));
0266 
0267   if ~ok; error('test_node_order fails - cant do 3D meshes'); end