ng_mk_common_model

 NG_MK_COMMON_MODEL: utility to create common models
 fmdl = ng_mk_cyl_model(mdl_type,
              mdl_shape, elec_pos, elec_shape);

 MDL_TYPE = 'CYL'
    mdl_shape = {height, [radius, [maxsz]]}
       if height = 0 -> calculate a 2D shape
       radius (OPT)  -> (default = 1)
       maxsz  (OPT)  -> max size of mesh elems (default = course mesh)
   
    ELECTRODE POSITIONS:
     elec_pos = [n_elecs_per_plane,z_planes] 
        OR
     elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)

 ELECTRODE SHAPES::
    elec_shape = [width,height, maxsz]  % Rectangular elecs
       OR
    elec_shape = [radius, 0, maxsz ]    % Circular elecs
       OR
    elec_shape = [0, 0, maxsz ]         % Point elecs

 MDL_TYPE = 'CIRC'
    Same as CYL interface, but uses netgen 2D interface. 
    For this reason cyl height must be zero
     elec_shape = [width, RFNUM]
    where RFNUM is an optional refinement parameter (larger = more) 

    For point electrodes, do not specify too small a maxh (ie < .1),
      otherwise netgen has a problem and crashes

 USAGE EXAMPLES:
 Simple 3D cylinder. Radius = 1. No electrodes
   fmdl= ng_mk_common_model('cyl',2,[8,.5],[0.05]);

0001 function fmdl = ng_mk_common_model(mdl_type, ...
0002                mdl_shape, elec_pos, elec_shape);
0003 % NG_MK_COMMON_MODEL: utility to create common models
0004 % fmdl = ng_mk_cyl_model(mdl_type,
0005 %              mdl_shape, elec_pos, elec_shape);
0006 %
0007 % MDL_TYPE = 'CYL'
0008 %    mdl_shape = {height, [radius, [maxsz]]}
0009 %       if height = 0 -> calculate a 2D shape
0010 %       radius (OPT)  -> (default = 1)
0011 %       maxsz  (OPT)  -> max size of mesh elems (default = course mesh)
0012 %
0013 %    ELECTRODE POSITIONS:
0014 %     elec_pos = [n_elecs_per_plane,z_planes]
0015 %        OR
0016 %     elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0017 %
0018 % ELECTRODE SHAPES::
0019 %    elec_shape = [width,height, maxsz]  % Rectangular elecs
0020 %       OR
0021 %    elec_shape = [radius, 0, maxsz ]    % Circular elecs
0022 %       OR
0023 %    elec_shape = [0, 0, maxsz ]         % Point elecs
0024 %
0025 % MDL_TYPE = 'CIRC'
0026 %    Same as CYL interface, but uses netgen 2D interface.
0027 %    For this reason cyl height must be zero
0028 %     elec_shape = [width, RFNUM]
0029 %    where RFNUM is an optional refinement parameter (larger = more)
0030 %
0031 %    For point electrodes, do not specify too small a maxh (ie < .1),
0032 %      otherwise netgen has a problem and crashes
0033 %
0034 % USAGE EXAMPLES:
0035 % Simple 3D cylinder. Radius = 1. No electrodes
0036 %   fmdl= ng_mk_common_model('cyl',2,[8,.5],[0.05]);
0037 
0038 % (C) Andy Adler, 2015.  Licenced under GPL v2 or v3
0039 % $Id: ng_mk_common_model.m 5112 2015-06-14 13:00:41Z aadler $
0040 
0041 if ischar(mdl_type) && strcmp(mdl_type,'UNIT_TEST'); do_unit_test; return; end
0042 
0043 if isstruct(mdl_type)
0044    error('No code yet to process options structs');
0045 else switch lower(mdl_type)
0046       case 'cyl'
0047          [body_geom, elec_geom, pp] = cyl_geom( mdl_shape, elec_pos, elec_shape);
0048          fmdl = ng_mk_geometric_models(body_geom, elec_geom);
0049          fmdl.body_geometry      = body_geom;
0050          fmdl.electrode_geometry = elec_geom;
0051       case 'circ'
0052          [shape, elec_pos, elec_shape] = circ_geom( mdl_shape, elec_pos, elec_shape);
0053          fmdl = ng_mk_2d_model(shape, elec_pos, elec_shape);
0054       otherwise
0055          error('mdl_type = (%s) not available (yet)', mdl_type);
0056 end; end
0057 
0058    if mdl_dim(fmdl)==3 && pp.is2D
0059       fmdl = mdl2d_from3d(fmdl);
0060    end
0061 
0062 % For 2D circle
0063 function [shape, elec_pos, elec_shape] = circ_geom( mdl_shape, elec_pos, elec_shape);
0064    if mdl_shape(1) ~=0; error('specifying "circ" implies height of 0'); end
0065    if     length(mdl_shape)==1
0066       rad = 1;
0067       maxh = 2*pi*1 / 16;
0068    elseif length(mdl_shape)==2
0069       rad = mdl_shape(2);
0070       maxh = 2*pi*rad / 16;
0071    else
0072       rad = mdl_shape(2);
0073       maxh= mdl_shape(3);
0074    end
0075 
0076 
0077    if     size(elec_pos,1) == 0;
0078       theta= [];
0079    elseif size(elec_pos,1) == 1;
0080       theta = linspace( 0, 2*pi, elec_pos(1)+1)'; theta(end)=[];
0081    else
0082       theta = pi/180*elec_pos(:,1);
0083    end
0084 
0085    elec_pos= rad*[sin(theta), cos(theta)];
0086 
0087    % Point electrodes need refinement specified
0088    if length(elec_shape) == 1 && elec_shape(1) == 0
0089       elec_shape(2) = 10; % recommended default from 2d code
0090    end
0091     
0092   % If the boundary point is near the electrode, then unnecessary
0093   % refinement will occur. This could be avoided by placing the
0094   % boundary nodes carefully, but that introduces edge cases.
0095    
0096    theta_pts = ceil(2*pi*rad/maxh) - length(theta);
0097    if theta_pts > 0 % if we have extra points to distribute
0098       theta_  = [theta; 2*pi+theta(1)];
0099       per_interval = round( theta_pts * diff(theta_)/2/pi );
0100       for i=1:length(per_interval)
0101          new_pts = linspace( theta_(i), theta_(i+1), per_interval(i)+2);
0102          theta = [theta; new_pts(2:end-1)']; 
0103       end
0104       theta = sort(theta);
0105 %  theta = linspace( 0, 2*pi, ceil(2*pi*rad/maxh))'; theta(end)=[];
0106    end
0107 
0108    % Must specify points starting like this: make then CCW
0109    shape = {rad*[sin(theta),-cos(theta)], maxh};
0110 
0111 
0112 function [body_geom, elec_geom, pp] = cyl_geom( cyl_shape, elec_pos, elec_shape);
0113    [body_geom, pp] = parse_shape(cyl_shape);
0114    elec_geom = parse_elecs(elec_pos, elec_shape, pp);
0115 
0116 
0117 function [body_geom,pp] = parse_shape(cyl_shape);
0118    tank_height = cyl_shape(1);
0119    tank_radius = 1;
0120    tank_maxh   = 0;
0121    is2D = 0;
0122 
0123    if length(cyl_shape)>1;
0124       tank_radius=cyl_shape(2);
0125    end
0126    if length(cyl_shape)>2; 
0127       tank_maxh  =cyl_shape(3);
0128    end
0129    if tank_height==0;
0130       is2D = 1;
0131 
0132       %Need some width to let netgen work, but not too much so
0133       % that it meshes the entire region
0134       tank_height = tank_radius/5; % initial extimate
0135       if tank_maxh>0
0136          tank_height = min(tank_height,2*tank_maxh);
0137       end
0138    end
0139 
0140    body_geom.cylinder.bottom_center = [0 0 0];
0141    body_geom.cylinder.top_center    = [0 0 tank_height];
0142    body_geom.cylinder.radius        = tank_radius;
0143    if tank_maxh > 0
0144       body_geom.max_edge_length     = tank_maxh;
0145    end
0146 
0147    pp.is2D = is2D; % put it here too.
0148    pp.tank_height= tank_height;
0149    pp.tank_radius= tank_radius;
0150    pp.body_geom  = body_geom;
0151 
0152 % ELECTRODE POSITIONS:
0153 %  elec_pos = [n_elecs_per_plane,z_planes]
0154 %     OR
0155 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0156 %
0157 % ELECTRODE SHAPES::
0158 %  elec_shape = [width,height, {maxsz}]  % Rectangular elecs
0159 %     OR
0160 %  elec_shape = [radius, {0, maxsz} ]  % Circular elecs
0161 %     maxsz  (OPT)  -> max size of mesh elems (default = courase mesh)
0162 %
0163 function elecs= parse_elecs(elec_pos, elec_shape, pp)
0164    is2D = pp.is2D;
0165    elec_geom = struct([]);
0166    hig = pp.tank_height;
0167    rad = pp.tank_radius;
0168    body_geom = pp.body_geom;
0169     
0170     n_elecs= size(elec_pos,1); 
0171     
0172     if n_elecs == 0
0173       elecs= {};
0174       return;
0175     end
0176    
0177    if is2D
0178       elec_pos(:,2) = hig/2;
0179    end
0180 
0181    % It never makes sense to specify only one elec
0182    % So elec_pos means the number of electrodes in this case
0183    if size(elec_pos,1) == 1
0184        % Parse elec_pos = [n_elecs_per_plane,z_planes]
0185       n_elecs= elec_pos(1); % per plane
0186       th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0187 
0188       on_elecs = ones(n_elecs, 1);
0189       el_th = []; 
0190       el_z  = []; 
0191       if length(elec_pos) == 1; % Height was forgotten, assume middle
0192          elec_pos(:,2) = hig/2;
0193       end
0194       for i=2:length(elec_pos)
0195         el_th = [el_th; th];
0196         el_z  = [el_z ; on_elecs*elec_pos(i)];
0197       end
0198 
0199    else
0200       el_th = elec_pos(:,1)*2*pi/360;
0201       el_z  = elec_pos(:,2);
0202    end
0203       
0204    n_elecs= size(el_z,1); 
0205 
0206    if size(elec_shape,1) == 1
0207       elec_shape = ones(n_elecs,1) * elec_shape;
0208    end
0209 
0210    elecs= {};
0211    for i= 1:n_elecs
0212      row = elec_shape(i,:); 
0213      elecs{i} = elec_spec( row, is2D, hig, rad, el_th(i), el_z(i) );
0214    end
0215 
0216 % use clockwise coordinate system
0217 
0218 function elec = elec_spec( row, is2D, hig, rad, el_th, el_z )
0219   xy_centre = rad*[sin(el_th),cos(el_th)];
0220   if     is2D
0221      if row(1) == 0;
0222         elec.point = [xy_centre, 0];
0223      else
0224         el_z = hig/2;
0225         row(2)=hig;
0226         elec.intersection.cylinder(1).top_center     = [0,0,2*row(1)];
0227         elec.intersection.cylinder(1).bottom_center  = [0,0,-2*row(1)];
0228         elec.intersection.cylinder(1).radius         = rad;
0229         elec.intersection.cylinder(1).complement_flag= 1;
0230         rad_dir = [sin(el_th);cos(el_th)];
0231         tan_dir = [0,1;-1,0]*rad_dir;
0232         elec.intersection.parallelepiped(1).vertex   = [0.97*xy_centre' - row(1)/2*tan_dir; el_z - row(2)/2];
0233         elec.intersection.parallelepiped(1).vector_a = [0; 0; 1];
0234         elec.intersection.parallelepiped(1).vector_b = [rad_dir; 0];
0235         elec.intersection.parallelepiped(1).vector_c = [tan_dir; 0];
0236         elec.intersection.parallelepiped(2).vertex   = [1.03*xy_centre' + row(1)/2*tan_dir; el_z + row(2)/2];
0237         elec.intersection.parallelepiped(2).vector_a =-[0; 0; 1];
0238         elec.intersection.parallelepiped(2).vector_b =-[rad_dir; 0];
0239         elec.intersection.parallelepiped(2).vector_c =-[tan_dir; 0];
0240 
0241         elec.enter_body_flag = 0;
0242      end
0243   else
0244      if row(1) == 0
0245         elec.point = [xy_centre, el_z];
0246      elseif length(row)<2 || row(2) == 0 % Circular electrodes
0247         elec.cylinder.top_center     = [1.03*xy_centre, el_z];
0248         elec.cylinder.bottom_center  = [0.97*xy_centre, el_z];
0249         elec.cylinder.radius         = row(1);
0250         elec.enter_body_flag = 0;
0251      elseif row(2)>0      % Rectangular electrodes
0252 % parallelepiped:     A parallelepiped is described by the following
0253 %                     subfields: vertex ([0; 0; 0]), vector_a ([1; 0; 0]),
0254 %                     vector_b ([0; 1; 0]), vector_c ([0; 0; 1]),
0255 %                     complement_flag (false).
0256         rad_dir = [sin(el_th);cos(el_th)];
0257         tan_dir = [0,1;-1,0]*rad_dir;
0258         elec.intersection.parallelepiped(1).vertex   = [0.97*xy_centre' - row(1)/2*tan_dir; el_z - row(2)/2];
0259         elec.intersection.parallelepiped(1).vector_a = [0; 0; 1];
0260         elec.intersection.parallelepiped(1).vector_b = [rad_dir; 0];
0261         elec.intersection.parallelepiped(1).vector_c = [tan_dir; 0];
0262         elec.intersection.parallelepiped(2).vertex   = [1.03*xy_centre' + row(1)/2*tan_dir; el_z + row(2)/2];
0263         elec.intersection.parallelepiped(2).vector_a =-[0; 0; 1];
0264         elec.intersection.parallelepiped(2).vector_b =-[rad_dir; 0];
0265         elec.intersection.parallelepiped(2).vector_c =-[tan_dir; 0];
0266      else
0267         error('negative electrode width not supported');
0268      end
0269   end
0270 
0271   if length(row)>=3 && row(3) > 0
0272      elec.max_edge_length     = row(3);
0273   else
0274      % Do nothing
0275   end
0276 
0277 
0278 function do_unit_test
0279   for tn = 1:do_test_number(0)
0280      fprintf('>>> ng_mk_cyl_models: TEST #%d\n',tn);
0281      fmdl= do_test_number(tn);
0282      subplot(3,3,rem(tn-1,9)+1)
0283      title(sprintf('test #%d',tn));
0284      show_fem(fmdl); drawnow
0285   end
0286 
0287 function fmdl= do_test_number(tn)
0288    eidors_msg('ng_mk_cyl_models: UNIT_TEST #%d',tn,1);
0289    switch tn
0290    case 1;
0291 % Simple 3D cylinder. Radius = 1. No electrodes
0292     fmdl= ng_mk_common_model('cyl',3,[0],[]); 
0293 
0294    case 2;
0295 % Simple 2D cylinder. Radius = 2. Set minsize to refine
0296     fmdl= ng_mk_common_model('cyl',[0,2,.2],[0],[]); 
0297 
0298    case 3;
0299 % 3D cylinder. Radius = 1. 2 planes of 8 elecs with radius 0.1
0300     fmdl= ng_mk_common_model('cyl',3,[8,1,2],[0.1]); 
0301 
0302    case 4;
0303 % 3D cylinder. Radius = 1. 6 circ elecs with elec refinement
0304     fmdl= ng_mk_common_model('cyl',3,[7,1],[0.2,0,0.05]); 
0305 
0306    case 5;
0307 % 3D cylinder. Radius = 1. 7 rect elecs with no refinement
0308     fmdl= ng_mk_common_model('cyl',3,[7,1],[0.2,0.3]); 
0309 
0310    case 6;
0311 % 2D cylinder. Radius = 1. 11 rect elecs with refinement
0312     fmdl= ng_mk_common_model('cyl',0,[11],[0.2,0,0.05]); 
0313 
0314    case 7;
0315 % 2D cylinder. Radius = 1.5. Refined(0.1). 11 elecs with refinement
0316     fmdl= ng_mk_common_model('cyl',[0,1,0.1],[11],[0.2,0,0.02]); 
0317 
0318    case 8;
0319 % 2D cylinder. elecs at 0, 90 and 120 degrees
0320     fmdl= ng_mk_common_model('cyl',0,[0;90;120],[0.2,0,0.03]); 
0321 
0322    case 9;
0323 % 2D cylinder. elecs at 0 (large,refined) and 120 (small) degrees
0324     fmdl= ng_mk_common_model('cyl',0,[0;120],[0.4,0,0.01;0.1,0,0.1]); 
0325 
0326    case 10;
0327 % 3D cylinder. elecs at 0, 30, 60, 90 in planes
0328     fmdl= ng_mk_common_model('cyl',3,[0,0.5;30,1;60,1.5;90,2.0],[0.2,0,0.1]); 
0329 
0330    case 11;
0331 % 3D cylinder. Various elecs at 0, 30, 60, 90 in planes
0332     el_pos = [0,0.5;30,1;60,1.5;90,2.0];
0333     el_sz  = [0.2,0,0.1;0.1,0,0.05;0.2,0.2,0.02;0.2,0.4,0.5];
0334     fmdl= ng_mk_common_model('cyl',3,el_pos,el_sz); 
0335 
0336 
0337    case 12;
0338     fmdl = ng_mk_common_model('circ', [0,1,.1], [4], [0.4,10]);
0339    case 13;
0340     fmdl = ng_mk_common_model('circ', 0, [9], [0.2,3]);
0341    case 14;
0342     fmdl = ng_mk_common_model('circ',[0,1,.10], [9], [0,1]);
0343    case 15;
0344     fmdl = ng_mk_common_model('circ',[0,1,.10], [9], [0]);
0345    case 16;
0346     fmdl = ng_mk_common_model('circ',[0,2,.05], [9], [0.1,100]);
0347 
0348    case 0; fmdl = 16; %%%% RETURN MAXIMUM
0349    otherwise;
0350      error('huh?')
0351    end
0352   
0353    if exist('img'); fmdl = img; end