# What Linear Equation Solver should be used for a problem with many dirichlet conditions?

I am solving a laplace equation on a finite-element mesh (tetrahedral, triagonal) and have many say 99% dirichlet conditions compared to the number of unknowns. Is there an efficient way to solve this kind of system?

To give an example for better understanding: Total number of nodes: 100,000 Nodes Dirichlet: 99,000 Number of Unknowns not Dirichlet: 1,000

It can be direct or iterative that doesnt matter.

It should be fast and efficient though.

Edit 1: The domain can be regular (box shaped) or irregular (eg. irregular shaped surface in 3D)

Edit 2: The dirichlet conditions are can be on the boundary and/or inside the domain as a single dirichlet condition is enough to make the solution unique.

• Are you saying that 99% of your nodes are situated on a boundary where Dirichlet conditions hold? – Wolfgang Bangerth Jul 25 '13 at 13:09
• How many unknowns are left after the elimination of Dirichlet conditions? What is the shape of the domain? – Guido Kanschat Jul 25 '13 at 13:09
• @Wolfgang Bangerth No, it is sufficient to have 1 dirichlet support with a domain to result in a unique solution. The dirichlet nodes are on the boundary as well as inside the domain. – Lukas Mosser Jul 25 '13 at 13:13
• @Guido Kanschat How would you eliminate the dirichlet conditions? E.g. Total number of nodes 100 000. Of those dirichlet: 100. Domain shape can be box to irregular. – Lukas Mosser Jul 25 '13 at 13:13
• @LukasMosser AMG is likely overkill for this problem. You can always lift the Dirichlet boundary conditions to the right hand side of a reduced system. You can derive this by (on paper) ordering unknowns so that all the Dirichlet entries come first in the 2x2 block system $\begin{pmatrix} I & 0 \\ B & A \end{pmatrix} \begin{bmatrix} x_D \\ x_I \end{bmatrix} = \begin{bmatrix} b_D \\ b_I \end{bmatrix}$, which becomes $A x_I = b_I - B b_D$. – Jed Brown Jul 26 '13 at 15:00

The best thing to do here is to actually reduce the size of your problem after you have assembled the full system matrices. Let's assume you've gotten to that point, so you have

$$\mathbf{\tilde{K}d}=\mathbf{\tilde{F}}$$

where $\mathbf{K}$ is a $N_{dof} \times N_{dof}$ matrix of stiffness coefficients, $\mathbf{d}$ is your full vector of nodal degrees of freedom (most of which are known in your case), and $\mathbf{F}$ is a $N_{dof} \times 1$ vector of "forces" which arose from integration of the source terms and Neumann boundary conditions.

Rather than solving this very large system, you want to reduce the system down to ONLY the unknowns in $\mathbf{d}$! The method I will describe here is very general, and is not limited to specially-shaped domains.

To do this, first (for the sake of simplicity) that the "top" half of your $\mathbf{d}$ vector contains your known values from the Dirichlet boundary conditions, and the bottom half contains your unknowns. Immediately, we see that we can simply ignore your first $N_{dof}/2$ equations. In practice, this means chopping off the first $N_{dof}/2$ rows of $\mathbf{\tilde{K}}$ and $\mathbf{\tilde{F}}$. Now we're left with a more manageable set of equations

$$\mathbf{\hat{K}d}=\mathbf{\hat{F}}$$

where $\mathbf{\hat{K}}$ is a $\frac{N_{dof}}{2} \times N_{dof}$ matrix, and $\mathbf{\hat{F}}$ is a vector of length $N_{dof}/2$. We haven't yet done anything with $\mathbf{d}$.

Now, we're going to reduce $\mathbf{d}$ down to only having the unknowns on the left side of the equation. Notice that the first $N_{dof}/2$ columns of $\mathbf{\hat{K}}$ are being multiplied with known values of $\mathbf{d}$. We will take advantage of this fact and move all of the known values to the right-hand side of the equation. That statement may have been a little cryptic, but the best way to illustrate it is with a little bit of code. I'm going to use a MATLAB-like syntax and take advantage of matrix slicing.

for ii=1:Ndof                            //loop through d
if(isDirichletDOF(ii))               //proceed if d(ii) is known
Fhat = Fhat - d(ii)*Khat(:,ii);
end
end


In the small snippet above, I only moved over columns associated with KNOWN values (ie, all of the dirichlet boundary values) and did not touch any of the others. Now that you've done this, you no longer need to worry about the known values of $\mathbf{d}$ any more, so just chop it down to a vector of only unknowns. Simultaneously, make a new matrix $\mathbf{K_{final}}$ from $\mathbf{\hat{K}}$ that uses only the columns that were not associated with known values. Now you finally have the system that you want to solve.

$$\mathbf{K_{final}d_{unknowns}} = \mathbf{\hat{F}_{modified}}$$

In your particular application, this will be MUCH MUCH faster than some other methods of applying boundary conditions. From here, you pick your favorite linear solver and go to town. I won't try to predict actual performance gains since that depends so much on the linear solver you use, but $\mathbf{K}$ is now $1/10,000$ the size that it used to be (if, for some reason you were not using a sparse matrix).

Obviously, in practice, you don't know your node numbering ahead of time (for non-trivial geometries) so you just need to do some bookkeeping to determine which rows to get rid of and which columns to move to the right hand side of the equation.