基于六面体网格的FEM算法——刚体矩阵存在病态的问题

上完学校的课，最近继续开始了脑电正问题建模的问题。在进行学习的过程中发现现有的关于FEM头模型搭建的代码都无法查看，于是开始着手自己编写相关代码。相关代码的编写逻辑见之前的博客——二维FEM建模及求解（入门级：数理方程到代码编写）。

但是在编写六面体网格的相关代码时出现了问题：构造出来的刚体矩阵存在极大的病态，计算了矩阵的条件数，为：1.2502e-25，结果也存在极大的问题。询问了老师，分析的结果是形函数的构造存在问题。在询问的过程中发现了自己最近研究中存在的问题——自己在一个人闭门造车，于是我把自己代码粘在这里，如果有大牛看见，希望不吝赐教。而如果后续解决了该问题，我也会在这里进行补充。

其中形函数的构造是按照该链接进行构造的。为：

代码的逻辑同 二维FEM建模及求解（入门级：数理方程到代码编写），另外需要补充的是，在求积分的过程中，将离散化区域的积分转化为了：均值 X 对应离散区域的面积。下面为对应函数代码：

function [K,u,rcondj] = FEM_hexahedral(cfg, mesh)
% 三维的六面体
% 边界假定为电势为0
% cfg：目前仅存在一个值cfg.conductivity(i)
% mesh包括：mesh.pos——顶点坐标；mesh.hex——六面体单元；mesh.tissue——组织标签
% mesh来自于fieldtrip代码，头模型为实例模型

%% 预处理
nds = mesh.pos;   % 顶点坐标
els = mesh.hex;   % 六面体单元
% 边界获取
% bcs = mesh.bnd((mesh.tissue==3),1:4);
bcs = 0;  % 目前尚未获取到边界元素，故暂时设置为0

% 电导率设置
cond = zeros(length(els),1);
for i = 1:length(cfg.conductivity)
    cond(mesh.tissue==i) = cfg.conductivity(i);
end

I = 1;  % 设置源

%% 计算
nnd = size(nds,1);u = zeros(nnd,1);K =zeros (nnd,nnd);f = zeros(nnd,1);
rcondj = zeros(size(els,1),1);
for j = 1:size(els,1)
    hexpos = nds(els(j,:),:);
    ndim = size(nds,2);  % 维度
    npos = size(els,2);  % 顶点数
    V = hex_volume(hexpos, ndim);  %计算六面体面积
    [kj,fj, rcondJ] = stima(hexpos,V,npos,cond(j),I);
    K(els(j,:),els(j,:)) = K(els(j,:),els(j,:)) + kj ;  % 仅对角线上有值,stima为第j个K
    f(els(j,:)) = f(els(j,:)) + fj;
    rcondj(j) = rcondJ;
end
freends = setdiff(1:nnd,bcs);     % 去除边界上的值，即边界条件u=0的描述。
u(freends) = K(freends,freends)\f(freends);

function volume = hex_volume(hexpos,ndim)
    % hexpos：六面体的8个顶点坐标，维度：8*3
    
    % 按照fieldtrip的顶点排序规则定义六面的各面
    facet = [1,2,3,4;5,6,7,8;1,2,5,6;
    2,3,6,7;3,4,7,8;4,1,8,5];
    
    % 计算中点坐标
    center = mean(hexpos);
    % 计算体积
    volume = 0;
    for k = 1:6
        vertices = hexpos(facet(k,:),:);
        v1 = det([ones(1,4);[vertices(1:3,:);center]'])/(prod(1:ndim));
        v2 = det([ones(1,4);[vertices(2:4,:);center]'])/(prod(1:ndim));
        volume = volume+abs(v1)+abs(v2);
        
    end
end

function [kj,fj, rcondJ] = stima(vertices,V,npos,cond,I)   
    % 计算第j个K
    % vertices：当前六面体的顶点
    % V：当前六面体的体积
    % npos：六面体顶点数
    % cond：当前六面体电导率
    % f：偶极子电流
    
    J = pinv([ones(npos,1),vertices,vertices(:,1).*vertices(:,2),vertices(:,1).*vertices(:,3),...
        vertices(:,2).*vertices(:,3),vertices(:,1).*vertices(:,2).*vertices(:,3)]);
    rcondJ = rcond(J);
    a1 = mean(sum(vertices(:,1:2),2));
    a2 = mean(sum(vertices(:,2:3),2));
    a3 = mean(sum([vertices(:,1),vertices(:,3)],2));
    a4 = mean(sum([vertices(:,1).*vertices(:,2),vertices(:,1).*vertices(:,3),vertices(:,2).*vertices(:,3)],2));
    b = [0;ones(3,1);a1;a2;a3;a4];
    c = [1;mean(vertices(:,1));mean(vertices(:,2));mean(vertices(:,3));mean(vertices(:,1).*vertices(:,2));
        mean(vertices(:,1).*vertices(:,3));mean(vertices(:,2).*vertices(:,3));
        mean(vertices(:,1).*vertices(:,2).*vertices(:,3))];
    B = b*b';
    kj = J'*B*J*V*cond;   % 试函数的转置与试函数的乘积再乘以element的面积就是K
    fj = I*J'*c*V;
end       
end

下面为对应函数调用调用代码：

cfg = [];
cfg.conductivity = [1.79 0.33 0.43 0.01 0.14];
[K,u,rcondj] = FEM_hexahedral(cfg, mesh);  % K为刚体矩阵，rcondj为内部构造形函数相应的矩阵
[row,col,v] = find(K);
figure(1)
plot(u)
title("对应于坐标序号的电压分布")
figure(2)
plot3(row,col,v,'.')
title("刚体矩阵的分布")

 mesh数据为：fieldtrip工具箱“subject01”被试的网格数据——降采样设置为10。-其它文档类资源-优快云下载