system_mat_fields

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

 system_mat_fields detects whether an electrode is a CEM by checking
  1) whether it has more than 1 node, and 2) if it's on the boundary
 If you want an internal CEM that's not on the boundary, then it
 has to be specified by using
    fmdl.system_mat_fields.CEM_boundary = [extra boundary ]

0001 function FC= system_mat_fields( fwd_model )
0002 % SYSTEM_MAT_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 %   FC:        s_mat= C' * S * conduct * C = FC' * conduct * FC;
0008 %
0009 % system_mat_fields detects whether an electrode is a CEM by checking
0010 %  1) whether it has more than 1 node, and 2) if it's on the boundary
0011 % If you want an internal CEM that's not on the boundary, then it
0012 % has to be specified by using
0013 %    fmdl.system_mat_fields.CEM_boundary = [extra boundary ]
0014 
0015 % (C) 2008-2018 Andy Adler. License: GPL version 2 or version 3
0016 % $Id: system_mat_fields.m 5959 2019-06-05 02:08:59Z aadler $
0017 
0018 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0019 
0020 copt.cache_obj = mk_cache_obj(fwd_model);
0021 copt.fstr = 'system_mat_fields';
0022 copt.log_level = 4;
0023 FC= eidors_cache(@calc_system_mat_fields,{fwd_model},copt );
0024 
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= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
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    i_cem = 0; % Index into electrodes that are CEM (but not pt)
0086    cem_boundary = pp.boundary;
0087    if isfield(fwd_model,'system_mat_fields')
0088       cem_boundary = [cem_boundary;fwd_model.system_mat_fields.CEM_boundary];
0089    end
0090    for i= 1:pp.n_elec
0091       eleci = fwd_model.electrode(i);
0092       % contact impedance zc is in [Ohm.m] for 2D or [Ohm.m^2] for 3D
0093       zc=  eleci.z_contact;
0094       if isfield(eleci,'faces')  && ~isempty(eleci.faces)
0095          if ~isempty(eleci.nodes)
0096             eidors_msg('Warning: electrode %d has both faces and nodes',i);
0097          end
0098          bdy_nds= eleci.faces;
0099          [~, bdy_area] = find_electrode_bdy( ...
0100              bdy_nds, fwd_model.nodes,[]);
0101       else
0102          [bdy_idx, bdy_area] = find_electrode_bdy( ...
0103              cem_boundary, fwd_model.nodes, eleci.nodes );
0104              % bdy_area is in [m] for 2D or [m^2] for 3D
0105 
0106          % For pt elec model, bdy_idx = [], so this doesn't run
0107          if isempty(bdy_idx); continue; end % no action for pt elecs
0108 
0109          bdy_nds= cem_boundary(bdy_idx,:);
0110       end
0111       i_cem = i_cem + 1;
0112          % 3D: [m^2]/[Ohm.m^2] = [S]
0113          % 2D: [m]  /[Ohm.m]   = [S]
0114       for j= 1:size(bdy_nds,1);
0115          FFdata= [FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0116          FFiidx= [FFiidx; FFi_block' + sidx];
0117          FFjidx= [FFjidx; FFi_block  + cidx];
0118 
0119          CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0120          CCjidx= [CCjidx; bdy_nds(j,:) ; (pp.n_node+i_cem)*ones(1,d0)];
0121          CCdata= [CCdata; [1;-1]*ones(1,d0)];
0122          sidx = sidx + d0;
0123          cidx = cidx + d0;
0124       end
0125       
0126    end
0127 
0128 function do_unit_test
0129    imdl=  mk_common_model('a2c2',16);
0130    FC = system_mat_fields( imdl.fwd_model);
0131    unit_test_cmp('sys_mat1', size(FC), [128,41]);
0132    unit_test_cmp('sys_mat2', FC(1:2,:), [[0,-1,1,0;-2,1,1,0], zeros(2,37)]/2, 1e-14);
0133 
0134    % THis is a 45 degree rotation of the previous
0135    imdl=  mk_common_model('a2c0',16);
0136    FC2= system_mat_fields( imdl.fwd_model);
0137    M = sqrt(.5)*[1,-1;1,1];
0138    unit_test_cmp('sys_mat3', M*FC2(1:2,:), [[0,-1,1,0;-2,1,1,0], zeros(2,37)]/2, 1e-14);
0139 
0140    % Check we can handle partial CEM/pt electrodes
0141    imdl = mk_common_model('b2s',4); fmdl = imdl.fwd_model;
0142    fmdl.electrode([2,4])= []; % remove to simplify
0143    fmdl.stimulation = stim_meas_list([1,2,1,2]);
0144 
0145    FF = system_mat_fields(fmdl);
0146    unit_test_cmp('sys_mat_b2s_1a', size(FF), [256,81]);
0147    vh = fwd_solve( mk_image(fmdl,1) ); 
0148    unit_test_cmp('sys_mat_b2s_1b', vh.meas, 2.182865407049302, 1e-12);
0149 
0150    fmdl.electrode(2).nodes = [4:6];
0151    FF = system_mat_fields(fmdl);
0152    unit_test_cmp('sys_mat_b2s_2a', size(FF), [260,82]);
0153    vh = fwd_solve( mk_image(fmdl,1) ); 
0154    unit_test_cmp('sys_mat_b2s_2b', vh.meas, 1.807627806615849, 1e-12);
0155 
0156    fmdl.electrode(1).nodes = [76:78];
0157    FF = system_mat_fields(fmdl);
0158    unit_test_cmp('sys_mat_b2s_3a', size(FF), [264,83]);
0159    vh = fwd_solve( mk_image(fmdl,1) ); 
0160    unit_test_cmp('sys_mat_b2s_3b', vh.meas, 1.432226638247073, 1e-12);
0161 
0162    imdl=  mk_common_model('a2C2',4); fmdl=imdl.fwd_model;
0163    fmdl.nodes(1,:) = [];
0164    fmdl.gnd_node = 2;
0165    fmdl.elems(1:4,:) = [];
0166    fmdl.elems = fmdl.elems - 1;
0167    fmdl.boundary = find_boundary(fmdl);
0168    FC = system_mat_fields( fmdl);
0169    unit_test_cmp('sys_mat-bdyCEM-1', size(FC),[128,44]);
0170    unit_test_cmp('sys_mat-bdyCEM-2', FC(121:end,41:end), ...
0171              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0172 
0173 % Add a CEM via a boundary
0174    fmdl.electrode(5) = struct('nodes',1:4,'z_contact',.01);
0175    FC = system_mat_fields( fmdl);
0176    unit_test_cmp('sys_mat-ctrCEM-1', size(FC),[136,45]);
0177    unit_test_cmp('sys_mat-ctrCEM-2', FC(121:128,41:44), ...
0178              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0179    unit_test_cmp('sys_mat-ctrCEM-2', FC(129:end,end), ...
0180              -4.204482076268572,1e-12)
0181 
0182 % Add a CEM via a boundary
0183    Nelec = 4;
0184    fmdl= getfield( mk_common_model('a2C2',Nelec), 'fwd_model');
0185 
0186    FC = system_mat_fields( fmdl );
0187    unit_test_cmp('sys_mat-b2cCEM-0', Nelec, num_elecs(fmdl));
0188    fmdl.electrode(3).nodes = []; % no longer a node elec
0189    fmdl.electrode(3).faces = [1,2;2,3;3,1];
0190    FC = system_mat_fields( fmdl);
0191 spy(FC)
0192 %for i=129:148; disp([i,find(FC(i,:))]); end
0193 %full(FC(129:end,1:3))
0194 
0195 
0196    unit_test_cmp('sys_mat-b2cCEM-1', size(FC), ...
0197        [num_elems(fmdl)*(size(fmdl.elems,2)-1) + 3*2 + 2*3, ...
0198         num_nodes(fmdl)+num_elecs(fmdl)]);
0199 
0200    fmdl= getfield( mk_common_model('a2C2',Nelec), 'fwd_model');
0201    fmdl.electrode(5) = struct('nodes',1:3,'z_contact',.01);
0202    fmdl.system_mat_fields.CEM_boundary = [1,2;2,3;3,1];
0203    FC = system_mat_fields( fmdl );
0204 
0205    d1=num_elems(fmdl)*(size(fmdl.elems,2)-1) + ...
0206         4*2 + 2*3;
0207    d2= num_nodes(fmdl)+num_elecs(fmdl);
0208    unit_test_cmp('sys_mat-b2cCEM-1', size(FC),[d1,d2]);
0209 disp([size(FC),d1,d2])
0210    % Check this!!!! why are internal nodes doubled!!
0211    unit_test_cmp('sys_mat-b2cCEM-3', FC(129:136,42:45), ...
0212              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0213    unit_test_cmp('sys_mat-ctrCEM-4', FC([137:138,141:142],end), ...
0214              -3.535533905932737, 1e-12);
0215    unit_test_cmp('sys_mat-ctrCEM-5', FC([139:140],end), ...
0216              -4.204482076268572, 1e-12);