Jacobian

Jacobian Computes the Jacobian for 2D EIT when elementwise basis is used
 Function [J]=Jacobian(g,H,Agrad,U,U0,rho,style);
 computes the Jacobian for 2D EIT when elementwise basis
 is used.

 INPUT

 Node = nodal data structure
 Element = element data structure
 Agrad = \int_{Element(ii) \nabla\phi_i\cdot\nabla\phi_j
 U = voltages corresponding to the injected currents
 U0 = voltages corresponding to the measurement field
 rho = resistivity or admittivity vector
 style = either 'real' for reconstructing resistivity or 'comp' 
 for reconstructin admittivity 

 OUTPUT

 J = Jacobian

0001 function [J]=Jacobian(Node,Element,Agrad,U,U0,rho,style);
0002 
0003 %Jacobian Computes the Jacobian for 2D EIT when elementwise basis is used
0004 % Function [J]=Jacobian(g,H,Agrad,U,U0,rho,style);
0005 % computes the Jacobian for 2D EIT when elementwise basis
0006 % is used.
0007 %
0008 % INPUT
0009 %
0010 % Node = nodal data structure
0011 % Element = element data structure
0012 % Agrad = \int_{Element(ii) \nabla\phi_i\cdot\nabla\phi_j
0013 % U = voltages corresponding to the injected currents
0014 % U0 = voltages corresponding to the measurement field
0015 % rho = resistivity or admittivity vector
0016 % style = either 'real' for reconstructing resistivity or 'comp'
0017 % for reconstructin admittivity
0018 %
0019 % OUTPUT
0020 %
0021 % J = Jacobian
0022 
0023 % M. Vauhkonen 13.8.1999,
0024 % University of Kuopio, Department of Applied Physics, PO Box 1627,
0025 % FIN-70211 Kuopio, Finland, email: Marko.Vauhkonen@uku.fi
0026 
0027 
0028 NNode=max(size(Node));
0029 NElement=max(size(Element));
0030 
0031 if nargin<3 
0032 Agrad=sparse(NNode^2,NElement); % Gradients of the basis functions integrated over
0033                        % each element.
0034 
0035  for ii=1:NElement
0036   ind=Element(ii).Topology;
0037   gg=reshape([Node(ind).Coordinate],2,max(size(ind)))'; % A 3x2 or 6x2 matrix of triangle vertices in (x,y) coord.
0038   if max(size(ind))==3
0039    anis=grinprodgaus(gg,1);
0040   else
0041    anis=grinprodgausquad(gg,1);
0042   end
0043   Aa=sparse(NNode,NNode);
0044   Aa(ind,ind)=anis;
0045   Agrad(:,ii)=Aa(:);    
0046  end
0047  J=Agrad;
0048 else
0049 if style=='comp'
0050   J=zeros(size(U,2)*size(U0,2),2*size(Agrad,2));
0051  for ii=1:size(Agrad,2);
0052     Agrad1=reshape(Agrad(:,ii),NNode,NNode);
0053     AgU1=Agrad1*U(1:NNode,:);
0054     AgU2=Agrad1*U(NNode+1:2*NNode,:);
0055     JJ=-U0.'*[AgU1;AgU2];
0056     JJ=JJ(:);
0057     J(:,ii)=JJ;
0058     JJ=-U0.'*[-AgU2;AgU1];
0059     J(:,ii+size(Agrad,2))=JJ(:);
0060  end
0061 elseif style=='real' 
0062  J=zeros(size(U,2)*size(U0,2),size(Agrad,2));
0063  for ii=1:size(Agrad,2);
0064   JJ=U0.'*1/rho(ii)^2*reshape(Agrad(:,ii),NNode,NNode)*U;
0065   JJ=JJ(:);
0066   J(:,ii)=JJ;
0067  end
0068 end
0069 end
0070 
0071