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