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

Question

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!

Answer 1

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.

Answer 2

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.