update_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;

0001 function FC1 = update_system_mat_fields( fwd_model0, fwd_model1 )
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 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0010 % $Id: update_system_mat_fields.m 5424 2017-04-25 17:45:19Z aadler $
0011 
0012 if ischar(fwd_model0) && strcmp(fwd_model0,'UNIT_TEST'); do_unit_test; return; end
0013 
0014 copt.cache_obj = mk_cache_obj(fwd_model0);
0015 copt.fstr = 'system_mat_fields';
0016 copt.log_level = 4;
0017 FC0 = eidors_cache(@system_mat_fields,{fwd_model0},copt ); % this should be a cache hit
0018 
0019 t = tic;
0020 copt.cache_obj = mk_cache_obj(fwd_model1);
0021 copt.fstr = 'update_system_mat_fields';
0022 dFC = eidors_cache(@calc_update_system_mat_fields,{fwd_model0, fwd_model1},copt );
0023 
0024 % updated system_mat_fields
0025 FC1 = FC0 + dFC;
0026 t = toc(t);
0027 
0028 % now we fake out the caching system by telling it that we are doing
0029 % system_mat_fields(fwd_model1) with the new model, calls to system_mat_fields
0030 % will get the cached value from this function
0031 fstr = 'system_mat_fields';
0032 cache_obj = mk_cache_obj(fwd_model1);
0033 eidors_obj('set-cache', cache_obj, fstr, {FC1}, t);
0034 
0035 
0036 % only cache stuff which is really relevant here
0037 function cache_obj = mk_cache_obj(fwd_model)
0038    cache_obj.elems       = fwd_model.elems;
0039    cache_obj.nodes       = fwd_model.nodes;
0040    try
0041    cache_obj.electrode   = fwd_model.electrode; % if we have it
0042    end
0043    cache_obj.type        = 'fwd_model';
0044    cache_obj.name        = ''; % it has to have one
0045 
0046 function dFC= calc_update_system_mat_fields( fwd_model0, fwd_model1 )
0047    p0= fwd_model_parameters( fwd_model0, 'skip_VOLUME' );
0048    p1= fwd_model_parameters( fwd_model1, 'skip_VOLUME' );
0049    d0= p0.n_dims+0;
0050    d1= p0.n_dims+1;
0051    e= p0.n_elem;
0052    n= p0.n_node;
0053 
0054    % find moved nodes, then the elements affected by those moved nodes
0055    [dn,~] = find(abs(fwd_model0.nodes - fwd_model1.nodes) > eps);
0056    dn = unique(dn);
0057    [de,~] = find(ismember(fwd_model1.elems,dn)); % our modified nodes touched these elements
0058    assert(all(all(fwd_model0.elems == fwd_model1.elems)), 'expected fmdl0 and fmdl1 to have the same element structure');
0059    assert(all(de <= e), 'invalid elem# found for delta nodes');
0060 
0061    FFjidx= floor([0:d0*e-1]'/d0)*d1*ones(1,d1) + ones(d0*e,1)*(1:d1);
0062    FFiidx= [1:d0*e]'*ones(1,d1);
0063    FFdata= zeros(d0*e,d1);
0064    dfact = (d0-1)*d0;
0065    for j=de(:)'
0066      a0=  inv([ ones(d1,1), p0.NODE( :, p0.ELEM(:,j) )' ]);
0067      a1=  inv([ ones(d1,1), p1.NODE( :, p0.ELEM(:,j) )' ]);
0068      idx= d0*(j-1)+1 : d0*j;
0069      FFdata(idx,1:d1)= a1(2:d1,:)/ sqrt(dfact*abs(det(a1))) - a0(2:d1,:)/ sqrt(dfact*abs(det(a0)));
0070    end %for j=1:ELEMs
0071 
0072 if 0 % Not complete electrode model
0073    FF= sparse(FFiidx,FFjidx,FFdata);
0074    CC= sparse((1:d1*e),p0.ELEM(:),ones(d1*e,1), d1*e, n);
0075 else
0076    % TODO this could be better: we could look to see which, if any electrodes
0077    % have been modified and only update those ones... currently the
0078    % implementation here is pretty brain dead on the assumption this part is
0079    % fast
0080    % TODO currently can't handle electrode node number changes
0081    [F2data0,F2iidx0,F2jidx0, C2data0,C2iidx0,C2jidx0] = compl_elec_mdl(fwd_model0,p0,dn);
0082    [F2data1,F2iidx1,F2jidx1, C2data1,C2iidx1,C2jidx1] = compl_elec_mdl(fwd_model1,p1,dn);
0083 
0084    FF= sparse([FFiidx(:);  F2iidx0(:);  F2iidx1(:)],...
0085               [FFjidx(:);  F2jidx0(:);  F2jidx1(:)],...
0086               [FFdata(:); -F2data0(:); +F2data1(:)]);
0087    
0088    CC= sparse([(1:d1*e)';     C2iidx0(:)], ...
0089               [p0.ELEM(:);    C2jidx0(:)], ...
0090               [ones(d1*e,1);  C2data0(:)]);
0091 end
0092 
0093 dFC= FF*CC;
0094 
0095 % Add parts for complete electrode model
0096 function [FFdata,FFiidx,FFjidx, CCdata,CCiidx,CCjidx] = ...
0097              compl_elec_mdl(fwd_model,pp,dn)
0098    d0= pp.n_dims;
0099    FFdata= zeros(0,d0);
0100    FFd_block= sqrtm( ( ones(d0) + eye(d0) )/6/(d0-1) ); % 6 in 2D, 12 in 3D
0101    FFiidx= zeros(0,d0);
0102    FFjidx= zeros(0,d0);
0103    FFi_block= ones(d0,1)*(1:d0);
0104    CCdata= zeros(0,d0);
0105    CCiidx= zeros(0,d0);
0106    CCjidx= zeros(0,d0);
0107   
0108    sidx= d0*pp.n_elem;
0109    cidx= (d0+1)*pp.n_elem;
0110    for i= 1:pp.n_elec
0111       eleci = fwd_model.electrode(i);
0112       % contact impedance zc is in [Ohm.m] for 2D or [Ohm.m^2] for 3D
0113       zc=  eleci.z_contact;
0114       if any(find(ismember(eleci.nodes,dn)))
0115          [bdy_idx, bdy_area] = find_electrode_bdy( ...
0116              pp.boundary, fwd_model.nodes, eleci.nodes );
0117              % bdy_area is in [m] for 2D or [m^2] for 3D
0118    
0119          for j= 1:length(bdy_idx);
0120             bdy_nds= pp.boundary(bdy_idx(j),:);
0121    
0122             % 3D: [m^2]/[Ohm.m^2] = [S]
0123             % 2D: [m]  /[Ohm.m]   = [S]
0124             FFdata= [FFdata; FFd_block * sqrt(bdy_area(j)/zc)];
0125             FFiidx= [FFiidx; FFi_block' + sidx];
0126             FFjidx= [FFjidx; FFi_block  + cidx];
0127    
0128             CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0129             CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0130             CCdata= [CCdata; [1;-1]*ones(1,d0)];
0131             sidx = sidx + d0;
0132             cidx = cidx + d0;
0133          end
0134       else
0135          [bdy_idx] = find_electrode_bdy( ...
0136              pp.boundary, fwd_model.nodes, eleci.nodes );
0137              % bdy_area is in [m] for 2D or [m^2] for 3D
0138    
0139          for j= 1:length(bdy_idx);
0140             bdy_nds= pp.boundary(bdy_idx(j),:);
0141    
0142             % 3D: [m^2]/[Ohm.m^2] = [S]
0143             % 2D: [m]  /[Ohm.m]   = [S]
0144             FFdata= [FFdata; FFd_block * 0];
0145             FFiidx= [FFiidx; FFi_block' + sidx];
0146             FFjidx= [FFjidx; FFi_block  + cidx];
0147    
0148             CCiidx= [CCiidx; FFi_block(1:2,:) + cidx];
0149             CCjidx= [CCjidx; bdy_nds ; (pp.n_node+i)*ones(1,d0)];
0150             CCdata= [CCdata; [1;-1]*ones(1,d0)];
0151             sidx = sidx + d0;
0152             cidx = cidx + d0;
0153          end
0154       end
0155       
0156    end
0157 
0158 function do_unit_test
0159    disp('running system_mat_fields UNIT_TEST');
0160    system_mat_fields('UNIT_TEST'); % we depend explicitly on system_mat_fields... check it works properly!
0161 
0162    eidors_cache('clear_name', 'system_mat_fields');
0163    eidors_cache('clear_name', 'update_system_mat_fields');
0164    disp('running update_system_mat_fields UNIT_TEST');
0165    imdl=  mk_geophysics_model('h2e',32);
0166    ndim = size(imdl.fwd_model.nodes,2);
0167    for i = 1:10
0168       fmdl0 = imdl.fwd_model;
0169       fmdl0.nodes(1,:) = fmdl0.nodes(1,:) + rand(1,ndim)*1e-10; % defeat cache
0170       t=tic; FC0 = system_mat_fields(fmdl0); t0(i)=toc(t);
0171    
0172       % perturb nodes
0173       fmdl1 = fmdl0;
0174       nn = fmdl1.electrode(1).nodes;
0175       nn = [nn(1); nn(end)];
0176       fmdl1.nodes(nn,:) = fmdl1.nodes(nn,:) + 1e-4+rand(2,ndim)*1e-10; % perturb
0177       t=tic; FC1a = update_system_mat_fields(fmdl0, fmdl1); t1(i)=toc(t);
0178       t=tic; FC1b = system_mat_fields(fmdl1); t2(i)=toc(t);
0179       eidors_cache('clear_name', 'system_mat_fields');
0180       t=tic; FC1c = system_mat_fields(fmdl1); t3(i)=toc(t);
0181    end
0182    fprintf('system_mat_fields(fmdl0) = %0.2f sec\n',mean(t0));
0183    fprintf('update_system_mat_fields(fmdl0,fmdl1) = %0.3f sec [faster?]\n',mean(t1));
0184    fprintf('system_mat_fields(fmdl1) = %0.3f sec [cache hit?]\n',mean(t2));
0185    fprintf('system_mat_fields(fmdl1+delta) = %0.3f sec [recalculate]\n',mean(t3));
0186    unit_test_cmp('delta FC is fast                  ',mean(t1) < mean(t0)/2,1);
0187    unit_test_cmp('cache trick is really fast        ',mean(t2./t0) < 0.015,1);
0188    unit_test_cmp('cache was cleared before recalc   ',mean(t3) > mean(t0)*0.9,1); % did actually clear cache
0189    unit_test_cmp('cache trick returns correct result',FC1a,FC1b);
0190    unit_test_cmp('delta FC gives same result        ',FC1a,FC1c,10*eps);
0191    fprintf('speed-up: %0.2f\n',mean(t0./t1));
0192    if(sum(sum((FC1a - FC1c).^2)) > 0)
0193       err = FC1a - FC1c
0194    end
0195 
0196    if 0
0197       disp('---- profiling ----');
0198       fmdl1.nodes(nn,:) = fmdl1.nodes(nn,:) + rand(2,ndim)*1e-3; % perturb
0199       t=tic; FC0 = system_mat_fields(fmdl0); t0=toc(t);
0200       profile clear;
0201       profile on;
0202       t = tic; FC1a = update_system_mat_fields(fmdl0, fmdl1); t1=toc(t);
0203       profile off;
0204       profview;
0205       profsave(profile('info'),'profile');
0206       fprintf('update_system_mat_fields(fmdl0,fmdl1) = %0.3f sec [profiled]\n',t1);
0207    end