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 6164 2021-11-03 08:56:55Z bgrychtol $
0013 
0014 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 
0017 if nargin==1; % We need to calculate
0018    if isfield(img,'elem_data')
0019       img.fwd_solve.get_all_meas = 1;
0020       vh = fwd_solve(img);
0021       vv = vh.volt(:,1);
0022    elseif isfield(img,'node_data');
0023       vv = img.node_data(:,1);
0024       error('show_current: cannot interpolate conductivity onto elements (yet)');
0025    else
0026       error('show_current: one parameter provided, and cannot solve for voltages');
0027    end
0028 end 
0029 
0030 copt.fstr = 'calc_elem_current';
0031 elemcur = eidors_cache(@do_calc_elem_current, {img, vv}, copt);
0032 
0033 
0034 function elemcur = do_calc_elem_current(img, vv)
0035     dims = size(img.fwd_model.nodes,2);
0036     Nel = size(img.fwd_model.elems,1);
0037     elemcur = zeros(Nel,dims);
0038     % Calc field as I = sigma E
0039     %V1 = V0 + Ex*x1 + Ey*y1   [ 1 x1 y1 ] [ V0 ]
0040     %V2 = V0 + Ex*x2 + Ey*y2 = [ 1 x2 y2 ]*[ Ex ]
0041     %V3 = V0 + Ex*x3 + Ey*y    [ 1 x3 y3 ] [ Ey ]
0042     oo = ones(dims+1,1);
0043     for i=1:Nel
0044       idx = img.fwd_model.elems(i,:);
0045       nod = img.fwd_model.nodes(idx,:);
0046       if dims ==2
0047          VE  = ([oo, nod])\fix_dim(vv(idx));
0048       else
0049          VE  = ([oo, nod])\vv(idx);
0050       end
0051       elemcur(i,:) = img.elem_data(i,1)*VE(2:end)';
0052     %  elemcur(i+1,:) = (reshape(img.elem_data(i,1,:,:),dims,dims)*VE(2:end))';
0053     end
0054 
0055 % In case it is the wrong vector shape
0056 function vv = fix_dim(vv)
0057     if size(vv,1) == 1
0058         vv = vv';
0059     end
0060 
0061 function do_unit_test
0062    fmdl.nodes = [0,0;0,1;1,0;1,1];
0063    fmdl.elems = [1,2,3;2,3,4];
0064    fmdl.electrode(1).nodes = [1,2]; fmdl.electrode(1).z_contact = 0.01;
0065    fmdl.electrode(2).nodes = [3,4]; fmdl.electrode(2).z_contact = 0.01;
0066    fmdl.gnd_node = 1;
0067    fmdl.stimulation(1).stim_pattern = [1;-1];
0068    fmdl.stimulation(1).meas_pattern = [1,-1];
0069    fmdl.solve = @fwd_solve_1st_order;
0070    fmdl.system_mat = @system_mat_1st_order;
0071    fmdl.type = 'fwd_model';
0072    fmdl.normalize_measurements= 0;
0073    img = mk_image(fmdl,[1,1]); 
0074    img.fwd_solve.get_all_meas = 1;
0075 
0076    e_curr = calc_elem_current(img);
0077    unit_test_cmp('simple_mdl:', e_curr, [-1,0;-1,0],1e-12);
0078 
0079 
0080    imdl= mk_common_model('d2c2',8);
0081    img = calc_jacobian_bkgnd( imdl );
0082    img.fwd_solve.get_all_meas = 1;
0083    vh = fwd_solve(img);
0084    show_current(img, vh.volt(:,1));
0085    e_curr = calc_elem_current(img);
0086    unit_test_cmp('d2c2:', e_curr([1,10,100],:), ...
0087           [2.422593061890268, -0.920998260630422;
0088            2.887551382948032, -1.051869372020626;
0089            1.349507157073426, -0.842871247870084],1e-12);