fwd_solve_1st_order

 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

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 6005 2019-06-29 13:26:42Z 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 all dirichlet patterns are the same, then calc in one go
0072    if length(dirichlet_nodes) == 1; rhs = 1:size(pp.QQ,2);
0073    else                           ; rhs = i; end
0074    v(idx,rhs)= left_divide( s_mat.E(idx,idx), ...
0075              neumann_nodes{i}(idx,:) - s_mat.E(idx,:)*dirichlet_values{i});
0076 end
0077 
0078 % If model has a ground node, check if current flowing in this node
0079 if has_gnd_node
0080    Ignd = s_mat.E(dirichlet_nodes{1},:)*v;
0081    Irel = Ignd./sum(abs(pp.QQ)); % relative
0082    if max(abs(Irel))>1e-6
0083       warning('current flowing through ground node. Check stimulation pattern')
0084    end
0085 end
0086 
0087 % calc voltage on electrodes
0088 
0089 % This is horribly inefficient, override
0090 % v_els= pp.N2E * v;
0091 idx = find(any(pp.N2E));
0092 v_els= pp.N2E(:,idx) * v(idx,:);
0093 
0094 
0095 % create a data structure to return
0096 data.meas= meas_from_v_els(v_els, fwd_model.stimulation);
0097 data.time= NaN; % unknown
0098 data.name= 'solved by fwd_solve_1st_order';
0099 try; if img.fwd_solve.get_all_meas == 1
0100    outmap = pp.mr_mapper(1:pp.n_node);
0101    data.volt(outmap,:) = v(1:pp.n_node,:); % but not on CEM nodes
0102 end; end
0103 try; if img.fwd_solve.get_all_nodes== 1
0104    data.volt(pp.mr_mapper,:) = v;                % all, including CEM nodes
0105 end; end
0106 try; if img.fwd_solve.get_elec_curr== 1
0107 %  data.elec_curr = pp.N2E * s_mat.E * v;
0108    idx = find(any(pp.N2E));
0109    data.elec_curr = pp.N2E(:,idx) * s_mat.E(idx,:) * v;
0110 end; end
0111 
0112 
0113 % has_gnd_node = flag if the model has a gnd_node => can warn if current flows
0114 function [dirichlet_nodes, dirichlet_values, neumann_nodes, has_gnd_node]= ...
0115             find_dirichlet_nodes( fwd_model, pp );
0116    fnanQQ = isnan(pp.QQ);
0117    lstims = size(pp.QQ,2);
0118    % Can't use any(...) because if does implicit all
0119    if any(any(fnanQQ))
0120       has_gnd_node = 0; % no ground node is specified
0121       % Are all dirichlet_nodes the same
0122 
0123       % Don't use all on sparse, it will make them full
0124       % Check if all rows are the same
0125       if ~any(any(fnanQQ(:,1)*ones(1,lstims) - fnanQQ,2))
0126          dirichlet_nodes{1} = find(fnanQQ(:,1));
0127          dirichlet_values{1} = sparse(size(pp.N2E,2), size(fnanQQ,2));
0128          dirichlet_values{1}(fnanQQ) = pp.VV(fnanQQ);
0129          neumann_nodes{1} = pp.QQ;
0130          neumann_nodes{1}(fnanQQ) = 0;
0131       else % one at a time
0132          for i=1:size(fnanQQ,2)
0133             fnanQQi= fnanQQ(:,i);
0134             if any(fnanQQi)
0135                dirichlet_nodes{i} = find(fnanQQi);
0136                dirichlet_values{i} = sparse(size(pp.N2E,2), 1);
0137                dirichlet_values{i}(fnanQQi) = pp.VV(fnanQQi,i);
0138                neumann_nodes{i} = pp.QQ(:,i);
0139                neumann_nodes{i}(fnanQQi) = 0;
0140             elseif isfield(pp,'gnd_node')
0141                dirichlet_nodes{i} = pp.gnd_node;
0142                dirichlet_values{i} = sparse(size(pp.N2E,2), 1);
0143                neumann_nodes{1}   = pp.QQ(:,i);
0144                has_gnd_node= 1;
0145             else
0146                error('no required ground node on model');
0147             end
0148          end
0149       end
0150    elseif isfield(pp,'gnd_node')
0151       dirichlet_nodes{1} = pp.gnd_node;
0152       dirichlet_values{1} = sparse(size(pp.N2E,2), size(fnanQQ,2));
0153       neumann_nodes{1}   = pp.QQ;
0154       has_gnd_node= 1;
0155    else
0156       error('no required ground node on model');
0157    end
0158 
0159 function pp = set_gnd_node(fwd_model, pp);
0160    if isfield(fwd_model,'gnd_node');
0161       pp.gnd_node = fwd_model.gnd_node;
0162    else
0163       % try to find one in the model center
0164       ctr =  mean(fwd_model.nodes,1);
0165       d2  =  sum(bsxfun(@minus,fwd_model.nodes,ctr).^2,2);
0166       [~,pp.gnd_node] = min(d2);
0167       eidors_msg('Warning: no ground node found: choosing node %d',pp.gnd_node(1),1);
0168    end
0169 
0170 function vv = meas_from_v_els( v_els, stim)
0171    try
0172 % Was 1.82s
0173 %        % measured voltages from v
0174 %    %   vv = zeros( pp.n_meas, 1 );
0175 %        idx=0;
0176 %        for i=1:length(stim)
0177 %           meas_pat= stim(i).meas_pattern;
0178 %           n_meas  = size(meas_pat,1);
0179 %           vv( idx+(1:n_meas) ) = meas_pat*v_els(:,i);
0180 %           idx= idx+ n_meas;
0181 %        end
0182 
0183 % This code replaced the previous - Nov 4, 2013
0184 % Now 0.437s
0185 % Why is it faster??
0186 
0187        [n_elec,n_stim] = size(v_els);
0188 
0189        copt.cache_obj = {stim};
0190        copt.fstr = 'v2meas';
0191        copt.log_level = 4;
0192        v2meas = eidors_cache(@get_v2meas, {n_elec,n_stim,stim}, copt);
0193        vv = v2meas' * v_els(:);
0194    catch err
0195       if strcmp(err.identifier, 'MATLAB:innerdim');
0196           error(['measurement pattern not compatible with number' ...
0197                   'of electrodes for stimulation patter %d'],i);
0198       else
0199           rethrow(err);
0200       end
0201    end
0202 
0203    
0204 function v2meas = get_v2meas(n_elec,n_stim,stim)
0205     v2meas = sparse(n_elec*n_stim,0);
0206     for i=1:n_stim
0207         meas_pat= stim(i).meas_pattern;
0208         n_meas  = size(meas_pat,1);
0209         v2meas((i-1)*n_elec + 1: i*n_elec,end+(1:n_meas)) = meas_pat';
0210     end
0211 
0212 
0213 function [E, m_idx, pp] = mdl_reduction(E, mr, img, pp);
0214    % if mr is a string we assume it's a function name
0215    if isa(mr,'function_handle') || ischar(mr)
0216       mr = feval(mr,img.fwd_model);
0217    end
0218    % mr is now a struct with fields: main_region, regions
0219    m_idx = mr.main_region;
0220    E = E(m_idx, m_idx);
0221    for i=1:length(mr.regions)
0222       invEi=   mr.regions(i).invE;
0223 % FIXME:!!! data_mapper has done the c2f. But we don't want that here.
0224 %  kludge is to reach into the fine model field. This is only ok because
0225 %  model_reduction is only valid if one parameter describes each field
0226       field = mr.regions(i).field; 
0227       field = find(img.fwd_model.coarse2fine(:,field));
0228       field = field(1); % they're all the same - by def of model_reduction
0229       sigma = img.elem_data(field);
0230       E = E - sigma*invEi;
0231    end
0232 
0233    % Adjust the applied current and measurement matrices
0234    pp.QQ = pp.QQ(m_idx,:);
0235    pp.VV = pp.VV(m_idx,:);
0236    pp.N2E= pp.N2E(:,m_idx);
0237    pp.mr_mapper = cumsum(m_idx); %must be logical
0238    pp.gnd_node = pp.mr_mapper(pp.gnd_node);
0239    if pp.gnd_node==0
0240       error('model_reduction removes ground node');
0241    end
0242    
0243         
0244 function unit_test_voltage_stims;
0245    stim = zeros(16,1); volt=stim; stim([1,4]) = NaN; volt([1,4]) = [1,2];
0246    stimulv.stim_pattern = stim;
0247    stimulv.volt_pattern = volt;
0248    stimulv.meas_pattern = [1,0,0,-1,zeros(1,12)];
0249 
0250    img = mk_image( mk_common_model('a2c2',16),1);
0251    img.fwd_model.stimulation = stimulv;
0252    img.fwd_solve.get_all_nodes = 1;
0253    vh = fwd_solve_1st_order(img);
0254    unit_test_cmp('a2c2 Vstim #1', vh.meas, diff(volt(1:2,1)), 1e-14);
0255    unit_test_cmp('a2c2 Vstim #2', vh.volt(27+[1,4],1), volt([1,4]), 1e-14);
0256    tst = [ ...
0257    1.503131926779798; 1.412534629974291; 1.529078332819747;
0258    1.354399248512161; 1.546241676995996];
0259    unit_test_cmp('a2c2 Vstim #3', vh.volt(1:5:25,1), tst, 1e-14);
0260 
0261    img.fwd_model.stimulation(2) = stimulv;
0262    img.fwd_model.stimulation(1).stim_pattern([1,4]) = 0;
0263    vh = fwd_solve_1st_order(img);
0264    unit_test_cmp('a2c2 Vstim #4', vh.volt(1:5:25,2), tst, 1e-14);
0265 
0266    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt;
0267    imgn.calc_colours.clim = 1; subplot(221); show_fem(imgn,1);
0268 
0269    img = mk_image( mk_common_model('a2C2',16),1);
0270    img.fwd_model.stimulation = stimulv;
0271    img.fwd_solve.get_all_nodes = 1;
0272    vh = fwd_solve_1st_order(img);
0273    unit_test_cmp('a2C2 Vstim #1', vh.meas, diff(volt(1:2)), 1e-14);
0274    unit_test_cmp('a2C2 Vstim #2', vh.volt(num_nodes(img)+[1,4]), volt([1,4]), 1e-14);
0275    tst = [ ...
0276    1.499999999999998; 1.302478674263331; 1.609665333411830; ...
0277    1.215039511028270; 1.691145536046686];
0278    unit_test_cmp('a2C2 Vstim #3', vh.volt(1:5:25), tst, 1e-13);
0279 
0280    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt(1:num_nodes(img));
0281    imgn.calc_colours.clim = 1; subplot(222); show_fem(imgn,1);
0282 
0283    stim = zeros(16,1); volt=stim; stim([1,4]) = NaN; stim(8)=1; volt([1,4]) = [1,1];
0284    img.fwd_model.stimulation(2).stim_pattern = stim;
0285    img.fwd_model.stimulation(2).volt_pattern = volt;
0286    img.fwd_model.stimulation(2).meas_pattern = [1,-1,zeros(1,14)];
0287    vh = fwd_solve_1st_order(img);
0288    unit_test_cmp('a2C2 Vstim #4', vh.volt(num_nodes(img)+[1,4],1), [1;2], 1e-14);
0289    unit_test_cmp('a2C2 Vstim #5', vh.volt(1:5:25,1), tst, 1e-13);
0290    unit_test_cmp('a2C2 Vstim #6', vh.volt(num_nodes(img)+[1,4],2), [1;1], 1e-14);
0291    tst = [ 1.029942389400905; 1.024198991581187; ...
0292            1.048244746016660; 1.006551737030278; 1.057453501332724];
0293    unit_test_cmp('a2C2 Vstim #7', vh.volt(1:5:25,2), tst, 1e-13); % needs weaker tolerance
0294 
0295    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt(1:num_nodes(img),2);
0296    subplot(223); show_fem(imgn,1);
0297 
0298    stim = zeros(16,1); volt=stim; stim([3,6]) = NaN;  volt([3,6]) = [1,2];
0299    img.fwd_model.stimulation(3).stim_pattern = stim;
0300    img.fwd_model.stimulation(3).volt_pattern = volt;
0301    img.fwd_model.stimulation(3).meas_pattern = [1,-1,zeros(1,14)];
0302    vh = fwd_solve_1st_order(img);
0303 
0304    imgn = rmfield(img,'elem_data'); imgn.node_data = vh.volt(1:num_nodes(img),3);
0305    imgn.calc_colours.clim = 1; subplot(224); show_fem(imgn,1);
0306 
0307    unit_test_cmp('a2C2 Vstim #7', vh.volt(num_nodes(img)+[3,6],3), [1;2], 1e-14);
0308 
0309 
0310 
0311 function do_unit_test
0312    unit_test_voltage_stims;
0313 
0314 
0315    img = mk_image( mk_common_model('b2c2',16),1);
0316    vh = fwd_solve_1st_order(img);
0317    tst = [ ...
0318     0.959567140078593; 0.422175237237900; 0.252450963869202; ...
0319     0.180376116490602; 0.143799778367518];
0320    unit_test_cmp('b2c2 TEST', vh.meas(1:5), tst, 1e-12);
0321 
0322    img.fwd_model = rmfield(img.fwd_model,'gnd_node');
0323    vh = fwd_solve_1st_order(img);
0324    unit_test_cmp('b2c2 gnd_node', vh.meas(1:5), tst, 1e-12);
0325 
0326    img.fwd_solve.get_elec_curr = 1;
0327    vh = fwd_solve_1st_order(img);
0328    pp = fwd_model_parameters( img.fwd_model); EC = pp.N2E*pp.QQ;
0329    unit_test_cmp('b2b2 (CEM) elec_curr', vh.elec_curr, EC, 1e-11);
0330 
0331    img.fwd_solve.get_all_meas = 1;
0332    vh = fwd_solve_1st_order(img);
0333     plot(vh.volt);
0334 
0335    img = mk_image( mk_common_model('b2C2',16),1);
0336    vh = fwd_solve_1st_order(img);
0337    tst = [ 0.385629619754662; 0.235061644846908; 0.172837756982388
0338            0.142197580506776; 0.126808900182258; 0.120605655110661];
0339    unit_test_cmp('b2C2 (CEM) TEST', vh.meas(15:20), tst, 1e-12);
0340 
0341    img.fwd_solve.get_elec_curr = 1;
0342    vh = fwd_solve_1st_order(img);
0343    pp = fwd_model_parameters( img.fwd_model); EC = pp.N2E*pp.QQ;
0344    unit_test_cmp('b2C2 (CEM) elec_curr', vh.elec_curr, EC, 1e-11);
0345 
0346    % bad stim patterns (flow through ground node)
0347    img.fwd_model.stimulation(1).stim_pattern(2) = 0;
0348    vh = fwd_solve_1st_order(img);
0349    lastw = lastwarn;
0350    unit_test_cmp('gnd_node warning', lastw, ...
0351     'current flowing through ground node. Check stimulation pattern');
0352 
0353 
0354    %2D resistor
0355    current = 4; measure=1;
0356    [R,img] = test_2d_resistor(current,measure);
0357    img.fwd_solve.get_all_nodes = 1;
0358    vs = fwd_solve_1st_order( img);
0359    va= measure*current*sum(R); % analytic
0360    unit_test_cmp('2D resistor test', va, vs.meas, 1e-12);
0361 
0362    unit_test_cmp('2D R z_contact', ...
0363                  [diff(vs.volt([13,1])), diff(vs.volt([14,12]))], ...
0364                  R(2)/2*current*[1,-1], 1e-12);
0365    unit_test_cmp('2D R voltages', vs.volt(1:3:10)-vs.volt(1), ...
0366                  R(1)*current*linspace(0,1,4)', 1e-12);
0367 
0368    [R,img] = test_2d_resistor_faces(current,measure);
0369    vs = fwd_solve_1st_order( img);
0370    unit_test_cmp('2D resistor faces', va, vs.meas, 1e-12);
0371 
0372    %3D resistor
0373    [R,img] = test_3d_resistor(current,measure);
0374    img.fwd_solve.get_all_nodes = 1;
0375    vs = fwd_solve_1st_order( img);
0376    va= current*sum(R);
0377    unit_test_cmp('3D resistor test', va, vs.meas, 1e-10);
0378    unit_test_cmp('3D R voltages', vs.volt(1:12:72)-vs.volt(1), ...
0379                  R(1)*current*linspace(0,1,6)', 1e-10);
0380    unit_test_cmp('3D R z_contact', ...
0381                  [diff(vs.volt([73,1])), diff(vs.volt([74,72]))], ...
0382                  R(2)/2*current*[1,-1], 1e-10);
0383 
0384    [R,img] = test_3d_resistor_faces(current,measure);
0385    vs = fwd_solve_1st_order( img);
0386    unit_test_cmp('3D resistor faces', va, vs.meas, 1e-10);
0387 
0388 
0389 function [R,img] = test_2d_resistor(current,measure)
0390    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0391    z_contact= .1; wid = 3; len = 12; 
0392 
0393    fmdl=mk_grid_model([],linspace(0,wid,3), linspace(0,len,4));
0394    fmdl.electrode(1).nodes = find(fmdl.nodes(:,2) ==   0);
0395    fmdl.electrode(2).nodes = find(fmdl.nodes(:,2) == len);
0396    [fmdl.electrode(:).z_contact] = deal(z_contact);
0397    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0398    img= mk_image(fmdl,conduc);
0399 
0400    Block_R = len / wid / conduc;
0401    Contact_R = z_contact/wid;
0402    R = [Block_R, 2*Contact_R];
0403 
0404 % define electrode using face rather than nodes
0405 function [R,img] = test_2d_resistor_faces(current,measure)
0406    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0407    z_contact= .1; wid = 3; len = 12; 
0408 
0409    fmdl=mk_grid_model([],linspace(0,wid,3), linspace(0,len,4));
0410    bdy = fmdl.boundary;
0411    bdy( any(reshape(fmdl.nodes(bdy,2),size(bdy))>0,2),:)=[];
0412    fmdl.electrode(1).nodes = [];
0413    fmdl.electrode(1).faces = bdy;
0414    fmdl.electrode(2).nodes = find(fmdl.nodes(:,2) == len);
0415    [fmdl.electrode(:).z_contact] = deal(z_contact);
0416    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0417    img= mk_image(fmdl,conduc);
0418 
0419    Block_R = len / wid / conduc;
0420    Contact_R = z_contact/wid;
0421    R = [Block_R, 2*Contact_R];
0422 
0423 function [R,img] = test_3d_resistor(current,measure);;
0424    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0425    z_contact= .1; wid = 2; len = 5; hig=3; 
0426 
0427    fmdl=mk_grid_model([],0:wid, 0:hig, 0:len);
0428    fmdl.electrode(1).nodes = find(fmdl.nodes(:,3) ==   0);
0429    fmdl.electrode(2).nodes = find(fmdl.nodes(:,3) == len);
0430    [fmdl.electrode(:).z_contact] = deal(z_contact);
0431    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0432    img= mk_image(fmdl,conduc);
0433 
0434    Block_R =  len / wid / hig / conduc;
0435    Contact_R = z_contact/(wid*hig);
0436    R = [Block_R, 2*Contact_R];
0437 
0438 % define electrode using face rather than nodes
0439 function [R,img] = test_3d_resistor_faces(current,measure);;
0440    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0441    z_contact= .1; wid = 2; len = 5; hig=3; 
0442 
0443    fmdl=mk_grid_model([],0:wid, 0:hig, 0:len);
0444 %  fmdl.electrode(1).nodes = find(fmdl.nodes(:,3) ==   0);
0445    bdy = fmdl.boundary;
0446    bdy( any(reshape(fmdl.nodes(bdy,3),size(bdy))>0,2),:)=[];
0447    fmdl.electrode(1).nodes = [];
0448    fmdl.electrode(1).faces = bdy;
0449    fmdl.electrode(2).nodes = find(fmdl.nodes(:,3) == len);
0450    [fmdl.electrode(:).z_contact] = deal(z_contact);
0451    fmdl.stimulation = stim_meas_list([1,2,1,2],2,current,measure);
0452    img= mk_image(fmdl,conduc);
0453 
0454    Block_R =  len / wid / hig / conduc;
0455    Contact_R = z_contact/(wid*hig);
0456    R = [Block_R, 2*Contact_R];
0457 
0458 function [R,img] = test_3d_resistor_old(current,measure);
0459    ll=5*1; % length
0460    ww=1*2; % width
0461    hh=1*3; % height
0462    conduc= .13;  % conductivity in Ohm-meters
0463    z_contact= 1e-1;
0464    scale = .46;
0465    nn=0;
0466    for z=0:ll; for x=0:ww; for y=0:hh
0467       nn=nn+1;
0468       mdl.nodes(nn,:) = [x,y,z];
0469    end; end; end
0470    mdl= eidors_obj('fwd_model','3D rectangle');
0471    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0472    mdl.nodes= mdl.nodes*scale;
0473    mdl= rmfield(mdl,'coarse2fine');
0474 
0475    mdl.boundary= find_boundary(mdl.elems);
0476    mdl.gnd_node = 1;
0477    elec_nodes= [1:(ww+1)*(hh+1)];
0478    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0479    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0480    stim.stim_pattern= [-1;1]*current;
0481    stim.meas_pattern= [-1,1]*measure;
0482    mdl.stimulation= stim;
0483    mdl.electrode= elec;
0484    mdl = mdl_normalize(mdl,0);
0485 
0486    mdl.solve = @fwd_solve_1st_order;
0487    mdl.system_mat = @system_mat_1st_order;
0488    img= eidors_obj('image','3D rectangle', ...
0489          'elem_data', ones(size(mdl.elems,1),1) * conduc, ...
0490          'fwd_model', mdl); 
0491 
0492    % analytical solution
0493    Block_R =  ll / ww / hh / scale/ conduc;
0494    Contact_R = z_contact/(ww*hh)/scale^2;
0495    % Contact R reflects z_contact / (width/scale)^2. Here we need to use
0496    %  the scale, since this is not reflected in the size of the
0497    %  FEM as created by the grid. This is different to the test_2d_resistor,
0498    %  where the FEM is created scaled, so that the ww
0499    %  don't need to be scaled by the scale parameter.
0500    R =[ Block_R , 2*Contact_R];