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 solutionu(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++?
Dirichlet boundary condition for sparse matrix - Solving Ax=b only for free nodes?
finite element;boundary conditions;sparse
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 meshmesh=squaremesh(9);% build stiffness matrixK=bilin_assembly(@(u,v,ux,uy,vx,vy,x,y) ux.*vx+uy.*vy,mesh);% build load vectorf=lin_assembly(@(v,vx,vy,x,y) 1.*v,mesh);% number of nodesn=size(mesh.p,2);% nodes on the dirichlet boundaryD=boundary_dofs(mesh);% free nodes are all nodes except dirichlet boundary nodesI=setdiff(1:n,D);% initialize u as zerou=zeros(n,1);% set dirichlet conditionX=mesh.p(1,D);u(D)=sin(10*X);% solve the problem in the free nodesu(I)=K(I,I)\(f(I)-K(I,D)*u(D));% plot the solutionfigure(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 commandu(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, runningtic;K(I,I);tocprintsElapsed 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.