fwd_solve_1st_order

PURPOSE ^

FWD_SOLVE_1ST_ORDER: data= fwd_solve_1st_order( img)

SYNOPSIS ^

function data =fwd_solve_1st_order(fwd_model, img)

DESCRIPTION ^

 FWD_SOLVE_1ST_ORDER: data= fwd_solve_1st_order( img)
 First order FEM forward solver
 Input:
    img       = image struct
 Output:
    data = measurements struct
 Options: (to return internal FEM information)
    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
    img.fwd_solve.get_elec_curr= 1 (data.elec_curr = current on electrodes)

 Model Reduction: use precomputed fields to reduce the size of
    the forward solution. Nodes which are 1) not used in the output
    (i.e. not electrodes) 2) all connected to the same conductivity via
    the c2f mapping are applicable.
 see: Model Reduction for FEM Forward Solutions, Adler & Lionheart, EIT2016

    img.fwd_model.model_reduction = @calc_model_reduction;
       where the functionputs a struct with fields: main_region, regions
       OR
    img.fwd_model.model_reduction.main_region = vector, and 
    img.fwd_model.model_reduction.regions = struct

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function data =fwd_solve_1st_order(fwd_model, img)
0002 % FWD_SOLVE_1ST_ORDER: data= fwd_solve_1st_order( img)
0003 % First order FEM forward solver
0004 % Input:
0005 %    img       = image struct
0006 % Output:
0007 %    data = measurements struct
0008 % Options: (to return internal FEM information)
0009 %    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
0010 %    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
0011 %    img.fwd_solve.get_elec_curr= 1 (data.elec_curr = current on electrodes)
0012 %
0013 % Model Reduction: use precomputed fields to reduce the size of
0014 %    the forward solution. Nodes which are 1) not used in the output
0015 %    (i.e. not electrodes) 2) all connected to the same conductivity via
0016 %    the c2f mapping are applicable.
0017 % see: Model Reduction for FEM Forward Solutions, Adler & Lionheart, EIT2016
0018 %
0019 %    img.fwd_model.model_reduction = @calc_model_reduction;
0020 %       where the functionputs a struct with fields: main_region, regions
0021 %       OR
0022 %    img.fwd_model.model_reduction.main_region = vector, and
0023 %    img.fwd_model.model_reduction.regions = struct
0024 
0025 % (C) 1995-2017 Andy Adler. License: GPL version 2 or version 3
0026 % $Id: fwd_solve_1st_order.m 5578 2017-06-21 02:52:30Z aadler $
0027 
0028 % correct input paralemeters if function was called with only img
0029 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0030 
0031 if nargin == 1
0032    img= fwd_model;
0033 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0034    warning('EIDORS:DeprecatedInterface', ...
0035       ['Calling FWD_SOLVE_1ST_ORDER with two arguments is deprecated and will cause' ...
0036        ' an error in a future version. First argument ignored.']);
0037 end
0038 fwd_model= img.fwd_model;
0039 
0040 img = data_mapper(img);
0041 if ~ismember(img.current_params, supported_params)
0042     error('EIDORS:PhysicsNotSupported', '%s does not support %s', ...
0043     'FWD_SOLVE_1ST_ORDER',img.current_params);
0044 end
0045 % all calcs use conductivity
0046 img = convert_img_units(img, 'conductivity');
0047 
0048 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0049 pp = set_gnd_node(fwd_model, pp);
0050 s_mat= calc_system_mat( img );
0051 
0052 if isfield(fwd_model,'model_reduction')
0053    [s_mat.E, main_idx, pp] = mdl_reduction(s_mat.E, ...
0054            img.fwd_model.model_reduction, img, pp );
0055 else
0056    pp.mr_mapper = 1:size(s_mat.E,1);
0057 end
0058 
0059 % Normally EIT uses current stimulation. In this case there is
0060 %  only a ground node, and this is the only dirichlet_nodes value.
0061 %  In that case length(dirichlet_nodes) is 1 and the loop runs once
0062 % If voltage stimulation is done, then we need to loop on the
0063 %  matrix to calculate faster.
0064 [dirichlet_nodes, dirichlet_values, neumann_nodes, has_gnd_node]= ...
0065          find_dirichlet_nodes( fwd_model, pp );
0066 
0067 v= full(horzcat(dirichlet_values{:})); % Pre fill in matrix
0068 for i=1:length(dirichlet_nodes)
0069    idx= 1:size(s_mat.E,1);
0070    idx( dirichlet_nodes{i} ) = [];
0071    if length(dirichlet_nodes) == 1; rhs = 1:size(pp.QQ,2);
0072    else                           ; rhs = i; end
0073    v(idx,rhs)= left_divide( s_mat.E(idx,idx), ...
0074              neumann_nodes{i}(idx,:) - s_mat.E(idx,:)*dirichlet_values{i});
0075 end
0076 
0077 % If model has a ground node (rather than voltage stim electrodes)
0078 if has_gnd_node
0079    Ignd = s_mat.E(dirichlet_nodes{1},:)*v;
0080    Irel = Ignd./sum(abs(pp.QQ)); % relative
0081    if max(abs(Irel))>1e-6
0082       warning('current flowing through ground node. Check stimulation pattern')
0083    end
0084 end
0085 
0086 % calc voltage on electrodes
0087 
0088 % This is horribly inefficient, override
0089 % v_els= pp.N2E * v;
0090 idx = find(any(pp.N2E));
0091 v_els= pp.N2E(:,idx) * v(idx,:);
0092 
0093 
0094 % create a data structure to return
0095 data.meas= meas_from_v_els(v_els, fwd_model.stimulation);
0096 data.time= NaN; % unknown
0097 data.name= 'solved by fwd_solve_1st_order';
0098 try; if img.fwd_solve.get_all_meas == 1
0099    outmap = pp.mr_mapper(1:pp.n_node);
0100    data.volt(outmap,:) = v(1:pp.n_node,:); % but not on CEM nodes
0101 end; end
0102 try; if img.fwd_solve.get_all_nodes== 1
0103    data.volt(pp.mr_mapper,:) = v;                % all, including CEM nodes
0104 end; end
0105 try; if img.fwd_solve.get_elec_curr== 1
0106 %  data.elec_curr = pp.N2E * s_mat.E * v;
0107    idx = find(any(pp.N2E));
0108    data.elec_curr = pp.N2E(:,idx) * s_mat.E(idx,:) * v;
0109 end; end
0110 
0111 
0112 % has_gnd_node = flag if the model has a gnd_node => can warn if current flows
0113 function [dirichlet_nodes, dirichlet_values, neumann_nodes, has_gnd_node]= ...
0114             find_dirichlet_nodes( fwd_model, pp );
0115    fnanQQ = isnan(pp.QQ);
0116    lstims = size(pp.QQ,2);
0117    if any(fnanQQ)
0118       has_gnd_node = 0; % no ground node is specified
0119       % Are all dirichlet_nodes the same
0120 
0121       % Don't use all on sparse, it will make them full
0122       if ~any(any(fnanQQ(:,1)*ones(1,lstims) - fnanQQ,2))
0123          dirichlet_nodes{1} = find(fnanQQ(:,1));
0124          dirichlet_values{1} = sparse(size(pp.N2E,2), size(fnanQQ,2));
0125          dirichlet_values{1}(fnanQQ) = pp.VV(fnanQQ);
0126          neumann_nodes{1} = pp.QQ;
0127          neumann_nodes{1}(fnanQQ) = 0;
0128       else % one at a time
0129          for i=1:size(fnanQQ,2)
0130             fnanQQi= fnanQQ(:,i);
0131             dirichlet_nodes{i} = find(fnanQQi);
0132             dirichlet_values{i} = sparse(size(pp.N2E,2), 1);
0133             dirichlet_values{i}(fnanQQi) = pp.VV(fnanQQi,i);
0134             neumann_nodes{i} = pp.QQ(:,i);
0135             neumann_nodes{i}(fnanQQi) = 0;
0136          end
0137       end
0138    elseif isfield(pp,'gnd_node')
0139       dirichlet_nodes{1} = pp.gnd_node;
0140       dirichlet_values{1} = sparse(size(pp.N2E,2), size(fnanQQ,2));
0141       neumann_nodes{1}   = pp.QQ;
0142       has_gnd_node= 1;
0143    else
0144       error('no required ground node on model');
0145    end
0146 
0147 function pp = set_gnd_node(fwd_model, pp);
0148    if isfield(fwd_model,'gnd_node');
0149       pp.gnd_node = fwd_model.gnd_node;
0150    else
0151       % try to find one in the model center
0152       ctr =  mean(fwd_model.nodes,1);
0153       d2  =  sum(bsxfun(@minus,fwd_model.nodes,ctr).^2,2);
0154       [~,pp.gnd_node] = min(d2);
0155       eidors_msg('Warning: no ground node found: choosing node %d',pp.gnd_node(1),1);
0156    end
0157 
0158 function vv = meas_from_v_els( v_els, stim)
0159    try
0160 % Was 1.82s
0161 %        % measured voltages from v
0162 %    %   vv = zeros( pp.n_meas, 1 );
0163 %        idx=0;
0164 %        for i=1:length(stim)
0165 %           meas_pat= stim(i).meas_pattern;
0166 %           n_meas  = size(meas_pat,1);
0167 %           vv( idx+(1:n_meas) ) = meas_pat*v_els(:,i);
0168 %           idx= idx+ n_meas;
0169 %        end
0170 
0171 % This code replaced the previous - Nov 4, 2013
0172 % Now 0.437s
0173 % Why is it faster??
0174 
0175        [n_elec,n_stim] = size(v_els);
0176 
0177        copt.cache_obj = {stim};
0178        copt.fstr = 'v2meas';
0179        copt.log_level = 4;
0180        v2meas = eidors_cache(@get_v2meas, {n_elec,n_stim,stim}, copt);
0181        vv = v2meas' * v_els(:);
0182    catch err
0183       if strcmp(err.identifier, 'MATLAB:innerdim');
0184           error(['measurement pattern not compatible with number' ...
0185                   'of electrodes for stimulation patter %d'],i);
0186       else
0187           rethrow(err);
0188       end
0189    end
0190 
0191    
0192 function v2meas = get_v2meas(n_elec,n_stim,stim)
0193     v2meas = sparse(n_elec*n_stim,0);
0194     for i=1:n_stim
0195         meas_pat= stim(i).meas_pattern;
0196         n_meas  = size(meas_pat,1);
0197         v2meas((i-1)*n_elec + 1: i*n_elec,end+(1:n_meas)) = meas_pat';
0198     end
0199 
0200 
0201 function [E, m_idx, pp] = mdl_reduction(E, mr, img, pp);
0202    % if mr is a string we assume it's a function name
0203    if isa(mr,'function_handle') || isstr(mr)
0204       mr = feval(mr,img.fwd_model);
0205    end
0206    % mr is now a struct with fields: main_region, regions
0207    m_idx = mr.main_region;
0208    E = E(m_idx, m_idx);
0209    for i=1:length(mr.regions)
0210       invEi=   mr.regions(i).invE;
0211 % FIXME:!!! data_mapper has done the c2f. But we don't want that here.
0212 %  kludge is to reach into the fine model field. This is only ok because
0213 %  model_reduction is only valid if one parameter describes each field
0214       field = mr.regions(i).field; 
0215       field = find(img.fwd_model.coarse2fine(:,field));
0216       field = field(1); % they're all the same - by def of model_reduction
0217       sigma = img.elem_data(field);
0218       E = E - sigma*invEi;
0219    end
0220 
0221    % Adjust the applied current and measurement matrices
0222    pp.QQ = pp.QQ(m_idx,:);
0223    pp.VV = pp.VV(m_idx,:);
0224    pp.N2E= pp.N2E(:,m_idx);
0225    pp.mr_mapper = cumsum(m_idx); %must be logical
0226    pp.gnd_node = pp.mr_mapper(pp.gnd_node);
0227    if pp.gnd_node==0
0228       error('model_reduction removes ground node');
0229    end
0230    
0231         
0232 function unit_test_voltage_stims;
0233    stim = zeros(16,1); volt=stim; stim([1,4]) = NaN; volt([1,4]) = [1,2];
0234    stimulv.stim_pattern = stim;
0235    stimulv.volt_pattern = volt;
0236    stimulv.meas_pattern = [1,0,0,-1,zeros(1,12)];
0237 
0238    img = mk_image( mk_common_model('a2c2',16),1);
0239    img.fwd_model.stimulation = stimulv;
0240    img.fwd_solve.get_all_nodes = 1;
0241    vh = fwd_solve_1st_order(img);
0242    unit_test_cmp('a2c2 Vstim #1', vh.meas, diff(volt(1:2)), 1e-14);
0243    unit_test_cmp('a2c2 Vstim #2', vh.volt(27+[1,4]), volt([1,4]), 1e-14);
0244    tst = [ ...
0245    1.503131926779798; 1.412534629974291; 1.529078332819747;
0246    1.354399248512161; 1.546241676995996];
0247    unit_test_cmp('a2c2 Vstim #3', vh.volt(1:5:25), tst, 1e-14);
0248 
0249    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt;
0250    imgn.calc_colours.clim = 1; subplot(221); show_fem(imgn,1);
0251 
0252    img = mk_image( mk_common_model('a2C2',16),1);
0253    img.fwd_model.stimulation = stimulv;
0254    img.fwd_solve.get_all_nodes = 1;
0255    vh = fwd_solve_1st_order(img);
0256    unit_test_cmp('a2C2 Vstim #1', vh.meas, diff(volt(1:2)), 1e-14);
0257    unit_test_cmp('a2C2 Vstim #2', vh.volt(num_nodes(img)+[1,4]), volt([1,4]), 1e-14);
0258    tst = [ ...
0259    1.499999999999998; 1.302478674263331; 1.609665333411830; ...
0260    1.215039511028270; 1.691145536046686];
0261    unit_test_cmp('a2C2 Vstim #3', vh.volt(1:5:25), tst, 1e-13);
0262 
0263    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt(1:num_nodes(img));
0264    imgn.calc_colours.clim = 1; subplot(222); show_fem(imgn,1);
0265 
0266    stim = zeros(16,1); volt=stim; stim([1,4]) = NaN; stim(8)=1; volt([1,4]) = [1,1];
0267    img.fwd_model.stimulation(2).stim_pattern = stim;
0268    img.fwd_model.stimulation(2).volt_pattern = volt;
0269    img.fwd_model.stimulation(2).meas_pattern = [1,-1,zeros(1,14)];
0270    vh = fwd_solve_1st_order(img);
0271    unit_test_cmp('a2C2 Vstim #4', vh.volt(num_nodes(img)+[1,4],1), [1;2], 1e-14);
0272    unit_test_cmp('a2C2 Vstim #5', vh.volt(1:5:25,1), tst, 1e-13);
0273    unit_test_cmp('a2C2 Vstim #6', vh.volt(num_nodes(img)+[1,4],2), [1;1], 1e-14);
0274    tst = [ 1.029942389400905; 1.024198991581187; ...
0275            1.048244746016660; 1.006551737030278; 1.057453501332724];
0276    unit_test_cmp('a2C2 Vstim #6', vh.volt(1:5:25,2), tst, 1e-14);
0277 
0278    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt(1:num_nodes(img),2);
0279    subplot(223); show_fem(imgn,1);
0280 
0281    stim = zeros(16,1); volt=stim; stim([3,6]) = NaN;  volt([3,6]) = [1,2];
0282    img.fwd_model.stimulation(3).stim_pattern = stim;
0283    img.fwd_model.stimulation(3).volt_pattern = volt;
0284    img.fwd_model.stimulation(3).meas_pattern = [1,-1,zeros(1,14)];
0285    vh = fwd_solve_1st_order(img);
0286 
0287    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt(1:num_nodes(img),3);
0288    imgn.calc_colours.clim = 1; subplot(224); show_fem(imgn,1);
0289 
0290    unit_test_cmp('a2C2 Vstim #7', vh.volt(num_nodes(img)+[3,6],3), [1;2], 1e-14);
0291 
0292 
0293 
0294 function do_unit_test
0295    unit_test_voltage_stims;
0296 
0297 
0298    img = mk_image( mk_common_model('b2c2',16),1);
0299    vh = fwd_solve_1st_order(img);
0300    tst = [ ...
0301     0.959567140078593; 0.422175237237900; 0.252450963869202; ...
0302     0.180376116490602; 0.143799778367518];
0303    unit_test_cmp('b2c2 TEST', vh.meas(1:5), tst, 1e-12);
0304 
0305    img.fwd_model = rmfield(img.fwd_model,'gnd_node');
0306    vh = fwd_solve_1st_order(img);
0307    unit_test_cmp('b2c2 gnd_node', vh.meas(1:5), tst, 1e-12);
0308 
0309    img.fwd_solve.get_elec_curr = 1;
0310    vh = fwd_solve_1st_order(img);
0311    pp = fwd_model_parameters( img.fwd_model); EC = pp.N2E*pp.QQ;
0312    unit_test_cmp('b2b2 (CEM) elec_curr', vh.elec_curr, EC, 1e-11);
0313 
0314    img.fwd_solve.get_all_meas = 1;
0315    vh = fwd_solve_1st_order(img);
0316     plot(vh.volt);
0317 
0318    img = mk_image( mk_common_model('b2C2',16),1);
0319    vh = fwd_solve_1st_order(img);
0320    tst = [ 0.385629619754662; 0.235061644846908; 0.172837756982388
0321            0.142197580506776; 0.126808900182258; 0.120605655110661];
0322    unit_test_cmp('b2C2 (CEM) TEST', vh.meas(15:20), tst, 1e-12);
0323 
0324    img.fwd_solve.get_elec_curr = 1;
0325    vh = fwd_solve_1st_order(img);
0326    pp = fwd_model_parameters( img.fwd_model); EC = pp.N2E*pp.QQ;
0327    unit_test_cmp('b2C2 (CEM) elec_curr', vh.elec_curr, EC, 1e-11);
0328 
0329    % bad stim patterns (flow through ground node)
0330    img.fwd_model.stimulation(1).stim_pattern(2) = 0;
0331    vh = fwd_solve_1st_order(img);
0332    lastw = lastwarn;
0333    unit_test_cmp('gnd_node warning', lastw, ...
0334     'current flowing through ground node. Check stimulation pattern');
0335 
0336 
0337    %2D resistor
0338    current = 4; measure=1;
0339    [R,img] = test_2d_resistor(current,measure);
0340    img.fwd_solve.get_all_nodes = 1;
0341    vs = fwd_solve_1st_order( img);
0342    va= measure*current*sum(R); % analytic
0343    unit_test_cmp('2D resistor test', va, vs.meas, 1e-12);
0344 
0345    unit_test_cmp('2D R z_contact', ...
0346                  [diff(vs.volt([13,1])), diff(vs.volt([14,12]))], ...
0347                  R(2)/2*current*[1,-1], 1e-12);
0348    unit_test_cmp('2D R voltages', vs.volt(1:3:10)-vs.volt(1), ...
0349                  R(1)*current*linspace(0,1,4)', 1e-12);
0350 
0351    %3D resistor
0352    [R,img] = test_3d_resistor(current,measure);
0353    img.fwd_solve.get_all_nodes = 1;
0354    vs = fwd_solve_1st_order( img);
0355    va= current*sum(R);
0356    unit_test_cmp('3D resistor test', va, vs.meas, 1e-10);
0357    unit_test_cmp('3D R voltages', vs.volt(1:12:72)-vs.volt(1), ...
0358                  R(1)*current*linspace(0,1,6)', 1e-10);
0359    unit_test_cmp('3D R z_contact', ...
0360                  [diff(vs.volt([73,1])), diff(vs.volt([74,72]))], ...
0361                  R(2)/2*current*[1,-1], 1e-10);
0362 
0363 
0364 function [R,img] = test_2d_resistor(current,measure)
0365    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0366    z_contact= .1; wid = 3; len = 12; 
0367 
0368    fmdl=mk_grid_model([],linspace(0,wid,3), linspace(0,len,4));
0369    fmdl.electrode(1).nodes = find(fmdl.nodes(:,2) ==   0);
0370    fmdl.electrode(2).nodes = find(fmdl.nodes(:,2) == len);
0371    [fmdl.electrode(:).z_contact] = deal(z_contact);
0372    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0373    img= mk_image(fmdl,conduc);
0374 
0375    Block_R = len / wid / conduc;
0376    Contact_R = z_contact/wid;
0377    R = [Block_R, 2*Contact_R];
0378 
0379 function [R,img] = test_3d_resistor(current,measure);;
0380    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0381    z_contact= .1; wid = 2; len = 5; hig=3; 
0382 
0383    fmdl=mk_grid_model([],0:wid, 0:hig, 0:len);
0384    fmdl.electrode(1).nodes = find(fmdl.nodes(:,3) ==   0);
0385    fmdl.electrode(2).nodes = find(fmdl.nodes(:,3) == len);
0386    [fmdl.electrode(:).z_contact] = deal(z_contact);
0387    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0388    img= mk_image(fmdl,conduc);
0389 
0390    Block_R =  len / wid / hig / conduc;
0391    Contact_R = z_contact/(wid*hig);
0392    R = [Block_R, 2*Contact_R];
0393 
0394 function [R,img] = test_3d_resistor_old(current,measure);
0395    ll=5*1; % length
0396    ww=1*2; % width
0397    hh=1*3; % height
0398    conduc= .13;  % conductivity in Ohm-meters
0399    z_contact= 1e-1;
0400    scale = .46;
0401    nn=0;
0402    for z=0:ll; for x=0:ww; for y=0:hh
0403       nn=nn+1;
0404       mdl.nodes(nn,:) = [x,y,z];
0405    end; end; end
0406    mdl= eidors_obj('fwd_model','3D rectangle');
0407    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0408    mdl.nodes= mdl.nodes*scale;
0409    mdl= rmfield(mdl,'coarse2fine');
0410 
0411    mdl.boundary= find_boundary(mdl.elems);
0412    mdl.gnd_node = 1;
0413    elec_nodes= [1:(ww+1)*(hh+1)];
0414    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0415    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0416    stim.stim_pattern= [-1;1]*current;
0417    stim.meas_pattern= [-1,1]*measure;
0418    mdl.stimulation= stim;
0419    mdl.electrode= elec;
0420    mdl = mdl_normalize(mdl,0);
0421 
0422    mdl.solve = @fwd_solve_1st_order;
0423    mdl.system_mat = @system_mat_1st_order;
0424    img= eidors_obj('image','3D rectangle', ...
0425          'elem_data', ones(size(mdl.elems,1),1) * conduc, ...
0426          'fwd_model', mdl); 
0427 
0428    % analytical solution
0429    Block_R =  ll / ww / hh / scale/ conduc;
0430    Contact_R = z_contact/(ww*hh)/scale^2;
0431    % Contact R reflects z_contact / (width/scale)^2. Here we need to use
0432    %  the scale, since this is not reflected in the size of the
0433    %  FEM as created by the grid. This is different to the test_2d_resistor,
0434    %  where the FEM is created scaled, so that the ww
0435    %  don't need to be scaled by the scale parameter.
0436    R =[ Block_R , 2*Contact_R];

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005