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 6481 2022-12-26 10:42:06Z 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    if length(unique(fwd_model.elems(:)))~=n
0043        warning('EIDORS:unused_nodes','Number of nodes inconsistent with Elements. Consider using "remove_unused_nodes"');
0044    end
0045 
0046    FFjidx= floor([0:d0*e-1]'/d0)*d1*ones(1,d1) + ones(d0*e,1)*(1:d1);
0047    FFiidx= [1:d0*e]'*ones(1,d1);
0048    FFdata = assemble_elements(d1,d0,p);
0049 
0050 if 0 % Not complete electrode model
0051    FF= sparse(FFiidx,FFjidx,FFdata);
0052    CC= sparse((1:d1*e),p.ELEM(:),ones(d1*e,1), d1*e, n);
0053 else
0054    [F2data,F2iidx,F2jidx, C2data,C2iidx,C2jidx] = ...
0055              compl_elec_mdl(fwd_model,p);
0056    FF= sparse([FFiidx(:); F2iidx(:)],...
0057               [FFjidx(:); F2jidx(:)],...
0058               [FFdata(:); F2data(:)]);
0059    
0060    CC= sparse([(1:d1*e)';    C2iidx(:)], ...
0061       [double(p.ELEM(:));    C2jidx(:)], ...
0062               [ones(d1*e,1); C2data(:)]);
0063 end
0064 
0065 FC= FF*CC;
0066 
0067 % Add parts for complete electrode model
0068 function [FFdata,FFiidx,FFjidx, CCdata,CCiidx,CCjidx] = ...
0069              compl_elec_mdl(fmdl,pp)
0070    d0= pp.n_dims;
0071    FFdata= zeros(0,d0);
0072    FFiidx= zeros(0,d0);
0073    FFjidx= zeros(0,d0);
0074    CCdata= zeros(0,d0);
0075    CCiidx= zeros(0,d0);
0076    CCjidx= zeros(0,d0);
0077   
0078    sidx= d0*pp.n_elem;
0079    cidx= (d0+1)*pp.n_elem;
0080    i_cem = 0; % Index into electrodes that are CEM (but not pt)
0081    pp.cem_boundary = pp.boundary;
0082    if isfield(fmdl,'system_mat_fields')
0083       pp.cem_boundary = [pp.cem_boundary;
0084           fmdl.system_mat_fields.CEM_boundary];
0085       pp.cem_boundary = unique(pp.cem_boundary, ...
0086           'rows','stable');
0087           % TODO, how does this get duplicated
0088    end
0089    for i= 1:pp.n_elec
0090       eleci = fmdl.electrode(i);
0091       % contact impedance zc is in [Ohm.m] for 2D or [Ohm.m^2] for 3D
0092       zc=  eleci.z_contact;
0093       [i_cem, sidx, cidx, blk]= compl_elec_mdl_i( ...
0094          eleci, zc, i_cem, sidx, cidx, fmdl, pp);
0095       FFdata= [FFdata;blk.FFdata];
0096       FFiidx= [FFiidx;blk.FFiidx];
0097       FFjidx= [FFjidx;blk.FFjidx];
0098       CCdata= [CCdata;blk.CCdata];
0099       CCiidx= [CCiidx;blk.CCiidx];
0100       CCjidx= [CCjidx;blk.CCjidx];
0101    end
0102 
0103 function [i_cem, sidx, cidx, blk]=compl_elec_mdl_i( ...
0104        eleci, zc, i_cem, sidx, cidx, fmdl, pp);
0105    d0= pp.n_dims;
0106    FFd_block= sqrtm( ( ones(d0) + eye(d0) )/6/(d0-1) ); % 6 in 2D, 12 in 3D
0107    FFi_block= ones(d0,1)*(1:d0);
0108    blk.FFdata= zeros(0,d0);
0109    blk.FFiidx= zeros(0,d0);
0110    blk.FFjidx= zeros(0,d0);
0111    blk.CCdata= zeros(0,d0);
0112    blk.CCiidx= zeros(0,d0);
0113    blk.CCjidx= zeros(0,d0);
0114 
0115 
0116    if isfield(eleci,'faces')  && ~isempty(eleci.faces)
0117       if ~isempty(eleci.nodes)
0118          eidors_msg('Warning: electrode %d has both faces and nodes',i);
0119       end
0120       bdy_nds= eleci.faces;
0121       [~, bdy_area] = find_electrode_bdy( ...
0122           bdy_nds, fmdl.nodes,[]);
0123    else
0124       [bdy_idx, bdy_area] = find_electrode_bdy( ...
0125         pp.cem_boundary, fmdl.nodes, eleci.nodes);
0126           % bdy_area is in [m] for 2D or [m^2] for 3D
0127 
0128       % For pt elec model, bdy_idx = [], so this doesn't run
0129       if isempty(bdy_idx); return; end % no action for pt elecs
0130 
0131       bdy_nds= pp.cem_boundary(bdy_idx,:);
0132    end
0133    i_cem = i_cem + 1;
0134       % 3D: [m^2]/[Ohm.m^2] = [S]
0135       % 2D: [m]  /[Ohm.m]   = [S]
0136 
0137 %    for j= 1:size(bdy_nds,1);
0138 %       blk.FFdata= [blk.FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0139 %       blk.FFiidx= [blk.FFiidx; FFi_block' + sidx];
0140 %       blk.FFjidx= [blk.FFjidx; FFi_block  + cidx];
0141 %
0142 %       blk.CCiidx= [blk.CCiidx; FFi_block(1:2,:) + cidx];
0143 %       blk.CCjidx= [blk.CCjidx; bdy_nds(j,:) ; (pp.n_node+i_cem)*ones(1,d0)];
0144 %       blk.CCdata= [blk.CCdata; [1;-1]*ones(1,d0)];
0145 %       sidx = sidx + d0;
0146 %       cidx = cidx + d0;
0147 %    end
0148 
0149    N = size(bdy_nds,1);
0150    blk.FFdata = kron(sqrt(bdy_area(:)/zc), FFd_block);
0151    blk.FFiidx = repmat(FFi_block',N,1) + kron(d0*(0:(N-1))',ones(d0)) + sidx;
0152    blk.FFjidx = repmat(FFi_block ,N,1) + kron(d0*(0:(N-1))',ones(d0)) + cidx;
0153    blk.CCiidx = repmat(FFi_block(1:2,:),N,1) + kron(d0*(0:(N-1))',ones(2,d0)) + cidx;
0154    blk.CCjidx = (pp.n_node+i_cem) * ones(2*N,d0);
0155    blk.CCjidx(1:2:end,:) = bdy_nds;
0156    blk.CCdata = ones(2*N, d0);
0157    blk.CCdata(2:2:end,:) = -1;
0158    sidx = sidx + N*d0;
0159    cidx = cidx + N*d0;
0160    
0161 
0162 function  FFdata = assemble_elements(d1,d0,p);
0163    dfact = (d0-1)*d0;
0164    if 0 % OLD CODE
0165        FFdata= zeros(d0*p.n_elem,d1);
0166        for j=1:p.n_elem;
0167          a=  inv([ ones(d1,1), p.NODE( :, p.ELEM(:,j) )' ]);
0168          idx= d0*(j-1)+1 : d0*j;
0169          FFdata(idx,1:d1)= a(2:d1,:)/ sqrt(dfact*abs(det(a)));
0170        end %for j=1:ELEMs
0171    end
0172    M3 = ones(d1,d1,p.n_elem);
0173    M3(2:end,:,:)= reshape(p.NODE(:,p.ELEM),d0,d1,[]);
0174 % To stabilize inverse, remove average
0175 % This is a multiple of row1 = 1
0176 % This affects row one only
0177    M3(2:end,:,:) = bsxfun(@minus, M3(2:end,:,:), ...
0178               mean(M3(2:end,:,:),2));
0179    [I,D] = vectorize_inverse(M3);
0180    Is = bsxfun(@times, I, sqrt(abs(D)/dfact));
0181    Is = permute(Is(:,2:d1,:),[2,3,1]);
0182    FFdata= reshape(Is,[],d1);
0183 
0184 
0185 function do_unit_test
0186    
0187    do_unit_test_2d_test
0188    do_unit_test_3d_test
0189 
0190 function do_unit_test_3d_test
0191    imdl= mk_common_model('n3r2',[16,2]); 
0192    fmdl= imdl.fwd_model;
0193    FC = system_mat_fields( fmdl );
0194    SZ(2) = num_nodes(fmdl) + num_elecs(fmdl);
0195    SZ(1) = 3*(num_elems(fmdl) + num_elecs(fmdl)*2);
0196      % Two faces per electrode
0197    unit_test_cmp('sys_mat3d', size(FC), SZ);
0198 
0199    E14 = [...
0200    0.329216512695190  -0.329216512695190  0                   0;
0201    0.032424996343701  -0.032424996343701 -0.506252451622855   0.506252451622855;
0202   -0.098764953808557                   0  0.098764953808557                   0;
0203    0.329216512695190  -0.329216512695190  0                   0];
0204    unit_test_cmp('sys_mat3dx1', FC(1:4,[1,33,64,65]), E14,1e-14);
0205    
0206 
0207 function do_unit_test_2d_test
0208 
0209    imdl=  mk_common_model('a2c2',16);
0210    FC = system_mat_fields( imdl.fwd_model);
0211    unit_test_cmp('sys_mat1', size(FC), [128,41]);
0212    unit_test_cmp('sys_mat2', FC(1:2,:), [[0,-1,1,0;-2,1,1,0], zeros(2,37)]/2, 1e-14);
0213 
0214    % THis is a 45 degree rotation of the previous
0215    imdl=  mk_common_model('a2c0',16);
0216    FC2= system_mat_fields( imdl.fwd_model);
0217    M = sqrt(.5)*[1,-1;1,1];
0218    unit_test_cmp('sys_mat3', M*FC2(1:2,:), [[0,-1,1,0;-2,1,1,0], zeros(2,37)]/2, 1e-14);
0219 
0220    % Check we can handle partial CEM/pt electrodes
0221    imdl = mk_common_model('b2s',4); fmdl = imdl.fwd_model;
0222    fmdl.electrode([2,4])= []; % remove to simplify
0223    fmdl.stimulation = stim_meas_list([1,2,1,2]);
0224 
0225    FF = system_mat_fields(fmdl);
0226    unit_test_cmp('sys_mat_b2s_1a', size(FF), [256,81]);
0227    vh = fwd_solve( mk_image(fmdl,1) ); 
0228    unit_test_cmp('sys_mat_b2s_1b', vh.meas, 2.182865407049302, 1e-12);
0229 
0230    fmdl.electrode(2).nodes = [4:6];
0231    FF = system_mat_fields(fmdl);
0232    unit_test_cmp('sys_mat_b2s_2a', size(FF), [260,82]);
0233    vh = fwd_solve( mk_image(fmdl,1) ); 
0234    unit_test_cmp('sys_mat_b2s_2b', vh.meas, 1.807627806615849, 1e-12);
0235 
0236    fmdl.electrode(1).nodes = [76:78];
0237    FF = system_mat_fields(fmdl);
0238    unit_test_cmp('sys_mat_b2s_3a', size(FF), [264,83]);
0239    vh = fwd_solve( mk_image(fmdl,1) ); 
0240    unit_test_cmp('sys_mat_b2s_3b', vh.meas, 1.432226638247073, 1e-12);
0241 
0242    imdl=  mk_common_model('a2C2',4);
0243    fmdl=rmfield(imdl.fwd_model,'stimulation');
0244    fmdl.nodes(1,:) = [];
0245    fmdl.gnd_node = 2;
0246    fmdl.elems(1:4,:) = [];
0247    fmdl.elems = fmdl.elems - 1;
0248    fmdl.boundary = find_boundary(fmdl);
0249    FC = system_mat_fields( fmdl);
0250    unit_test_cmp('sys_mat-bdyCEM-1', size(FC),[128,44]);
0251    unit_test_cmp('sys_mat-bdyCEM-2', FC(121:end,41:end), ...
0252              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0253 
0254 % Add a CEM via a boundary
0255    fmdl.electrode(5) = struct('nodes',1:4,'z_contact',.01);
0256    FC = system_mat_fields( fmdl);
0257    unit_test_cmp('sys_mat-ctrCEM-1', size(FC),[136,45]);
0258    unit_test_cmp('sys_mat-ctrCEM-2', FC(121:128,41:44), ...
0259              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0260    unit_test_cmp('sys_mat-ctrCEM-2', FC(129:end,end), ...
0261              -4.204482076268572,1e-12)
0262 
0263 % Add a CEM via a boundary
0264    Nelec = 4;
0265    imdl= mk_common_model('a2C2',Nelec);
0266    fmdl=rmfield(imdl.fwd_model,'stimulation');
0267 
0268    FC = system_mat_fields( fmdl );
0269    unit_test_cmp('sys_mat-b2cCEM-0', Nelec, num_elecs(fmdl));
0270    fmdl.electrode(3).nodes = []; % no longer a node elec
0271    fmdl.electrode(3).faces = [1,2;2,3;3,1];
0272    FC = system_mat_fields( fmdl);
0273 %spy(FC)
0274 %for i=129:148; disp([i,find(FC(i,:))]); end
0275 %full(FC(129:end,1:3))
0276 
0277 
0278    unit_test_cmp('sys_mat-b2cCEM-1', size(FC), ...
0279        [num_elems(fmdl)*(size(fmdl.elems,2)-1) + 3*2 + 2*3, ...
0280         num_nodes(fmdl)+num_elecs(fmdl)]);
0281 
0282    imdl= mk_common_model('a2C2',Nelec);
0283    fmdl= rmfield(imdl.fwd_model,'stimulation');
0284    fmdl.electrode(5) = struct('nodes',1:3,'z_contact',.01);
0285    fmdl.system_mat_fields.CEM_boundary = [1,2;2,3;3,1];
0286    FC = system_mat_fields( fmdl );
0287 
0288    d1=num_elems(fmdl)*(size(fmdl.elems,2)-1) + ...
0289         4*2 + 2*3;
0290    d2= num_nodes(fmdl)+num_elecs(fmdl);
0291    unit_test_cmp('sys_mat-b2cCEM-1', size(FC),[d1,d2]);
0292    unit_test_cmp('sys_mat-b2cCEM-3', FC(129:136,42:45), ...
0293              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0294    unit_test_cmp('sys_mat-ctrCEM-4', FC([137:138,141:142],end), ...
0295              -3.535533905932737, 1e-12);
0296    unit_test_cmp('sys_mat-ctrCEM-5', FC([139:140],end), ...
0297              -4.204482076268572, 1e-12);