我怀疑你会失望地发现，一旦你能够实现二次形函数，你得到的解决方案是不正确的。这是因为对于 bilaplacian 算子 ( )，在每个单元格上具有二次多项式的形状函数是不够的，但您还需要具有跨单元格没有“扭结”的特殊类型 (Hermite) 形状函数接口，但在面处是连续可微的。这些比您正在寻找的二次形状函数更难实现，但对于产生您试图求解的方程的正确解是必要的。$\nabla^4$

这是我想到的版本。它有点简洁，但在 github 上有一个注释版本（使用那个，它实际上定义了输入/输出并解释了各种中间计算的逻辑）。基于一些初始基准测试，它比你的循环/分支繁重的代码高出大约 300 倍（Matlab 2017a 在一个大的网格上），因为它把所有时间都花在了矢量化调用上。

这是真实代码的链接： https ://github.com/tjolsen/Mesh_Utilities/blob/master/linear_to_quad_tris/linear_to_quadratic.m

function [QuadTris,Points] = linear_to_quadratic(t,p)

nt = size(t,1);
t_idx = (1:nt)';

all_edges = [t(:,1) t(:,2) t_idx t_idx 3*ones(nt,1);
t(:,2) t(:,3) t_idx t_idx 1*ones(nt,1);
t(:,3) t(:,1) t_idx t_idx 2*ones(nt,1)];

forward_edges = all_edges(:,2) > all_edges(:,1);
all_edges(forward_edges,4) = 0;
all_edges(~forward_edges,3) = 0;
all_edges(~forward_edges,[1,2]) = all_edges(~forward_edges,[2,1]);

[unique_edges,~,ic] = unique(all_edges(:,[1,2]), 'rows');
Edges = [unique_edges, 
         accumarray(ic, all_edges(:,3),[],@sum),
         accumarray(ic, all_edges(:,4),[],@sum)];
t2f = accumarray([all_edges(:,3)+all_edges(:,4),all_edges(:,5)], 
                 sign(all_edges(:,3) - all_edges(:,4)).*ic, [], @sum);

newPoints = [(p(Edges(:,1),1)+p(Edges(:,2),1))/2, 
             (p(Edges(:,1),2)+p(Edges(:,2),2))/2];
Points = [p; newPoints];

np = size(p,1);
QuadTris = [t, np+(abs(t2f))];

end

遵循@TylerOlson 的建议，这里是一个 Octave/Matlab 脚本，它生成一个在每个三角形中有六个点的网格：三个顶点和三个中点。如果有人知道简化这种不雅代码的方法，我将不胜感激。

clear
%% Persson's routine to make simple triangles in a rectangular grid
nx=4;ny=3;
[x,y]=ndgrid(linspace(-1,1,nx),linspace(-1,1,ny)); % forms x and y lists
p=[x(:),y(:)]; 
t=[1,nx+2,nx+1;1,2,nx+2]; 
t=kron(t,ones(nx-1,1))+kron(ones(size(t)),(0:nx-2)');
t=kron(t,ones(ny-1,1))+kron(ones(size(t)),(0:ny-2)'*nx);
np=length(p(:,1));nt=length(t(:,1));
%%make a list of all sides from t;  There will be redundant entries
s=zeros(3*nt,6);
for kt=1:nt
  for ks=1:3
    if ks==1
      s(3*(kt-1)+ks,1)=min(t(kt,1:2));
      s(3*(kt-1)+ks,2)=max(t(kt,1:2));
      s(3*(kt-1)+ks,3)=kt;
    elseif ks==2
      s(3*(kt-1)+ks,1)=min(t(kt,2:3));
      s(3*(kt-1)+ks,2)=max(t(kt,2:3));
      s(3*(kt-1)+ks,3)=kt;
    else
      s(3*(kt-1)+ks,1)=min(t(kt,[1,3]));
      s(3*(kt-1)+ks,2)=max(t(kt,[1,3]));
      s(3*(kt-1)+ks,3)=kt;
    endif
  endfor
endfor
%% once this has been done for all triangles, work on s to eliminate dupes
for kk=1:max(s(:,1))
  order1=find(s(:,1)==kk);%s(order1,:)
  [junk,order2]=sort(s(order1,2));
  dupe=find(junk(1:end-1)==junk(2:end));
  n=s(order1(order2),:);
  for kl=1:length(dupe)
    n(dupe(end+1-kl),4)=n(dupe(end+1-kl)+1,3);% add on triangle no of dupe
  endfor
  [junk,order3,mojunk]=unique(n(:,2),"first");
  n=n(order3,:);
  if kk==1
    snodupe=n;
  else
    snodupe=[snodupe;n];
  endif
endfor
%% add to position table the position of new nodes
for kn=1:length(snodupe)
    p(np+kn,1)=0.5*(p(snodupe(kn,1),1)+p(snodupe(kn,2),1));
    p(np+kn,2)=0.5*(p(snodupe(kn,1),2)+p(snodupe(kn,2),2));
endfor
%% build a new connectivity table, t6
t6=zeros(nt,6);
for kt=1:nt
  %t(kt,:)
  for ks=1:3
    if ks==1
      ss=t(kt,1);
      se=t(kt,2);
    elseif ks==2
      ss=t(kt,2);
      se=t(kt,3);
    else
      ss=t(kt,3);
      se=t(kt,1);
    endif
    %% where is this side in snodupe table?
    %% to lookup in snodupe need ends to be in order
    s1=min([ss se]);s2=max([ss se]);
    ind1=find(snodupe(:,1)==s1);
    ind2=find(snodupe(ind1,2)==s2);
    sideno=ind1(1)-1+ind2;
    t6(kt,2*ks-1)=ss;
    t6(kt,2*ks)=np+sideno;
  endfor
endfor