# Do I need to impose boundary conditions in the Jacobian matrix?

In the framework of Finite Element Method, when the Newton method is used, we solve $J(x^k) \delta x = -f(x^k)$, and the increment $\delta x$ would not change some entries from $x^k$ related to Dirichlet boundary conditions. So, I was wondering about the needs to impose homogenous dirichlet conditions in the Jacobian matrix $J(x^k)$, such that $\delta x$ would have some zero entries. My question are: is this really need? Should I leave the Jacobian untouched and post-process $x^{k+1}$?

Thanks!

• @GeoffOxberry . Does it mean that if the Jacobian and internal force is calculated under the consideration of Dirichlet boundary and then the RHS of those non Dirichlet node don't need to be modified even though i zero out the Dirichlet row and set the RHS related to Dirichlet boundary to zero? Thanks. Apr 10, 2016 at 18:20

This remark is a bit too long for a comment, but regarding DanielShapero's remark with respect to changing the Jacobian matrix, for every degree of freedom $i$ corresponding to a Dirichlet boundary condition (BC), a common approach is to set $i$th row of the Jacobian matrix equal to the $i$th row of the identity matrix, meaning your (quasi-)Newton iterates should be linear systems that can be permuted into the form

\begin{align} \left[\begin{array}{cc}A & B \\ 0 & I\end{array}\right]\left[\begin{array}{c}\delta x_{interior} \\ \delta x_{BC}\end{array}\right] = \left[\begin{array}{c}-f_{interior}(x^{k}) \\ -f_{BC}(x^{k})\end{array}\right]. \end{align}

As long as the initial guess to the nonlinear solve sets the Dirichlet BC degrees of freedom to the values specified in those boundary conditions, then at every iteration of the nonlinear solver, the Dirichlet boundary conditions should be satisfied, and $\delta{x}_{BC} = 0$ because $-f_{BC}(x^{k}) = 0$.

You can do something similar for functional iteration-type nonlinear solvers.

EDIT: If you zero out the rows of the permuted Jacobian matrix to get

\begin{align} \left[\begin{array}{cc}A & B \\ 0 & 0\end{array}\right]\left[\begin{array}{c}\delta x_{interior} \\ \delta x_{BC}\end{array}\right] = \left[\begin{array}{c}-f_{interior}(x^{k}) \\ -f_{BC}(x^{k})\end{array}\right], \end{align}

then the resulting rank-deficient system can be solved using least-squares direct methods (QR, SVD) or Krylov methods (which will calculate a Drazin inverse and converge to the correct solution). LU-based methods should return an error due to a zero pivot, and can't be used without reformulating the system in terms of the interior degrees of freedom only, which is a minor drawback of zeroing out the rows. Solving the reformulation $A \delta{x}_{interior} = -f_{interior}(x^{k})$ is probably preferable at that point.

• Thanks, @GeoffOxberry the comprehensive answer. Regarding the consequences of this approach, did you face any problem related to load balance, communication etc. in parallel architectures? E.g. when dealing with distributed unstructured meshes, this permutation you suggested can be expensive, in terms of communications. Oct 9, 2014 at 22:13
• Besides, PETSc function 'MatZeroRows' is very convenient. (sorry the many 'edits', I was not aware the ENTER key publishes the comment). Oct 9, 2014 at 22:25
• @FredericoTeixeira: The permutation is unnecessary for calculations, and should not be done in production, so there should not be issues due to communication or load balancing. I only use permutation in the explanation to clarify the notation so that the matrices reduce to block form. You can use MatZeroRows to set the diagonal entries to 1 and zero the remaining elements; just set the diag argument to 1. Oct 9, 2014 at 22:40

There are multiple ways of implementing Dirichlet boundary conditions when using the finite element method.

For each unknown $i$ of the system belonging to the Dirichlet boundary, you can zero out row $i$ of the Jacobian and entry $i$ of the right-hand side. At that rate, you can reformulate the problem for just the unknowns in the interior, which reduces its size.

Alternatively, you can enforce the Dirichlet boundary conditions at every step of the nonlinear solve. However, this assumes that you're the one responsible for writing the solver; if you want to use someone else's software to solve the nonlinear system for you, you may have to resort to the first approach.

This book has a whole section on Dirichlet boundary conditions which covers both of these options in more detail and a few others besides, e.g. penalty methods.

• I don't agree with the construction you propose in your second paragraph, and I've suggested an alternative. It was too long to post as a comment, so I wrote it as an answer. Oct 9, 2014 at 18:45
• I usually use the approach you outline in your answer, but the OP formulated the problem as solving for an update vector $\delta x$ rather than solving directly for the new iterate $x^{k + 1}$. Oct 9, 2014 at 20:27
• @DanielShapero, thanks your comments and suggestions! The formulation in terms of $\delta x$ is very convinient for PETSc non-linear solvers. Oct 9, 2014 at 22:24