Dirichlet condition in finite element method

Question

I'm trying to understand the approach described in these questions:

How to efficiently implement Dirichlet boundary conditions in global sparse finite element stiffnes matrices

How to apply non zero Dirichlet boundary condition in finite elements?

While I understand what's happening, I don't understand why the approach is justified.

My understanding of what's happening:

Throughout this post, I will use the notation of $a_{0,..,n}$ denoting the $n \times 1$ matrix of elements $(a_0, ..., a_n)$.

Here is what I see as the a description of the actual solution - suppose we have the system $Ax = 0$ where $A$ is a $n \times n$ matrix, and suppose $x_0 =1 $ by the boundary condition, and all the other $x_i$ are free. Denote $A$ as \begin{bmatrix} c_0 & c_{1,...,n}^T\\ c_{1,...,n} & K \end{bmatrix}

Then $A(x_0e_0 + \Sigma_{i>0}x_i e_i) = A(x_0e_0) + A(\Sigma_{i>0}x_i e_i) = 0$ and so $A(\Sigma_{i>0}x_i e_i) = -A(x_0e_0) = -A(e_0)$, in other words: \begin{equation*} \begin{bmatrix} c_{1,...,n}^T\\ K \end{bmatrix} * (x_{1,...,n}) = -(c_{0,...,n}) = - \begin{bmatrix} c_0\\ c_{1,...,n} \end{bmatrix} \end{equation*} and so we essentially end up solving a system $\tilde{A}\tilde{x} = \tilde{b}$ where $\tilde{A}$ is an $n \times (n-1)$ matrix.

However, the approaches in the linked questions reach a different equation: \begin{equation*} K * (x_{1,...,n}) = -(c_{1,...,n}) \end{equation*}

My question:

Essentially, in the linked approaches, the first row is zeroed out, and I don't quite understand how that makes sense. You can zero the first column, sure, as long as you move it to the right hand side, but I'm not sure how it makes sense to zero out the whole row too - if there are multiple solutions to $K * (x_{1,...,n}) = -(c_{1,...,n})$, we need to find the one that ends up giving $-(c_{0,...,n})$.

How is this step justified?

Answer 1

From a different perspective, regardless of where a system $A\cdot x = 0$ comes from (does not have to be FEM), if you change your mind and would like to prescribe part of $x$, then you'd necessarily also have to introduce a new unknown reaction "force" $f_0$ to compensate: $$ \begin{bmatrix}A_{00} & A_{f0}\\A_{0f} & A_{ff}\end{bmatrix} \cdot \begin{bmatrix}x_0 \\ x_f\end{bmatrix} = \begin{bmatrix}f_0 \\ 0\end{bmatrix} $$

You just don't get to decide both the LHS and RHS.

One could from there just algebraically rearrange the system to actually solve for the free unknowns $x_f$ and $f_0$

$$ A_{ff} \cdot x_f = -A_{f0}\cdot x_0 \\ f_0 = A_{00} \cdot x_0 + A_{0f} \cdot x_f $$

None of $A$ is zeroed out here, but if you don't care about $f_0$, you of course never need to use $A_{00}$ or $A_{0f}$.

Some simple FEM codes often used in teaching often employ this algebraic approach for solving linear problems (and computing the reaction forces as a bonus).

I wouldn't recommend this approach for real production codes, as assembling only the actual $A_{ff}$ directly is cheap and easy, and computing the reaction forces $f_0$ just comes as part of computing the residual in a nonlinear problem anyway.

Answer 2

From the finite element theory, if you specify Dirichlet boundary condition on part of your domain, then the test function is zero on that part. This removes the effect of any force applied both surface and body forces as they are multiplied with the test function in the weak form and hence there is no force pulling the nearby nodes from the Dirichlet boundary, but the displacement specified or some value of your function specified will only impact the nearby nodes.

Putting it other way from a mechanical point of view, let us say you are pulling a 1D bar, on the one end you specify a Dirichlet boundary condition of 1mm and you also specify a force of 10N on that same node on that end. In this case the displacement boundary condition will be the only factor impacting the displacement of the nearby nodes and the force will have no effect because it is an over constraint and hence can be neglected.

From the linear algebra point of view, there are another set of vectors in your system of equations, the test function, let me denote $y_{0,...,n}$ as the test function vector as opposed to the trial function $x_{0,...,n}$. From the weak form you will end up with this set of equations,

$\begin{bmatrix} y_0 & y_{1,...,n} \end{bmatrix}^\top\begin{bmatrix} c_0 & c_{1,...,n}^T\\ c_{1,...,n} & K \end{bmatrix} \begin{bmatrix} x_0 \\ x_{1,...,n} \end{bmatrix} = 0$

Hence if our test function is zero $(y_0 = 0)$, then it multiplies to the top row with zero and hence in the resulting system of equations we can neglect that row.

The test function is also described as perturbation $(\delta x)$ to the unknown field that we solve in some engineering courses. And hence it makes sense to say that the perturbation to our unknown is zero where we know the unknowns (they are not unknowns anymore)

Answer 3

After some studies, I want to provide an example of a case where it can be shown that the submatrix $K$ as described in my question is actually positive-definite, not just positive-semidefinite.

Consider the case of solving the heat equation/Poisson's equation on a connected simplicial complex. Then the problem can be described using $Lx = b$, where $L$ is the discrete Laplacian.

It can be shown that $L$ has a one-dimensional kernel, formed by the vector $(1,1,...,1)$. This corresponds to the fact that a solution to the heat equation without a boundary condition is unique up to a constant.

Now using Kirchhoff's theorem, we can show that by removing any row and column, we acquire a matrix (this corresponds to the $K$ submatrix) that has a non-zero determinant, as we assume that the mesh is connected, and so some spanning tree exists. In other words, removing a row and column gives us a positive-definite matrix. This means the equation $K * (x_{1,...,n}) = -(c_{1,...,n})$ has a unique solution.

If we remove more rows and columns from a positive-definite matrix, then by Cauchy's interlacing theorem such a matrix remains positive definite. Thus the approach of removing $m$ rows and columns and solving the $(n-m) \times (n-m)$ and then solving for a right-hand side of size $(n-m)$ is appropriate because there is only a unique solution anyway.

---

Whether such an approach/explanation could be applied to other FEM problems that are not described by the Laplacian, I do not know.