Dirichlet boundary condition for sparse matrix - Solving Ax=b only for free nodes?

Question

I am solving Biot equation with sparse matrix in MATLAB. I have no problem with global sparse matrices assembly, but when I assign Dirichlet boundary condition, it is so slow.

From this topic, How to efficiently implement Dirichlet boundary conditions in global sparse finite element stiffnes matrices, @James also mentioned that it is so slow if global matrix is modified. My matrix is 25000x25000, and it took around 100s to assign boundary conditon.

In MATLAB, many authors state that, by using setdiff function, only free nodes are considered when solving the equation. Only right hand vector is modified when applying dirichlet boundary condition and global matrix is kept unchanged.

Hier is the code in MATLAB:

% Dirichlet conditions 
u = sparse(size(coordinates,1),1);
u(unique(dirichlet)) = u_d(coordinates(unique(dirichlet),:));
b = b - A * u;

% Computation of the solution
u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes);

What is the difference between two approaches? And is there any similar function as setdiff in MATLAB in C++?

Answer 1

I frequently solve simple FEM problems with 100 000 to 1 000 000 unknowns in Matlab and setting Dirichlet boundary conditions is blazing fast.

Here is a code (although you cannot run it without the related functions but you get the idea):

% 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);

The whole thing took me 3.5 seconds to run with n=263169. The result looks like this:

One cannot really see what's happening without removing the element edge lines:

What the command

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

actually does is that it forms a smaller sparse matrix by including only the columns and rows defined in set I. In fact, running

tic;K(I,I);toc

prints

Elapsed time is 0.018481 seconds.

so this is a very fast operation indeed.

Let's say that in C++ you store your sparse matrix in compressed row storage format. Then you can easily drop rows by looking at the contents of your row index array. Next perform a transpose and drop rows again. This way you end up with a smaller sparse matrix that you can use in your computations.

Edit: Forgot to mention, that in C++ standard library you have set_difference as well.

Answer 2

Just see this question. The way I did it in c++ is that instead of implementing the Dirichlet Condition afterwards, I do not put the element into the row/column where Dirichlet Condition apply.

This is just a condition loop inside the main loop when you do the triplets of each element, based on either the index belongs to those row/columns.

And you can use a hash-map (unsorted_set/unsorted_map) to store the fixed nodes. Total cost is still O(1).