jacobian_elem2nodes

 JACOBIAN_ELEM2NODES: calculate Jacobian on Nodes from Elem solver
 Calculate Jacobian Matrix for EIT Alg of Adler & Guardo 1996
 J         = Jacobian matrix
 img.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 % img.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.m 6422 2022-11-29 21:24:49Z aadler $
0013 
0014 if nargin == 1
0015    img= fwd_model;
0016 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0017    warning('EIDORS:DeprecatedInterface', ...
0018       ['Calling JACOBIAN_elem2nodes with two arguments is deprecated and will cause' ...
0019        ' an error in a future version. First argument ignored.']);
0020 end
0021 fwd_model= img.fwd_model;
0022 
0023 fem_fmdl= fwd_model.jacobian_elem2nodes.fwd_model;
0024 EtoN = mapper_nodes_elems( fem_fmdl);
0025 
0026 if ~isfield(img,'elem_data')
0027   img.elem_data= calc_elem_from_node_image(EtoN, img.node_data);
0028 end
0029 img.fwd_model = fem_fmdl;
0030 J= calc_jacobian(img)*EtoN';
0031 
0032 % Create an image on the elements with a fwd_model on the elemtents
0033 function e_d = calc_elem_from_node_image(EtoN, node_data); 
0034    n_elems= size(EtoN,2);
0035    mu=1e-3; % regularization hyperparameter (squared)
0036    % jnkflag needed to make lsqr shut up
0037 
0038    A = [EtoN;mu*speye(n_elems)];
0039    B = [node_data;zeros(n_elems,1)];
0040    if 0 % previously used iterative
0041       [e_d, jnkflag] = lsqr(A,B, 1e-8,1000);
0042    else
0043       e_d = left_divide(A,B);
0044    end
0045 
0046 return;
0047 % Test code in nodal_jacobian01
0048 
0049 
0050 % Three different ways to invert iEtoN. We only need something fairly
0051 % approximate. This is test code to check
0052 [Nn,Ne]= size(EtoN);  
0053 n1= ones(Nn,1); mu=1e-3;
0054 t=cputime;
0055    ed1= (EtoN'/(EtoN*EtoN'+mu^2*speye(Nn)))*n1;
0056 disp([cputime-t, std(ed1)]);
0057 t=cputime;
0058    ed2= ((EtoN'*EtoN+mu^2*speye(Ne))\EtoN')*n1;
0059 disp([cputime-t, std(ed2)]);
0060 t=cputime;
0061 % Use model saying [E;m*I]*i= [n;0]
0062    ed3= lsqr([EtoN;mu*speye(Ne)],[n1;zeros(Ne,1)],1e-8,1000);
0063 disp([cputime-t, std(ed3)]);