system_mat_2p5d_fields

 SYSTEM_MAT_2P5D_FIELDS: fields (elem to nodes) fraction of system mat
 FC= system_mat_fields( fwd_model )
 input: 
   fwd_model = forward model
 output:
   FT:        s_mat = C' * T * conduct * C = FT' * conduct * FT;
   where FC'*S*FC + k^2*FT'*S*FT, with k=real#, forms part of the Fourier integral

0001 function FT= system_mat_2p5d_fields( fwd_model )
0002 % SYSTEM_MAT_2P5D_FIELDS: fields (elem to nodes) fraction of system mat
0003 % FC= system_mat_fields( fwd_model )
0004 % input:
0005 %   fwd_model = forward model
0006 % output:
0007 %   FT:        s_mat = C' * T * conduct * C = FT' * conduct * FT;
0008 %   where FC'*S*FC + k^2*FT'*S*FT, with k=real#, forms part of the Fourier integral
0009 
0010 % (C) 2008, 2015, 2016 Andy Adler, Alistair Boyle. License: GPL version 2 or version 3
0011 % $Id: system_mat_2p5d_fields.m 5424 2017-04-25 17:45:19Z aadler $
0012 
0013 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0014 
0015 copt.cache_obj = mk_cache_obj(fwd_model);
0016 copt.fstr = 'system_mat_2p5d_fields';
0017 copt.log_level = 4;
0018 FT = eidors_cache(@calc_system_mat_2p5d_fields,{fwd_model},copt );
0019 
0020 % only cache stuff which is really relevant here
0021 function cache_obj = mk_cache_obj(fwd_model)
0022    cache_obj.elems       = fwd_model.elems;
0023    cache_obj.nodes       = fwd_model.nodes;
0024    try
0025    cache_obj.electrode   = fwd_model.electrode; % if we have it
0026    end
0027    cache_obj.type        = 'fwd_model';
0028    cache_obj.name        = ''; % it has to have one
0029 
0030 function FT = calc_system_mat_2p5d_fields( fwd_model )
0031    p= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0032    d0= p.n_dims+0;
0033    d1= p.n_dims+1;
0034    e= p.n_elem;
0035    n= p.n_node;
0036    assert(d0 == 2, 'only supports 2D meshes for 2.5D system matrices');
0037 
0038    dfact = (d0-1)*d0;
0039    SSe2 = sqrtm((ones(3) + eye(3)) / 12);
0040    area=zeros(e,1);
0041    for j=1:e
0042       a = inv([ ones(d1,1), p.NODE( :, p.ELEM(:,j) )' ]);
0043       area(j) = 1 / sqrt(dfact*abs(det(a)));
0044    end %for j=1:ELEMs
0045 
0046    % NOTE: CEM electrodes are taken care of by the system_mat_fields produced
0047    % FC matrix where we sum FC'*S*FC + k^2*FT'*S*FT for 2.5D solutions, so we
0048    % just need to have the correctly sized FT matrix, so that the FC and FT
0049    % products may be added directly
0050    n_cem = 0; % count CEM electrodes, add 1 column & row per CEM (all zeros)
0051    for i=1:length(fwd_model.electrode)
0052       n_cem = n_cem + (length(fwd_model.electrode(i).nodes) > 1);
0053    end
0054    [FFiidx,FFjidx,FFdata] = find(kron(spdiag(area),SSe2));
0055    FF= sparse(FFiidx(:), FFjidx(:), FFdata(:));
0056    CC= sparse((1:d1*e)', p.ELEM(:), ones(d1*e,1), d1*e, n+n_cem);
0057 
0058    FT= FF*CC;
0059 
0060 function do_unit_test
0061    imdl=  mk_common_model('a2c2',16);
0062    FT = system_mat_2p5d_fields( imdl.fwd_model);
0063    unit_test_cmp('sys_mat1 pem', size(FT), [192,41]);
0064    ss = 58.787753826797278;
0065    unit_test_cmp('sys_mat2 pem', FT(1:2,:), [[4,1,1,0;1,4,1,0], zeros(2,37)]/ss, 1e-14);
0066 
0067    imdl=  mk_geophysics_model('h2a', 16);
0068    FC = system_mat_fields( imdl.fwd_model);
0069    FT = system_mat_2p5d_fields( imdl.fwd_model);
0070    unit_test_cmp('sys_mat1 cem', size(FT, 2), size(FC, 2));
0071    unit_test_cmp('sys_mat2 cem', FT(:, end-15:end), FT(:, end-15:end)*0);