get_elem_volume

 GET_ELEM_VOLUME: VOL = get_elem_volume(fwd_model, map_node )
 Calculate volume (or area) of each element in model

 If the model has a 'coarse2fine' element, then the
 returned VOL applies to the coarse matrix (unless map_node <0)

 if map_node < 0 or map_node=='no_c2f', do not apply coarse2fine (if it exists)

 if abs(map_node) == 1, then calculated volumes are the volume fraction for each node
 BUGS: can't currently apply coarse2fine on map_node.

0001 function VOL = get_elem_volume( fwd_model, map_node )
0002 % GET_ELEM_VOLUME: VOL = get_elem_volume(fwd_model, map_node )
0003 % Calculate volume (or area) of each element in model
0004 %
0005 % If the model has a 'coarse2fine' element, then the
0006 % returned VOL applies to the coarse matrix (unless map_node <0)
0007 %
0008 % if map_node < 0 or map_node=='no_c2f', do not apply coarse2fine (if it exists)
0009 %
0010 % if abs(map_node) == 1, then calculated volumes are the volume fraction for each node
0011 % BUGS: can't currently apply coarse2fine on map_node.
0012 
0013 % (C) 2009 Andy Adler. License: GPL version 2 or version 3
0014 % $Id: get_elem_volume.m 6157 2021-11-02 13:02:02Z aadler $
0015 
0016 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0017 
0018 % If it looks like a fwd model
0019 if ~isfield(fwd_model,'type');
0020   if isfield(fwd_model,'nodes') && isfield(fwd_model,'elems')
0021      fwd_model.type = 'fwd_model';
0022   end
0023 end
0024    
0025 
0026 switch fwd_model.type
0027   case 'fwd_model'; % do nothing, we're ok
0028   case 'rec_model'; % do nothing, we're ok
0029   case 'inv_model'; fwd_model = fwd_model.fwd_model;
0030   case 'image';     fwd_model = fwd_model.fwd_model;
0031   otherwise;
0032      error('get_elem_volume: expecting fwd_model, got %s', ...
0033            fwd_model.type)
0034 end
0035 
0036 if nargin==1; map_node= 0; end
0037 if ischar(map_node) && strcmp(map_node,'no_c2f'); map_node=-inf; end
0038 
0039 % calculate element volume and surface area
0040 NODE = fwd_model.nodes';
0041 ELEM = fwd_model.elems';
0042 
0043 copt.cache_obj = {NODE, ELEM,map_node};
0044 copt.fstr = 'elem_volume';
0045 copt.boost_priority = -4;
0046 
0047 VOL = eidors_cache(@calc_volume, {NODE, ELEM}, copt);
0048 if isfield(fwd_model,'coarse2fine') && map_node>=0
0049    VOL= fwd_model.coarse2fine' * VOL;
0050 end
0051 
0052 % Calculate the mapping of each element onto the associated node
0053 % Map(i,j) = 1/Ne if elem j has node i
0054 if abs(map_node)==1
0055    [d,e]= size(ELEM);
0056    map = sparse( ELEM, ones(d,1)*(1:e), 1/d, size(NODE,2),e);
0057    VOL = map * VOL;
0058 end
0059 
0060 function VOL = calc_volume(NODE,ELEM)
0061     [d,e]= size(ELEM);
0062 
0063     VOL=zeros(e,1);
0064     ones_d = ones(1,d);
0065     d1fac = prod( 1:d-1 );
0066     if d > size(NODE,1)
0067         for i=1:e
0068             this_elem = NODE(:,ELEM(:,i));
0069             VOL(i)= abs(det([ones_d;this_elem])) / d1fac;
0070         end
0071     elseif d == 3 % 3D nodes in 2D mesh
0072         for i=1:e
0073             this_elem = NODE(:,ELEM(:,i));
0074             d12= det([ones_d;this_elem([1,2],:)])^2;
0075             d13= det([ones_d;this_elem([1,3],:)])^2;
0076             d23= det([ones_d;this_elem([2,3],:)])^2;
0077             VOL(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0078         end
0079     elseif d == 2 % 3D nodes in 1D mesh (ie resistor mesh)
0080         for i=1:e
0081             this_elem = NODE(:,ELEM(:,i));
0082             d12= det([ones_d;this_elem([1],:)])^2;
0083             d13= det([ones_d;this_elem([2],:)])^2;
0084             d23= det([ones_d;this_elem([3],:)])^2;
0085             VOL(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0086         end
0087     else
0088         error('mesh size not understood when calculating volumes')
0089     end
0090     
0091 
0092 
0093 function do_unit_test
0094   imdl = mk_common_model('a2c2',8);
0095   tt = 0.03125*ones(4,1);
0096   out = get_elem_volume(imdl.fwd_model);
0097   unit_test_cmp('fmdl:',  out(1:4), tt, 1e-10);
0098   out = get_elem_volume(imdl);
0099   unit_test_cmp('imdl:',  out(1:4), tt, 1e-10);
0100   out = get_elem_volume(mk_image(imdl,1));
0101   unit_test_cmp('image:', out(1:4), tt, 1e-10);
0102 
0103   fmdl = imdl.fwd_model;
0104   unit_test_cmp('fmdl (size 1):',  size(out,1), [num_elems(fmdl)]);
0105   out0= get_elem_volume(fmdl,0);
0106   unit_test_cmp('fmdl (size 1):',  size(out0,1), [num_elems(fmdl)]);
0107 
0108   out1= get_elem_volume(fmdl,1);
0109   unit_test_cmp('fmdl (size 1):',  size(out1,1),[num_nodes(fmdl)]);
0110   unit_test_cmp('fmdl (nodes 1):',  sum(out), sum(out1), 1e-10);
0111   ttn(1) = 0.041666666666667; ttn(2:5) = 0.088388347648318;
0112   unit_test_cmp('fmdl (nodes 2):',  out1(1:5), ttn(:), 1e-10);
0113 
0114   fmdl.coarse2fine = zeros(num_elems(fmdl),2);
0115   fmdl.coarse2fine(1:4,1) = 1;
0116   fmdl.coarse2fine(5:end,2) = 1;
0117   outc= get_elem_volume(fmdl,-2); % don't use
0118   unit_test_cmp('fmdl (c2f 1):',  out0, outc);
0119   outc= get_elem_volume(fmdl,0); % use c2f
0120   unit_test_cmp('fmdl (c2f 2):',  outc(1), sum(tt), 1e-10);
0121   unit_test_cmp('fmdl (c2f 3):',  sum(outc), sum(out0), 1e-10);
0122   unit_test_cmp('fmdl (c2f 4):',  size(outc), [2,1]);
0123 
0124   outc= get_elem_volume(fmdl,-1); % use c2f and map_node
0125   unit_test_cmp('fmdl (nodes 2):',  outc(1:5), ttn(:), 1e-10);