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

A(k,k)=1

*k*‘th row is set to one, and all other entries are set to 0. The

*k*‘th entry of the right-hand side vector b is set to 1

b(k)=1

*k*‘th equation is trivial and reads

1*x(k)=1

*k*‘th matrix row have been explicitly set to 0. The symmetry can be regained by setting all the rows in the

*k*‘th matrix column to 0. This is done by eliminating the known x(k) from all the equations it is part of. For every equation

*i*, for which A(i, k) is non-zero the right-hand side vector entry b(k) is modified as

b(k) = b(k) – A(i,k)*x(k)

A(i,k) = 0

*Bc_ind*and boundary values

*x(Bc_ind) = Bc_val*. We denote the indices of the unconstrained nodes as Free_ind. The boundary conditions can be implemented by directly applying the above described modifications to A and b for every boundary node separately.

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.

[vslider name=”BCtest”]