FemMatrix

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

 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]=FemMatrix(Node,Element,z);
0002 
0003 %FemMatrix Computes the blocks of the system matrix for 2D EIT with linear and quadratic basis
0004 % Function [Agrad,Kb,M,S,C]=FemMatrix(Node,Element,z);
0005 % computes the matrices needed in the finite element
0006 % approximation of the 2D EIT forward problem.
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.P. Kaipio
0021 % 25.4.1994. Modified 5.9.1994 by M. Vauhkonen for EIT.
0022 % Modified 13.8.1999 and 23.3.2000 for the EIDORS by M. Vauhkonen,
0023 % University of Kuopio, Department of Applied Physics, PO Box 1627,
0024 % FIN-70211 Kuopio, Finland, email: Marko.Vauhkonen@uku.fi
0025 
0026 Nel=max(size(z));                           %The number of electrodes.
0027 NNode=max(size(Node));                      %The number of nodes
0028 NElement=max(size(Element));                %The number of elements
0029 M=sparse(NNode,Nel);
0030 Kb=sparse(NNode,NNode);
0031 Agrad=sparse(NNode^2,NElement);
0032 s=zeros(Nel,1);
0033 g=reshape([Node.Coordinate],2,NNode)';      %Nodes
0034 
0035 for ii=1:NElement 
0036   A=sparse(NNode,NNode);
0037   ind=(Element(ii).Topology);               % Indices f the ii'th element
0038   gg=g(ind,:);                              % A 3x2 or 6x2 matrix of triangle vertices in (x,y) coord.
0039   if max(size(gg))==3                       % First order basis
0040    grint=grinprodgaus(gg,1);
0041   else
0042    grint=grinprodgausquad(gg,1);            % Second order basis
0043   end
0044 
0045  if any([Element(ii).Face{:,3}]),           % Checks if the triangle ii is a triangle that is
0046                                             % under the electrode.
0047     [In,Jn,InE]=find([Element(ii).Face{:,3}]); 
0048     bind=Element(ii).Face{Jn,1};            % Nodes on the boundary
0049     ab=g(bind(:),:);
0050 
0051     if max(size(bind))==2                   % First order basis
0052      bb1=bound1([ab]);Bb1=zeros(max(size(ind)),1);
0053      bb2=bound2([ab]);Bb2=zeros(max(size(ind)));
0054 
0055      s(InE)=s(InE)+1/z(InE)*2*bb1; % 2*bb1 = length of the electrode.
0056 
0057      eind=[find(bind(1)==ind),find(bind(2)==ind)];   
0058     else                    % Second order basis
0059       bb1=boundquad1([ab]);Bb1=zeros(max(size(ind)),1);
0060       bb2=boundquad2([ab]);Bb2=zeros(max(size(ind)));
0061 
0062       s(InE)=s(InE)+1/z(InE)*electrlen([ab]); 
0063 
0064       eind=[find(bind(1)==ind),find(bind(2)==ind),find(bind(3)==ind)];
0065     end 
0066 
0067     Bb1(eind)=bb1;   
0068     M(ind,InE)=M(ind,InE)-1/z(InE)*Bb1;
0069     
0070     Bb2(eind,eind)=bb2;
0071     A(ind,ind)=grint;
0072     Agrad(:,ii)=A(:);
0073     Kb(ind,ind)=Kb(ind,ind)+1/z(InE)*Bb2;    
0074  else                                       % The triangle isn't under the electrode.
0075     A(ind,ind) = grint; 
0076     Agrad(:,ii)=A(:); 
0077  end
0078 end  
0079 
0080 S=sparse(diag(s));
0081 
0082 
0083 [II1,C]=Current(Nel,NNode,'adj');
0084 C=C(:,1:Nel-1);                             % For the voltage reference
0085 C=sparse(C(:,1:Nel-1));                             
0086 
0087 S=C'*S*C;
0088 M=M*C;
0089 
0090 
0091 
0092 
0093 
0094 
0095