calc_elem_current

 calc_elem_current: calculate current vector in each FEM element

 e_curr = calc_elem_current( img, vv )
   img -> img object 
   volt-> voltage on nodes if not specified, img is solved
      via fwd_solve
   e_curr = current in each element [N_elems x N_dims]

0001 function elemcur = calc_elem_current( img, vv )
0002 % calc_elem_current: calculate current vector in each FEM element
0003 %
0004 % e_curr = calc_elem_current( img, vv )
0005 %   img -> img object
0006 %   volt-> voltage on nodes if not specified, img is solved
0007 %      via fwd_solve
0008 %   e_curr = current in each element [N_elems x N_dims]
0009 
0010 
0011 % (C) 2017 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: calc_elem_current.m 5426 2017-04-26 00:31:55Z aadler $
0013 
0014 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 dims = size(img.fwd_model.nodes,2);
0017 
0018 if nargin==1; % We need to calculate
0019    if isfield(img,'elem_data')
0020       img.fwd_solve.get_all_meas = 1;
0021       vh = fwd_solve(img);
0022       vv = vh.volt(:,1);
0023    elseif isfield(img,'node_data');
0024       vv = img.node_data(:,1);
0025       error('show_current: cannot interpolate conductivity onto elements (yet)');
0026    else
0027       error('show_current: one parameter provided, and cannot solve for voltages');
0028    end
0029 end 
0030 
0031 Nel = size(img.fwd_model.elems,1);
0032 elemcur = zeros(Nel,dims);
0033 % Calc field as I = sigma E
0034 %V1 = V0 + Ex*x1 + Ey*y1   [ 1 x1 y1 ] [ V0 ]
0035 %V2 = V0 + Ex*x2 + Ey*y2 = [ 1 x2 y2 ]*[ Ex ]
0036 %V3 = V0 + Ex*x3 + Ey*y    [ 1 x3 y3 ] [ Ey ]
0037 oo = ones(dims+1,1);
0038 for i=1:Nel
0039   idx = img.fwd_model.elems(i,:);
0040   nod = img.fwd_model.nodes(idx,:);
0041   if dims ==2
0042      VE  = ([oo, nod])\fix_dim(vv(idx));
0043   else
0044      VE  = ([oo, nod])\vv(idx);
0045   end
0046   elemcur(i,:) = img.elem_data(i,1)*VE(2:end)';
0047 %  elemcur(i+1,:) = (reshape(img.elem_data(i,1,:,:),dims,dims)*VE(2:end))';
0048 end
0049 
0050 % In case it is the wrong vector shape
0051 function vv = fix_dim(vv)
0052     if size(vv,1) == 1
0053         vv = vv';
0054     end
0055 
0056 function do_unit_test
0057    fmdl.nodes = [0,0;0,1;1,0;1,1];
0058    fmdl.elems = [1,2,3;2,3,4];
0059    fmdl.electrode(1).nodes = [1,2]; fmdl.electrode(1).z_contact = 0.01;
0060    fmdl.electrode(2).nodes = [3,4]; fmdl.electrode(2).z_contact = 0.01;
0061    fmdl.gnd_node = 1;
0062    fmdl.stimulation(1).stim_pattern = [1;-1];
0063    fmdl.stimulation(1).meas_pattern = [1,-1];
0064    fmdl.solve = @fwd_solve_1st_order;
0065    fmdl.system_mat = @system_mat_1st_order;
0066    fmdl.type = 'fwd_model';
0067    fmdl.normalize_measurements= 0;
0068    img = mk_image(fmdl,[1,1]); 
0069    img.fwd_solve.get_all_meas = 1;
0070 
0071    e_curr = calc_elem_current(img);
0072    unit_test_cmp('simple_mdl:', e_curr, [-1,0;-1,0],1e-12);
0073 
0074 
0075    imdl= mk_common_model('d2c2',8);
0076    img = calc_jacobian_bkgnd( imdl );
0077    img.fwd_solve.get_all_meas = 1;
0078    vh = fwd_solve(img);
0079    show_current(img, vh.volt(:,1));
0080    e_curr = calc_elem_current(img);
0081    unit_test_cmp('d2c2:', e_curr([1,10,100],:), ...
0082           [2.422593061890268, -0.920998260630422;
0083            2.887551382948032, -1.051869372020626;
0084            1.349507157073426, -0.842871247870084],1e-12);