FemMatrixNode

FemMatrixNode Computes the blocks of the system matrix for 2D EIT with linear and quadratic basis. Conductivity in linear basis.
 Function [Agrad,Kb,M,S,C]=FemMatrixNode(Node,Element,z);
 computes the matrices needed in the finite element
 approximation of the 2D EIT forward problem. Conductivity in linear basis.

 INPUT

 Node = nodal data structure
 Element = element data structure
 z = a vector of (complex) contact impedances

 OUTPUT

 Agrad = the gradient part of the system matrix
 Kb,M and S = other blocks of the system matrix
 C = voltage reference matrix

0001 function [Agrad,Kb,M,S,C]=FemMatrixNode(Node,Element,z);
0002 
0003 %FemMatrixNode Computes the blocks of the system matrix for 2D EIT with linear and quadratic basis. Conductivity in linear basis.
0004 % Function [Agrad,Kb,M,S,C]=FemMatrixNode(Node,Element,z);
0005 % computes the matrices needed in the finite element
0006 % approximation of the 2D EIT forward problem. Conductivity in linear basis.
0007 %
0008 % INPUT
0009 %
0010 % Node = nodal data structure
0011 % Element = element data structure
0012 % z = a vector of (complex) contact impedances
0013 %
0014 % OUTPUT
0015 %
0016 % Agrad = the gradient part of the system matrix
0017 % Kb,M and S = other blocks of the system matrix
0018 % C = voltage reference matrix
0019 
0020 % M. Vauhkonen 11.5.1994, modified from the version of J. Kaipio
0021 % 25.4.1994. Modified 5.9.1994 by M. Vauhkonen for EIT.
0022 % Modified for EIDORS by M. Vauhkonen 11.5.2000
0023 % University of Kuopio, Department of Applied Physics, PO Box 1627,
0024 % FIN-70211 Kuopio, Finland, email: Marko.Vauhkonen@uku.fi
0025 
0026 
0027 Nel=max(size(z));                           %The number of electrodes.
0028 NNode=max(size(Node));                      %The number of nodes
0029 NElement=max(size(Element));                %The number of elements
0030 M=sparse(NNode,Nel);
0031 Kb=sparse(NNode,NNode);
0032 s=zeros(Nel,1);
0033 mE=max(size(Element(1).Topology));
0034 if mE==3
0035  Agrad=sparse(NNode^2,NNode);
0036  H=reshape([Element.Topology],3,NElement)';
0037  mH=max(max(H));
0038 else
0039  H=reshape([Element.Topology],6,NElement)';
0040  mH=max(max(H(:,1:2:6)));
0041  Agrad=sparse(NNode^2,mH);
0042  clear H
0043 end
0044 g=reshape([Node.Coordinate],2,NNode)';      %Nodes
0045 
0046 for jj=1:mH
0047 Aa=sparse(NNode,NNode);
0048 El=Node(jj).ElementConnection;
0049  for ii=1:max(size(El))
0050    ind=Element(El(ii)).Topology; % Indices of the element
0051    gg=g(ind,:);
0052    if max(size(gg))==3
0053     indsig=ind;
0054     I=find(jj==indsig);
0055     anis=grinprodgausnode(gg,I);
0056     Aa(ind,ind)=Aa(ind,ind)+anis;
0057    else
0058     indsig=ind(1:2:6);
0059     I=find(jj==indsig);
0060     anis=grinprodgausnodequad(gg,I);
0061     Aa(ind,ind)=Aa(ind,ind)+anis;
0062    end
0063  end
0064  Agrad(:,jj)=Aa(:);
0065 end
0066 
0067 for ii=1:NElement
0068   ind=(Element(ii).Topology);               % The indices to g of the ii'th triangle.
0069   gg=g(ind,:);                              % A 3x2 or 6x2 matrix of triangle nodes in (x,y) coord.
0070   
0071  if any([Element(ii).Face{:,3}]),           %Checks if the triangle ii is the triangle that is
0072                                             % under the electrode.
0073     [In,Jn,InE]=find([Element(ii).Face{:,3}]);
0074     bind=Element(ii).Face{Jn,1};            % Nodes on the boundary
0075     ab=g(bind(:),:);
0076 
0077     if max(size(bind))==2                   % First order basis
0078      bb1=bound1([ab]);Bb1=zeros(max(size(ind)),1);
0079      bb2=bound2([ab]);Bb2=zeros(max(size(ind)));
0080 
0081      s(InE)=s(InE)+1/z(InE)*2*bb1; % 2*bb1 = length of the electrode.
0082 
0083      eind=[find(bind(1)==ind),find(bind(2)==ind)];
0084     else                                    % Second order basis
0085       bb1=boundquad1([ab]);Bb1=zeros(max(size(ind)),1);
0086       bb2=boundquad2([ab]);Bb2=zeros(max(size(ind)));
0087 
0088       s(InE)=s(InE)+1/z(InE)*electrlen([ab]);
0089 
0090       eind=[find(bind(1)==ind),find(bind(2)==ind),find(bind(3)==ind)];
0091     end
0092 
0093     Bb1(eind)=bb1;
0094     M(ind,InE)=M(ind,InE)-1/z(InE)*Bb1;
0095 
0096     Bb2(eind,eind)=bb2;
0097     Kb(ind,ind)=Kb(ind,ind)+1/z(InE)*Bb2;
0098   else                                      %The triangle isn't under the electrode.
0099   end
0100 end  
0101 
0102 S=sparse(diag(s));
0103 
0104 [II1,C]=Current(Nel,NNode,'adj');
0105 C=C(:,1:Nel-1);                             % For the voltage reference
0106 C=sparse(C(:,1:Nel-1));                             
0107 
0108 S=C'*S*C;
0109 M=M*C;
0110 
0111 
0112 
0113 
0114 
0115 
0116