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 6926 2024-05-31 15:34:13Z bgrychtol $
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 
0095 param.normalize_measurements = mdl_normalize('default');
0096 
0097 param.jacobian = @eidors_default;
0098 param.solve = @eidors_default;
0099 param.system_mat = @eidors_default;
0100 
0101 param = eidors_obj('fwd_model',param);
0102 
0103 return;
0104 
0105 % parse the elec_spec parameter
0106 % elec_spec = { 'planes', n_elecs, elec_planes }
0107 %      - puts plane each of n_elecs at planes specified by elec_planes
0108 % eg. elec_spec  =  {'planes', 16, [2,6,8]}
0109 %
0110 % elec_spec = { 'zigzag', n_elecs, elec_planes }
0111 %      - puts plane of n_elecs 'zigzagged' electrodes onto planes specified
0112 %        1st elec on plane 2, 2nd elec on plane 6, 3rd on plane 2, etc
0113 % eg. elec_spec  =  {'zigzag', 16, [2,6]}
0114 function elec_nodes= electrode_pattern( point_elec_nodes, elec_spec )
0115    elec_nodes= [];
0116    lpe = size(point_elec_nodes,2);
0117    nlev= size(point_elec_nodes,1);
0118    for i=1:3:length(elec_spec)-2
0119       spec = elec_spec{i};
0120       if      strcmp( spec,'planes' )
0121           n_elec= elec_spec{i+1};
0122           levs =  elec_spec{i+2};
0123 
0124           eidx= (0:n_elec-1);
0125           idx= round(eidx*lpe/n_elec) + 1;
0126           nodes= point_elec_nodes( levs, idx )';
0127           elec_nodes= [ elec_nodes; nodes(:) ];
0128       elseif  strcmp( spec,'zigzag' )
0129           n_elec= elec_spec{i+1};
0130           levs =  elec_spec{i+2};
0131           if any(levs > size(point_elec_nodes,1))
0132              error('requested electrode plane larger than FEModel');
0133           end
0134 
0135           eidx= (0:n_elec-1);
0136           idx= round(eidx*lpe/n_elec)*nlev + ...
0137                levs( rem( eidx, length(levs))+1);
0138           nodes= point_elec_nodes( idx );
0139           elec_nodes= [ elec_nodes; nodes(:) ];
0140       else
0141         error('elec_spec parameter not understood');
0142       end
0143    end
0144 
0145 % Create a simple 2D regular mesh, based on N circular rings
0146 %   and n_elec electrodes
0147 function [ELEM, NODE, bdy_nodes, point_elec_nodes, NODE_order] =  ...
0148           mk_2D_model( N );
0149   ELEM=[];
0150   NODE= [0;0];
0151   NODE_order= [1];
0152   int=1;
0153   for k=1:N
0154     phi= (0:4*k-1)*pi/2/k;
0155     NODE= [NODE k/N*[sin(phi);cos(phi)]];
0156 
0157 % NODE_order for extruded 3D model      3 1 2 3 1
0158 %                                     1 2 3 1 2 3
0159     NOq= rem(k+(0:k),3)+1;
0160     NODE_order= [NODE_order, NOq([1:k, k+1:-1:2, 1:k, k+1:-1:2])];
0161 
0162     ext= 2*(k*k-k+1);
0163     idxe=[0:k-1; 1:k];
0164     idxi=[0:k-1]; 
0165     elem= [ ext+idxe, ext+2*k+[-idxe idxe], ...
0166                      ext+rem(4*k-idxe,4*k), ...
0167             ext+idxe, ext+2*k+[-idxe idxe], ...
0168                      ext+rem(4*k-idxe,4*k);
0169             int+idxi, int+2*(k-1)+[-idxi, idxi], ... 
0170             int+rem(4*(k-1)-idxi, 4*(k-1)+(k==1) ) ...
0171             ext+4*k+1+idxi, ...
0172             ext+6*k+ [1-idxi 3+idxi], ...
0173             ext+8*k+3-idxi ];
0174     for j=1:k
0175       r1= rem(j+k-1,3)+1;
0176       r2= rem(j+k,3)+1;
0177       r3= 6-r1-r2;
0178       elem([r1 r2 r3],j+k*(0:7) )= elem(:,j+k*(0:7));
0179     end
0180 
0181     ELEM=[ ELEM elem(:,1:(8-4*(k==N))*k) ];
0182     int=ext;
0183   end %for k=1:N
0184 
0185   bdy_nodes= [ (ext  :ext+N*4-1) ; ...
0186                (ext+1:ext+N*4-1), ext ];
0187   point_elec_nodes= (ext):(ext+N*4-1) ;
0188  
0189 
0190 % 'extrude' a 2D model defined by ELEM and NODE into a 3D model
0191 % levels are defined by 'niveaux',
0192 % 2D parameters are ELEM, NODE, and bdy
0193 %
0194 % FIXME: The boundary calculated in 3D is no good. Instead
0195 %   it needs to be fixed using find_boundary, later
0196 function [ELEM, NODE, BDY, elec_nodes] = mk_3D_model( ...
0197      elem0, node0, niveaux, bdy, elec_nodes0, node_order );
0198 
0199   elem0= node_reorder( elem0, node_order);
0200 
0201   d= size(elem0,1);       %dimentions+1
0202   n= size(node0,2);       %NODEs
0203   e= size(elem0,2);       %ELEMents
0204 
0205 %                   D     U
0206   elem_odd= [elem0([3,2,1,1],:), ... % 1 up 1 2 3 down
0207              elem0([3,2,2,1],:), ... % 1 2 up 2 3 down 
0208              elem0([3,3,2,1],:)];    % 1 2 3 up 3 down
0209   elem_even=[elem0([1,2,3,3],:), ... % 3 up 1 2 3 down
0210              elem0([1,2,2,3],:), ... % 3 2 up 2 1 down 
0211              elem0([1,1,2,3],:)];    % 3 2 1 up 1 down
0212 
0213   NODE= [node0; niveaux(1)*ones(1,n) ];
0214   ELEM= [];
0215   bl= size(bdy,2);
0216 % Interlaced bdy idx
0217 
0218   bdy_order =node_order(bdy);
0219   bdy_up= find(bdy_order>[1;1]*min(bdy_order));
0220   bdy_dn= find(bdy_order<[1;1]*max(bdy_order));
0221   
0222   bdy_odd = [bdy; bdy(bdy_up')];
0223   bdy_even= [bdy; bdy(bdy_dn')];
0224   BDY = [];
0225  
0226   ln= length(niveaux);
0227   for k=2:ln
0228     NODE=[NODE  [node0; niveaux(k)*ones(1,n)] ];
0229     if rem(k,2)==1
0230         elem= elem_odd;
0231         bdy_e0= bdy_even;
0232         bdy_e1= bdy_odd;
0233     else
0234         elem= elem_even;
0235         bdy_e1= bdy_even;
0236         bdy_e0= bdy_odd;
0237     end
0238     el_add = (k-2)*n+[[zeros(3,e);n*ones(1,e)], ...
0239                       [zeros(2,e);n*ones(2,e)], ...
0240                       [zeros(1,e);n*ones(3,e)]];
0241     ELEM= [ELEM,elem + el_add];
0242     BDY= [BDY, bdy_e0+(k-2)*n+[zeros(2,bl);n*ones(1,bl)], ...
0243                bdy_e1+(k-2)*n+[n*ones(2,bl);zeros(1,bl)] ];
0244   end %for k
0245 
0246   % Now add top and bottom boundary
0247   BDY= [elem0, BDY, elem0+n*(ln-1) ];
0248 
0249   % elec_nodes is all nodes for all layers
0250   elec_nodes= ones(ln,1) * elec_nodes0 + ...
0251               (0:ln-1)'  * n*ones(1, length(elec_nodes0) );
0252 
0253 
0254 %param.electrode = mk_electrodes( elec_nodes );
0255 % Create the electrode structure from elec_nodes
0256 % Currently implements point electrodes with
0257 %   contact impedance of near zero.
0258 function elec_struct = mk_electrodes( elec_nodes)
0259    for i= 1:length( elec_nodes )
0260       elec_struct(i).nodes     = elec_nodes(i);
0261       elec_struct(i).z_contact = 0.001; % corresponds to 1 ohm
0262    end
0263    % Need to do it this way to be compatible accross versions
0264    if ~exist('elec_struct');
0265        elec_struct= [];
0266    end
0267 
0268 function elem=  node_reorder( elem0, node_order);
0269   e= size(elem0,2);       %ELEMents
0270 
0271   no_test=  node_order(elem0);
0272   no_test=  (0:e-1)'*[3,3,3]+no_test';
0273   elem=     elem0(no_test');
0274 
0275   no_test = node_order(elem);
0276   ok= ~norm(no_test - [1;2;3]*ones(1,e));
0277 
0278   if ~ok; error('test_node_order fails - cant do 3D meshes'); end
0279 
0280 function do_unit_test
0281   subplot(3,3,1)
0282   mdl= mk_circ_tank(2, [], 2 );
0283   show_fem(mdl, [0,1,1]);
0284   unit_test_cmp('2D mdl', length(mdl.elems), 16);
0285 
0286   subplot(3,3,2)
0287   mdl= mk_circ_tank(4, [], 16 );
0288   show_fem(mdl);
0289   unit_test_cmp('2D mdl', length(mdl.nodes), 41);
0290 
0291   subplot(3,3,3)
0292   mdl= mk_circ_tank(4, linspace(-1,1,3), {'planes', 8, 2} );
0293   show_fem(mdl);
0294   unit_test_cmp('2D mdl', length(mdl.elems), 384);
0295 
0296   subplot(3,3,4)
0297   mdl= mk_circ_tank(4, linspace(-1,1,5), {'zigzag', 4, [2,4]} );
0298   show_fem(mdl);
0299   unit_test_cmp('2D mdl', length(mdl.elems), 768);
0300 
0301   try 
0302      mdl=  mk_circ_tank(2, [], 18 );  error
0303      unit_test_cmp('test for error', 1,0);
0304   catch
0305      unit_test_cmp('test for error', 1,1);
0306   end
0307 
0308   subplot(3,3,5)
0309   mdl= mk_circ_tank(4, [], 0);
0310   show_fem(mdl);
0311   title 'no electodes'
0312 
0313   subplot(3,3,6)
0314   mdl= mk_circ_tank(4, linspace(-1,1,5), 0);
0315   show_fem(mdl);
0316   title 'no electodes'