我经常在 Matlab 中解决具有 100 000 到 1 000 000 个未知数的简单 FEM 问题，并且设置 Dirichlet 边界条件非常快。

这是一个代码（虽然你不能在没有相关功能的情况下运行它，但你明白了）：

% get mesh
mesh=squaremesh(9);

% build stiffness matrix
K=bilin_assembly(@(u,v,ux,uy,vx,vy,x,y) ux.*vx+uy.*vy,mesh);

% build load vector
f=lin_assembly(@(v,vx,vy,x,y) 1.*v,mesh);

% number of nodes
n=size(mesh.p,2);

% nodes on the dirichlet boundary
D=boundary_dofs(mesh);

% free nodes are all nodes except dirichlet boundary nodes
I=setdiff(1:n,D);

% initialize u as zero
u=zeros(n,1);

% set dirichlet condition
X=mesh.p(1,D);
u(D)=sin(10*X);

% solve the problem in the free nodes
u(I)=K(I,I)\(f(I)-K(I,D)*u(D));

% plot the solution
figure(1);
trisurf(mesh.t',mesh.p(1,:),mesh.p(2,:),u);

在 n=263169 的情况下，整个过程花费了我 3.5 秒的时间。结果如下所示：

如果不删除元素边缘线，就无法真正看到正在发生的事情：

什么命令

u(I)=K(I,I)\(f(I)-K(I,D)*u(D));

实际上，它通过仅包含集合 I 中定义的列和行来形成一个更小的稀疏矩阵。事实上，运行

tic;K(I,I);toc

印刷

Elapsed time is 0.018481 seconds.

所以这确实是一个非常快速的操作。

假设在 C++ 中，您以压缩行存储格式存储稀疏矩阵。然后，您可以通过查看行索引数组的内容轻松删除行。接下来执行转置并再次删除行。这样，您最终会得到一个较小的稀疏矩阵，您可以在计算中使用它。

编辑：忘了提，在 C++ 标准库中你也有set_difference。

看看这个问题。我在 C++ 中这样做的方式是，我没有在之后实现 Dirichlet 条件，而是不将元素放入应用 Dirichlet 条件的行/列中。

这只是主循环内的一个条件循环，当您根据索引属于那些行/列时对每个元素执行三元组时。

您可以使用哈希映射（unsorted_set/unsorted_map）来存储固定节点。总成本仍然是 O(1)。