JacobianNode

JacobianNode Computes the Jacobian for 2D EIT when node basis is used
 Function [J]=JacobianNode(g,H,Agrad,U,U0,rho,style);
 computes the Jacobian for 2D EIT when node 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]=JacobianNode(Node,Element,Agrad,U,U0,rho,style);
0002 
0003 %JacobianNode Computes the Jacobian for 2D EIT when node basis is used
0004 % Function [J]=JacobianNode(g,H,Agrad,U,U0,rho,style);
0005 % computes the Jacobian for 2D EIT when node 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 
0032 
0033 if nargin<3
0034 
0035 mE=max(size(Element(1).Topology));
0036 if mE==3
0037  Agrad=sparse(NNode^2,NNode);
0038  H=reshape([Element.Topology],3,NElement)';
0039  mH=max(max(H)); 
0040 else
0041  H=reshape([Element.Topology],6,NElement)';
0042  mH=max(max(H(:,1:2:6))); 
0043  Agrad=sparse(NNode^2,mH);
0044  clear H
0045 end
0046 g=reshape([Node.Coordinate],2,NNode)';      %Nodes
0047 
0048 for jj=1:mH
0049 Aa=sparse(NNode,NNode);
0050 El=Node(jj).ElementConnection;
0051  for ii=1:max(size(El))
0052    ind=Element(El(ii)).Topology; % Indices of the element
0053    gg=g(ind,:);
0054    if max(size(gg))==3
0055     indsig=ind;
0056     I=find(jj==indsig);
0057     anis=grinprodgausnode(gg,I);
0058     Aa(ind,ind)=Aa(ind,ind)+anis;
0059    else
0060      indsig=ind(1:2:6);
0061      I=find(jj==indsig);
0062      anis=grinprodgaus2ndnode(gg,I);
0063      Aa(ind,ind)=Aa(ind,ind)+anis;         
0064    end  
0065  end  
0066  Agrad(:,jj)=Aa(:);
0067 end
0068 J=Agrad;
0069 else
0070 if style=='comp'
0071   J=zeros(size(U,2)*size(U0,2),2*size(Agrad,2));
0072  for ii=1:size(Agrad,2);
0073     Agrad1=reshape(Agrad(:,ii),NNode,NNode);
0074     AgU1=Agrad1*U(1:NNode,:);
0075     AgU2=Agrad1*U(NNode+1:2*NNode,:);
0076     JJ=-U0.'*[AgU1;AgU2];
0077     JJ=JJ(:);
0078     J(:,ii)=JJ;
0079     JJ=-U0.'*[-AgU2;AgU1];
0080     J(:,ii+size(Agrad,2))=JJ(:);
0081  end
0082 elseif style=='real'
0083   J=zeros(size(U,2)*size(U0,2),size(Agrad,2));
0084    for ii=1:size(Agrad,2);
0085     JJ=-U0.'*reshape(Agrad(:,ii),NNode,NNode)*U;
0086     JJ=JJ(:);
0087     J(:,ii)=JJ;
0088    end
0089  end
0090 end
0091 
0092 
0093 
0094 
0095 
0096 
0097 
0098