jacobian_elem2nodes

PURPOSE ^

JACOBIAN_ELEM2NODES: calculate Jacobian on Nodes from Elem solver

SYNOPSIS ^

function J= jacobian_elem2nodes( fwd_model, img)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 3819 2013-04-13 07:36:27Z bgrychtol $
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    [e_d, jnkflag] = lsqr([EtoN;mu*speye(n_elems)], ...
0038                          [node_data;zeros(n_elems,1)], ...
0039                           1e-8,1000); % should be enough, accuracy is ok
0040 
0041 return;
0042 % Three different ways to invert iEtoN. We only need something fairly
0043 % approximate. This is test code to check
0044 [Nn,Ne]= size(EtoN);  
0045 n1= ones(Nn,1); mu=1e-3;
0046 t=cputime;
0047    ed1= (EtoN'/(EtoN*EtoN'+mu^2*speye(Nn)))*n1;
0048 disp([cputime-t, std(ed1)]);
0049 t=cputime;
0050    ed2= ((EtoN'*EtoN+mu^2*speye(Ne))\EtoN')*n1;
0051 disp([cputime-t, std(ed2)]);
0052 t=cputime;
0053 % Use model saying [E;m*I]*i= [n;0]
0054    ed3= lsqr([EtoN;mu*speye(Ne)],[n1;zeros(Ne,1)],1e-8,1000);
0055 disp([cputime-t, std(ed3)]);

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005