Consider the linear system of equations
A*x = b
We consider efficient application of the Dirichlet boundary conditions, i.e., the x value in the boundary nodes is explicitly defined. Assume node with index k is a boundary node. To set the value of x(k)=1 the k‘th row of the system matrix A is modified in the following way
A(k,i)=0 for i=1:n and i!=k
b(k) = b(k) – A(i,k)*x(k)
A(i,k) = 0
Free_ind = 1:n; Free_ind(Bc_ind) = ; % for all boundary nodes for id=1:length(Bc_ind) % deal with equation k k = Bc_ind(id); % eliminate the k'th variable rows = A(Free_ind,k); b(Free_ind) = b(Free_ind) - rows*Bc_val(id); % modify A and b A(:,k) = 0; A(k,:) = 0; A(k,k) = 1; b(k) = Bc_val(id); end
This implementation is very inefficient since every boundary node is processed independently. In particular, reseting row (and column) entries of a sparse matrix to zero introduces costly changes to the sparsity pattern. Avoiding multiple modifications of the sparsity pattern of A is expected to minimize the related overhead.
The code fragment below presents two different approaches to apply the Dirichlet boundary conditions in MATLAB efficiently. Instead of modifying matrix entries, the rows and columns corresponding to the boundary nodes can be all together removed from the system. The update to the right hand side vector is recast into a (dense) matrix-vector multiplication. For symmetric matrices only the upper- or lower-triangular part of A needs to be stored and the implementation is modified accordingly.
Free = 1:n; Free(Bc_ind) = ; b = zeros(n,1); b(Bc_ind) = Bc_val; % Approach Ia: b = b - (A(:,Bc_ind) + cs_transpose(A(Bc_ind,:)))*Bc_val'; % Approach Ib: b = b - (A(:,Bc_ind) + (A(Bc_ind,:))')*Bc_val'; % Approach IIa: b = b - A*x - (x'*A)'; % Approach IIb: b = b - A*x - A'*x; % Elimination of trivial equations 1 A = A(Free_ind,Free_ind); % Elimination of trivial equations 2 A(:,Bc_ind) = ; A(Bc_ind,:) = ;
Execution time of the above implementations on an example desktop architecture is compared in Figure TODO. The system size is approximately one million degrees of freedom. The number of Bc_ind nodes in the system varies. When the fraction of the boundary nodes is smaller than 1% we observe no substantial difference in the performance of the four methods (Ia,Ib,IIa,IIb) to modify the right-hand side vector b. Most practical two-dimensional models are in this regime. For higher fractions of up to 40% the approaches IIa and IIb perform equally well as for the lower fractions. Implementations Ia and Ib are significantly less efficient in this regime.
Figure TODO presents the execution time of the two studied methods to eliminate the the rows and columns corresponding to the boundary nodes from the system matrix A. When the fraction of the boundary nodes is smaller than 1% the Elimination 1 is slightly more efficient than Elimination 2. For higher fractions of boundary nodes Elimination 1 becomes increasingly more efficient, while the performance of Elimination 2 remains unchanged.
In summary, approaches IIa or IIb should be used to modify the right-hand side vector b, and Elimination 1 should be used to removing the trivial equations for the boundary nodes Bc_ind from the system matrix A.