jacobian.m

来自「MATLAB二维电阻抗断层成像算法!用于医学成像,里面包括有限元剖分正问题,及反」· M 代码 · 共 63 行

function [J]=Jacobian(Node,Element,Agrad,U,U0,rho,style);

% computes the Jacobian for 2D EIT when elementwise basis is used.
%
% INPUT
% Node = 剖分节点结构数组
% Element =剖分单元结构数组
% Agrad = \int_{Element(ii) \nabla\phi_i\cdot\nabla\phi_j
% U = 电极处的电压
% U0 = 计算场域内的电压
% rho = 电阻率或电导率向量
% style = 重构的电阻率或电导率为实数或复数 
%
% OUTPUT
% J = Jacobian

NNode=max(size(Node));       %剖分的节点数
NElement=max(size(Element)); %剖分的单元数

if nargin<3 
Agrad=sparse(NNode^2,NElement);%初始化单元的整体雅可比矩阵(梯度阵)

 %*******************************************************
 for ii=1:NElement          %对每个单元进行处理
  ind=Element(ii).Topology; %提取该单元的拓朴结构（三个相关节点）
  gg=reshape([Node(ind).Coordinate],2,max(size(ind)))'; 
  % 三个相关节点的坐标阵3*2
  if max(size(ind))==3      %判断是否为一阶单元
   anis=grinprodgaus(gg,1); %调用该单元三节点的梯度计算
  else
   anis=grinprodgausquad(gg,1);% 二阶单元处理
  end
  Aa=sparse(NNode,NNode);   %初始化雅可比矩阵的子矩阵  
  Aa(ind,ind)=anis;         %赋值雅可比矩阵的子矩阵
  Agrad(:,ii)=Aa(:);        %由子矩阵封装构成雅可比整体矩阵
 end
 J=Agrad;
 %*******************************************************
else
if style=='comp'            %复数计算则如下
  J=zeros(size(U,2)*size(U0,2),2*size(Agrad,2));
 for ii=1:size(Agrad,2);
    Agrad1=reshape(Agrad(:,ii),NNode,NNode);
    AgU1=Agrad1*U(1:NNode,:);
    AgU2=Agrad1*U(NNode+1:2*NNode,:);
    JJ=-U0.'*[AgU1;AgU2];
    JJ=JJ(:);
    J(:,ii)=JJ;
    JJ=-U0.'*[-AgU2;AgU1];
    J(:,ii+size(Agrad,2))=JJ(:);
 end
elseif style=='real'       %实数计算则如下
 J=zeros(size(U,2)*size(U0,2),size(Agrad,2));
 for ii=1:size(Agrad,2);  
  JJ=U0.'*1/rho(ii)^2*reshape(Agrad(:,ii),NNode,NNode)*U;
  JJ=JJ(:);
  J(:,ii)=JJ;
 end
end
end