function [Agrad,Kb,M,S,C]=FemMatrix(Node,Element,z);%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 matrixNel=max(size(z));                           %電極數NNode=max(size(Node));                      %節點數NElement=max(size(Element));                %單元數M=sparse(NNode,Nel);            %以稀疏陣初始化MKb=sparse(NNode,NNode);         %以稀疏陣初始化KbAgrad=sparse(NNode^2,NElement); %以稀疏陣初始化Agrads=zeros(Nel,1);                 %電極的電導率列向量g=reshape([Node.Coordinate],2,NNode)';      %節點坐標矩陣for ii=1:NElement               %對每個單元進行  A=sparse(NNode,NNode);        %Agrad的預處理矩陣  ind=(Element(ii).Topology);   %第ii單元的拓樸結構 節點號  gg=g(ind,:);                  %第ii單元的拓樸節點坐標矩陣 3x2 or 6x2  if max(size(gg))==3           % 一階單元   grint=grinprodgaus(gg,1);    % 一階單元梯度計算  else   grint=grinprodgausquad(gg,1);% 二階單元梯度計算  end if any([Element(ii).Face{:,3}]),      %如果ii單元在電極下(Element.Face中的信息)則：    [In,Jn,InE]=find([Element(ii).Face{:,3}]);     bind=Element(ii).Face{Jn,1};% 單元在邊界上的邊的節點    ab=g(bind(:),:);            % 單元在邊界上的邊的節點坐標    if max(size(bind))==2       % 一階單元邊？     bb1=bound1([ab]);Bb1=zeros(max(size(ind)),1);%計算bound1     bb2=bound2([ab]);Bb2=zeros(max(size(ind)));  %計算bound2     s(InE)=s(InE)+1/z(InE)*2*bb1; % 2*bb1 = length of the electrode.             eind=[find(bind(1)==ind),find(bind(2)==ind)];       else				    % Second order basis      bb1=boundquad1([ab]);Bb1=zeros(max(size(ind)),1);      bb2=boundquad2([ab]);Bb2=zeros(max(size(ind)));      s(InE)=s(InE)+1/z(InE)*electrlen([ab]);       eind=[find(bind(1)==ind),find(bind(2)==ind),find(bind(3)==ind)];    end     Bb1(eind)=bb1;       M(ind,InE)=M(ind,InE)-1/z(InE)*Bb1;%計算M矩陣        Bb2(eind,eind)=bb2;    A(ind,ind)=grint;    Agrad(:,ii)=A(:);                  %封裝構成Agrad矩陣    Kb(ind,ind)=Kb(ind,ind)+1/z(InE)*Bb2; %計算Kb矩陣 else                                     %如果單元不在邊界上則：計算。。.    A(ind,ind) = grint;     Agrad(:,ii)=A(:);  endend  S=sparse(diag(s));		 %計算S矩陣           [II1,C]=Current(Nel,NNode,'adj'); %為計算電壓參考C矩陣C=C(:,1:Nel-1);                   C=sparse(C(:,1:Nel-1));                             S=C'*S*C;M=M*C;