jacobian_elem2nodes

 JACOBIAN_ELEM2NODES: calculate Jacobian on Nodes from Elem solver
 Calculate Jacobian Matrix for EIT Alg of Adler & Guardo 1996
 J         = Jacobian matrix
 fwd_model = forward model defined on nodes (elems may not be defined)
 jacobian_elem2nodes.fwd_model 
           = full forward model (with nodes and elements)

 img = image background for jacobian calc

0001 function J= jacobian_elem2nodes( fwd_model, img);
0002 % JACOBIAN_ELEM2NODES: calculate Jacobian on Nodes from Elem solver
0003 % Calculate Jacobian Matrix for EIT Alg of Adler & Guardo 1996
0004 % J         = Jacobian matrix
0005 % fwd_model = forward model defined on nodes (elems may not be defined)
0006 % jacobian_elem2nodes.fwd_model
0007 %           = full forward model (with nodes and elements)
0008 %
0009 % img = image background for jacobian calc
0010 
0011 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: jacobian_elem2nodes.html 2819 2011-09-07 16:43:11Z aadler $
0013 
0014 fem_fmdl= fwd_model.jacobian_elem2nodes.fwd_model;
0015 EtoN = mapper_nodes_elems( fem_fmdl);
0016 
0017 if ~isfield(img,'elem_data')
0018   img.elem_data= calc_elem_from_node_image(EtoN, img.node_data);
0019 end
0020 
0021 J= calc_jacobian(fem_fmdl, img)*EtoN';
0022 
0023 % Create an image on the elements with a fwd_model on the elemtents
0024 function e_d = calc_elem_from_node_image(EtoN, node_data); 
0025    n_elems= size(EtoN,2);
0026    mu=1e-3; % regularization hyperparameter (squared)
0027    % jnkflag needed to make lsqr shut up
0028    [e_d, jnkflag] = lsqr([EtoN;mu*speye(n_elems)], ...
0029                          [node_data;zeros(n_elems,1)], ...
0030                           1e-8,1000); % should be enough, accuracy is ok
0031 
0032 return;
0033 % Three different ways to invert iEtoN. We only need something fairly
0034 % approximate. This is test code to check
0035 [Nn,Ne]= size(EtoN);  
0036 n1= ones(Nn,1); mu=1e-3;
0037 t=cputime;
0038    ed1= (EtoN'/(EtoN*EtoN'+mu^2*speye(Nn)))*n1;
0039 disp([cputime-t, std(ed1)]);
0040 t=cputime;
0041    ed2= ((EtoN'*EtoN+mu^2*speye(Ne))\EtoN')*n1;
0042 disp([cputime-t, std(ed2)]);
0043 t=cputime;
0044 % Use model saying [E;m*I]*i= [n;0]
0045    ed3= lsqr([EtoN;mu*speye(Ne)],[n1;zeros(Ne,1)],1e-8,1000);
0046 disp([cputime-t, std(ed3)]);