np_fwd_solve

 NP_FWD_SOLVE: data= np_fwd_solve( fwd_model, img)
 Fwd solver for Nick Polydorides EIDORS3D code
 Input:
    fwd_model = forward model
    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)

0001 function data= np_fwd_solve( fwd_model, img)
0002 % NP_FWD_SOLVE: data= np_fwd_solve( fwd_model, img)
0003 % Fwd solver for Nick Polydorides EIDORS3D code
0004 % Input:
0005 %    fwd_model = forward model
0006 %    img       = image struct
0007 % Output:
0008 %    data = measurements struct
0009 % Options: (to return internal FEM information)
0010 %    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
0011 %    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
0012 
0013 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0014 % $Id: np_fwd_solve.m 5394 2017-04-12 15:10:30Z aadler $
0015 
0016 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0017 
0018 warning('EIDORS:deprecated','NP_FWD_SOLVE is deprecated as of 07-Jun-2012. Use FWD_SOLVE_1ST_ORDER instead.');
0019 
0020 if nargin==1 %normally takes one parameter
0021    img = fwd_model;
0022    fwd_model = img.fwd_model;
0023 end
0024 
0025 p= np_fwd_parameters( fwd_model );
0026 
0027 %Set the tolerance for the forward solver
0028 tol = 1e-5;
0029 
0030 s_mat= calc_system_mat( fwd_model, img );
0031 
0032 Vfwd = forward_solver(s_mat.E, p.I, tol, s_mat.perm);
0033 
0034 Velec=Vfwd( p.n_node+(1:p.n_elec),:);
0035 voltH = zeros( p.n_meas, 1 );
0036 idx=0;
0037 for i=1:p.n_stim
0038    meas_pat= fwd_model.stimulation(i).meas_pattern;
0039    n_meas  = size(meas_pat,1);
0040    voltH( idx+(1:n_meas) ) = meas_pat*Velec(:,i);
0041    idx= idx+ n_meas;
0042 end
0043 
0044 % create a data structure to return
0045 data.meas= voltH;
0046 data.time= NaN; % unknown
0047 data.name= 'solved by np_fwd_solve';
0048 try; if img.fwd_solve.get_all_meas == 1
0049    data.volt = Vfwd(1:p.n_node,:); % but not on CEM nodes
0050 end; end
0051 try; if img.fwd_solve.get_all_nodes== 1
0052    data.volt = Vfwd;               % all, including CEM nodes
0053 end; end
0054 
0055 function do_unit_test
0056    test_3d_resistor
0057 
0058 function test_3d_resistor
0059    ll=5; % length
0060    ww=1; % width
0061    hh=1; % height
0062    conduc= .13;  % conductivity in Ohm-meters
0063    current= 4;  % Amps
0064    z_contact= 9e-1;
0065    scale = .46;
0066    nn=0;
0067    for z=0:ll; for x=0:ww; for y=0:hh
0068       nn=nn+1;
0069       mdl.nodes(nn,:) = [x,y,z];
0070    end; end; end
0071 
0072    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0073    mdl.nodes= mdl.nodes*scale;
0074    mdl= rmfield(mdl,'coarse2fine');
0075    mdl.boundary= find_boundary(mdl.elems);
0076    mdl.gnd_node = 1;
0077    mdl.normalize_measurements = 0;
0078    elec_nodes= [1:(ww+1)*(hh+1)];
0079    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0080    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0081    stim.stim_pattern= [-1;1]*current;
0082    stim.meas_pattern= [-1,1];
0083    mdl.stimulation= stim;
0084    mdl.electrode= elec;
0085 %  show_fem(mdl);
0086 
0087    % analytical solution
0088    R = ll / ww / hh / scale/ conduc + 2*z_contact/scale^2;
0089 
0090    V= current*R;
0091 %  fprintf('Solver %s: %f\n', 'analytic', V);
0092 
0093    mdl.solve = @np_fwd_solve;
0094    mdl.system_mat = @np_calc_system_mat;
0095 %  mdl.misc.perm_sym= '{n}'; %%% Check it works without
0096    img= mk_image(mdl, conduc);
0097 
0098    fsol= fwd_solve(img);
0099 %  fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0100    unit_test_cmp('np_fwd_solve vs analytic', fsol.meas, V, 1e-11);
0101 
0102 
0103 
0104    mdl.solve = @fwd_solve_1st_order;
0105    mdl.system_mat = @system_mat_1st_order;
0106 
0107    img= mk_image(mdl, conduc);
0108 
0109    fsol= fwd_solve(img);
0110 %  fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0111    unit_test_cmp('new solver 2d vs analytic', fsol.meas, V, 1e-11);