aa_system_mat_fields

 AA_SYSTEM_MAT_FIELDS: fields (elem to nodes) fraction of system mat
 FC= aa_system_mat_fields( fwd_model )
 input: 
   fwd_model = forward model
 output:
   FC:        s_mat= C' * S * conduct * C = FC' * conduct * FC;

0001 function FC= aa_system_mat_fields( fwd_model )
0002 % AA_SYSTEM_MAT_FIELDS: fields (elem to nodes) fraction of system mat
0003 % FC= aa_system_mat_fields( fwd_model )
0004 % input:
0005 %   fwd_model = forward model
0006 % output:
0007 %   FC:        s_mat= C' * S * conduct * C = FC' * conduct * FC;
0008 
0009 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0010 % $Id: aa_system_mat_fields.html 2819 2011-09-07 16:43:11Z aadler $
0011 
0012 cache_obj = mk_cache_obj(fwd_model);
0013 FC = eidors_obj('get-cache', cache_obj, 'aa_system_mat_fields');
0014 if ~isempty(FC)
0015    eidors_msg('aa_system_mat_fields: using cached value', 4);
0016    return
0017 end
0018 
0019 FC= calc_system_mat_fields( fwd_model );
0020 
0021 eidors_cache('boost_priority',1); % Moderate Priority boost
0022 eidors_obj('set-cache', cache_obj, 'aa_system_mat_fields', FC);
0023 eidors_msg('aa_system_mat_fields: setting cached value', 4);
0024 eidors_cache('boost_priority',-1);
0025 
0026 % only cache stuff which is really relevant here
0027 function cache_obj = mk_cache_obj(fwd_model);
0028    cache_obj.elems       = fwd_model.elems;
0029    cache_obj.nodes       = fwd_model.nodes;
0030    try
0031    cache_obj.electrode   = fwd_model.electrode; % if we have it
0032    end
0033    cache_obj.type        = 'fwd_model';
0034    cache_obj.name        = ''; % it has to have one
0035 
0036 function FC= calc_system_mat_fields( fwd_model );
0037    p= aa_fwd_parameters( fwd_model );
0038    d0= p.n_dims+0;
0039    d1= p.n_dims+1;
0040    e= p.n_elem;
0041    n= p.n_node;
0042 
0043    FFjidx= floor([0:d0*e-1]'/d0)*d1*ones(1,d1) + ones(d0*e,1)*(1:d1);
0044    FFiidx= [1:d0*e]'*ones(1,d1);
0045    FFdata= zeros(d0*e,d1);
0046    dfact = (d0-1)*d0;
0047    for j=1:e
0048      a=  inv([ ones(d1,1), p.NODE( :, p.ELEM(:,j) )' ]);
0049      idx= d0*(j-1)+1 : d0*j;
0050      FFdata(idx,1:d1)= a(2:d1,:)/ sqrt(dfact*abs(det(a)));
0051    end %for j=1:ELEMs
0052 
0053 if 0 % Not complete electrode model
0054    FF= sparse(FFiidx,FFjidx,FFdata);
0055    CC= sparse((1:d1*e),p.ELEM(:),ones(d1*e,1), d1*e, n);
0056 else
0057    [F2data,F2iidx,F2jidx, C2data,C2iidx,C2jidx] = ...
0058              compl_elec_mdl(fwd_model,p);
0059    FF= sparse([FFiidx(:); F2iidx(:)],...
0060               [FFjidx(:); F2jidx(:)],...
0061               [FFdata(:); F2data(:)]);
0062    
0063    CC= sparse([(1:d1*e)';    C2iidx(:)], ...
0064               [p.ELEM(:);   C2jidx(:)], ...
0065               [ones(d1*e,1); C2data(:)]);
0066 end
0067 
0068 FC= FF*CC;
0069 
0070 % Add parts for complete electrode model
0071 function [FFdata,FFiidx,FFjidx, CCdata,CCiidx,CCjidx] = ...
0072              compl_elec_mdl(fwd_model,pp);
0073    d0= pp.n_dims;
0074    FFdata= zeros(0,d0);
0075    FFd_block= sqrtm( ( ones(d0) + eye(d0) )/6/(d0-1) ); % 6 in 2D, 12 in 3D
0076    FFiidx= zeros(0,d0);
0077    FFjidx= zeros(0,d0);
0078    FFi_block= ones(d0,1)*(1:d0);
0079    CCdata= zeros(0,d0);
0080    CCiidx= zeros(0,d0);
0081    CCjidx= zeros(0,d0);
0082   
0083    sidx= d0*pp.n_elem;
0084    cidx= (d0+1)*pp.n_elem;
0085    for i= 1:pp.n_elec
0086       eleci = fwd_model.electrode(i);
0087       zc=  eleci.z_contact;
0088 %     ffb = find_bdy_idx( bdy, fwd_model.electrode(i).nodes);
0089       [bdy_idx, bdy_area] = find_electrode_bdy( ...
0090           pp.boundary, fwd_model.nodes, eleci.nodes );
0091 
0092       for j= 1:length(bdy_idx);
0093          bdy_nds= pp.boundary(bdy_idx(j),:);
0094 
0095          FFdata= [FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0096          FFiidx= [FFiidx; FFi_block' + sidx];
0097          FFjidx= [FFjidx; FFi_block  + cidx];
0098 
0099          CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0100          CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0101          CCdata= [CCdata; [1;-1]*ones(1,d0)];
0102          sidx = sidx + d0;
0103          cidx = cidx + d0;
0104       end
0105       
0106    end