Home > eidors > solvers > forward > system_mat_fields.m

system_mat_fields

PURPOSE ^

SYSTEM_MAT_FIELDS: fields (elem to nodes) fraction of system mat

SYNOPSIS ^

function FC= system_mat_fields( fwd_model )

DESCRIPTION ^

 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 ]

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5785 2018-05-21 01:05:57Z 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       [bdy_idx, bdy_area] = find_electrode_bdy( ...
0095           cem_boundary, fwd_model.nodes, eleci.nodes );
0096           % bdy_area is in [m] for 2D or [m^2] for 3D
0097 
0098       % For pt elec model, bdy_idx = [], so this doesn't run
0099       if isempty(bdy_idx); continue; end % no action for pt elecs
0100       i_cem = i_cem + 1;
0101 
0102       for j= 1:length(bdy_idx);
0103          bdy_nds= cem_boundary(bdy_idx(j),:);
0104 
0105          % 3D: [m^2]/[Ohm.m^2] = [S]
0106          % 2D: [m]  /[Ohm.m]   = [S]
0107          FFdata= [FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0108          FFiidx= [FFiidx; FFi_block' + sidx];
0109          FFjidx= [FFjidx; FFi_block  + cidx];
0110 
0111          CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0112          CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i_cem)*ones(1,d0)];
0113          CCdata= [CCdata; [1;-1]*ones(1,d0)];
0114          sidx = sidx + d0;
0115          cidx = cidx + d0;
0116       end
0117       
0118    end
0119 
0120 function do_unit_test
0121    imdl=  mk_common_model('a2c2',16);
0122    FC = system_mat_fields( imdl.fwd_model);
0123    unit_test_cmp('sys_mat1', size(FC), [128,41]);
0124    unit_test_cmp('sys_mat2', FC(1:2,:), [[0,-1,1,0;-2,1,1,0], zeros(2,37)]/2, 1e-14);
0125 
0126    % THis is a 45 degree rotation of the previous
0127    imdl=  mk_common_model('a2c0',16);
0128    FC2= system_mat_fields( imdl.fwd_model);
0129    M = sqrt(.5)*[1,-1;1,1];
0130    unit_test_cmp('sys_mat3', M*FC2(1:2,:), [[0,-1,1,0;-2,1,1,0], zeros(2,37)]/2, 1e-14);
0131 
0132    % Check we can handle partial CEM/pt electrodes
0133    imdl = mk_common_model('b2s',4); fmdl = imdl.fwd_model;
0134    fmdl.electrode([2,4])= []; % remove to simplify
0135    fmdl.stimulation = stim_meas_list([1,2,1,2]);
0136 
0137    FF = system_mat_fields(fmdl);
0138    unit_test_cmp('sys_mat_b2s_1a', size(FF), [256,81]);
0139    vh = fwd_solve( mk_image(fmdl,1) ); 
0140    unit_test_cmp('sys_mat_b2s_1b', vh.meas, 2.182865407049302, 1e-12);
0141 
0142    fmdl.electrode(2).nodes = [4:6];
0143    FF = system_mat_fields(fmdl);
0144    unit_test_cmp('sys_mat_b2s_2a', size(FF), [260,82]);
0145    vh = fwd_solve( mk_image(fmdl,1) ); 
0146    unit_test_cmp('sys_mat_b2s_2b', vh.meas, 1.807627806615849, 1e-12);
0147 
0148    fmdl.electrode(1).nodes = [76:78];
0149    FF = system_mat_fields(fmdl);
0150    unit_test_cmp('sys_mat_b2s_3a', size(FF), [264,83]);
0151    vh = fwd_solve( mk_image(fmdl,1) ); 
0152    unit_test_cmp('sys_mat_b2s_3b', vh.meas, 1.432226638247073, 1e-12);
0153 
0154    imdl=  mk_common_model('a2C2',4); fmdl=imdl.fwd_model;
0155    fmdl.nodes(1,:) = [];
0156    fmdl.gnd_node = 2;
0157    fmdl.elems(1:4,:) = [];
0158    fmdl.elems = fmdl.elems - 1;
0159    fmdl.boundary = find_boundary(fmdl);
0160    FC = system_mat_fields( fmdl);
0161    unit_test_cmp('sys_mat-bdyCEM-1', size(FC),[128,44]);
0162    unit_test_cmp('sys_mat-bdyCEM-2', FC(121:end,41:end), ...
0163              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0164 
0165 % Add a CEM via a boundary
0166    fmdl.electrode(5) = struct('nodes',1:4,'z_contact',.01);
0167    FC = system_mat_fields( fmdl);
0168    unit_test_cmp('sys_mat-ctrCEM-1', size(FC),[136,45]);
0169    unit_test_cmp('sys_mat-ctrCEM-2', FC(121:128,41:44), ...
0170              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0171    unit_test_cmp('sys_mat-ctrCEM-2', FC(129:end,end), ...
0172              -4.204482076268572,1e-12)
0173 
0174 % Add a CEM via a boundary
0175    imdl=  mk_common_model('a2C2',4); fmdl = imdl.fwd_model;
0176    fmdl.electrode(5) = struct('nodes',1:3,'z_contact',.01);
0177    fmdl.system_mat_fields.CEM_boundary = [1,2;2,3;3,1];
0178    FC = system_mat_fields( fmdl );
0179    unit_test_cmp('sys_mat-b2cCEM-1', size(FC),[142,46]);
0180    unit_test_cmp('sys_mat-b2cCEM-2', FC(129:136,42:45), ...
0181              -13.967473716321374*kron(eye(4),[1;1]),1e-12)
0182    unit_test_cmp('sys_mat-ctrCEM-3', FC([137:138,141:142],end), ...
0183              -3.535533905932737, 1e-12);
0184    unit_test_cmp('sys_mat-ctrCEM-4', FC([139:140],end), ...
0185              -4.204482076268572, 1e-12);

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005