hsfemmatrix.m

MATLAB二维电阻抗断层成像算法!用于医学成像,里面包括有限元剖分正问题,及反问题的算法.并且附有网络剖分数据表!

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 matrix

Nel=max(size(z));                           %电极数
NNode=max(size(Node));                      %节点数
NElement=max(size(Element));                %单元数
M=sparse(NNode,Nel);            %以稀疏阵初始化M
Kb=sparse(NNode,NNode);         %以稀疏阵初始化Kb
Agrad=sparse(NNode^2,NElement); %以稀疏阵初始化Agrad
s=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(:); 
 end
end  

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;