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, 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, 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 5664 2017-12-12 15:14:20Z nolwenn85 $
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 
0038 % calculate element volume and surface area
0039 NODE = fwd_model.nodes';
0040 ELEM = fwd_model.elems';
0041 
0042 copt.cache_obj = {NODE, ELEM,map_node};
0043 copt.fstr = 'elem_volume';
0044 copt.boost_priority = -4;
0045 
0046 VOL = eidors_cache(@calc_volume, {NODE, ELEM}, copt);
0047 if isfield(fwd_model,'coarse2fine') && map_node>=0
0048    VOL= fwd_model.coarse2fine' * VOL;
0049 end
0050 
0051 % Calculate the mapping of each element onto the associated node
0052 % Map(i,j) = 1/Ne if elem j has node i
0053 if abs(map_node)==1
0054    [d,e]= size(ELEM);
0055    map = sparse( ELEM, ones(d,1)*(1:e), 1/d, size(NODE,2),e);
0056    VOL = map * VOL;
0057 end
0058 
0059 function VOL = calc_volume(NODE,ELEM)
0060     [d,e]= size(ELEM);
0061 
0062     VOL=zeros(e,1);
0063     ones_d = ones(1,d);
0064     d1fac = prod( 1:d-1 );
0065     if d > size(NODE,1)
0066         for i=1:e
0067             this_elem = NODE(:,ELEM(:,i));
0068             VOL(i)= abs(det([ones_d;this_elem])) / d1fac;
0069         end
0070     elseif d == 3 % 3D nodes in 2D mesh
0071         for i=1:e
0072             this_elem = NODE(:,ELEM(:,i));
0073             d12= det([ones_d;this_elem([1,2],:)])^2;
0074             d13= det([ones_d;this_elem([1,3],:)])^2;
0075             d23= det([ones_d;this_elem([2,3],:)])^2;
0076             VOL(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0077         end
0078     elseif d == 2 % 3D nodes in 1D mesh (ie resistor mesh)
0079         for i=1:e
0080             this_elem = NODE(:,ELEM(:,i));
0081             d12= det([ones_d;this_elem([1],:)])^2;
0082             d13= det([ones_d;this_elem([2],:)])^2;
0083             d23= det([ones_d;this_elem([3],:)])^2;
0084             VOL(i)= sqrt(d12 + d13 + d23 ) / d1fac;
0085         end
0086     else
0087         error('mesh size not understood when calculating volumes')
0088     end
0089     
0090 
0091 
0092 function do_unit_test
0093   imdl = mk_common_model('a2c2',8);
0094   tt = 0.03125*ones(4,1);
0095   out = get_elem_volume(imdl.fwd_model);
0096   unit_test_cmp('fmdl:',  out(1:4), tt, 1e-10);
0097   out = get_elem_volume(imdl);
0098   unit_test_cmp('imdl:',  out(1:4), tt, 1e-10);
0099   out = get_elem_volume(mk_image(imdl,1));
0100   unit_test_cmp('image:', out(1:4), tt, 1e-10);
0101 
0102   fmdl = imdl.fwd_model;
0103   unit_test_cmp('fmdl (size 1):',  size(out,1), [num_elems(fmdl)]);
0104   out0= get_elem_volume(fmdl,0);
0105   unit_test_cmp('fmdl (size 1):',  size(out0,1), [num_elems(fmdl)]);
0106 
0107   out1= get_elem_volume(fmdl,1);
0108   unit_test_cmp('fmdl (size 1):',  size(out1,1),[num_nodes(fmdl)]);
0109   unit_test_cmp('fmdl (nodes 1):',  sum(out), sum(out1), 1e-10);
0110   ttn(1) = 0.041666666666667; ttn(2:5) = 0.088388347648318;
0111   unit_test_cmp('fmdl (nodes 2):',  out1(1:5), ttn(:), 1e-10);
0112 
0113   fmdl.coarse2fine = zeros(num_elems(fmdl),2);
0114   fmdl.coarse2fine(1:4,1) = 1;
0115   fmdl.coarse2fine(5:end,2) = 1;
0116   outc= get_elem_volume(fmdl,-2); % don't use
0117   unit_test_cmp('fmdl (c2f 1):',  out0, outc);
0118   outc= get_elem_volume(fmdl,0); % use c2f
0119   unit_test_cmp('fmdl (c2f 2):',  outc(1), sum(tt), 1e-10);
0120   unit_test_cmp('fmdl (c2f 3):',  sum(outc), sum(out0), 1e-10);
0121   unit_test_cmp('fmdl (c2f 4):',  size(outc), [2,1]);
0122 
0123   outc= get_elem_volume(fmdl,-1); % use c2f and map_node
0124   unit_test_cmp('fmdl (nodes 2):',  outc(1:5), ttn(:), 1e-10);