# Why is the test function space in FEM chosen with homogeneous boundary conditions?

It is so confusing, especially when I learns discontinuous galerkin method in broken Sobolev space and weak Dirichlet boundary condition.

If the trial function is chosen with homogeneous boundary conditions, how can I get a linear discrete equation corresponding to test basis functions with DOFs on the boundary? I mean if $v\in V_0 \subset V$, then the space of $u\in V$ would be "larger" than $V_0$. The final matrix would be rectangular rather than square. Extra conditions must be imposed to make it square and solvable.

My understanding now is that the homogeneous boundary condition is helpful in derivation of weak form from original differential equation, such as Poisson equation.

edit: When it comes to matrix, a $V_{0,h}\subset V_0$ and $V_{h} \subset V$ must be chosen first.

edit2: I am confused, sometimes the book may say bases of $V_0$ are identical to $V,$ but the boundary condition is different.

This is more of a conceptual argument but may be helpful in understanding why the math works out this way. As a bonus, no matrix will be involved (because the problem is actually already with the weak formulation before any discretization).

What FEM does is to numerically approximate the solution $y\in V$ to some abstract operator equation $Ay=f$ for $A:V\to V^*$ for $f\in V^*$ and some abstract vector space $V$ and its dual $V^*$ (it's not really necessary for this argument to understand what a dual space is; it's enough to accept that it's the space for which the following works).

The first step is to see that $Ay-f = 0$ if and only if $\langle Ay-f,v\rangle = 0$ for all $v\in V$. (Intuitively, $0$ is the only vector that is orthogonal to every other vector.) Note that it's the same space $V$ for $y$ and $v$. It is crucial for this equivalence that $V$ is a vector space, in particular, that for $v_1,v_2 \in V$ you also have $v_1+v_2\in V$.

The next step is specific for PDEs. If $A$ is a differential operator such as the Laplacian $-\Delta$, it's easy to see that the solution $y$ is not unique unless you fix constants. The usual (but not the only) way is to fix the constant by fixing the value at the boundary. To include this in the above approach means only looking in the space of all functions that have the same values at the boundary (both for $y$ and for $v$, remember!) Now it should be fairly obvious that this doesn't work unless that value is $0$, because otherwise if you add two functions with the same boundary values, their boundary value will be different.

Of course, that's a too nasty restriction to take lying down. Luckily there's a trick: Pick a function $y_b$ that has the desired values at the boundary and satisfies $A y_b=0$, compute a solution $y_0$ to $Ay=f$ with zero boundary values as described above (including using test functions $v$ with zero boundary values), and set $y=y_0+y_b$. This is precisely what you are doing in computations, with some short cuts to avoid explicitly constructing $y_b$ first.

You mention the weak Dirichlet boundary conditions: Here you do not prescribe boundary conditions, but "penalize" non-zero values so that the solution will actually be (close to) zero. Since you do not have boundary conditions for the solution, you also shouldn't have boundary conditions for the test functions. (And for DG, there's additional difficulties since the weak form is an operator equation $A:V\to W^*$ for $W\neq V$, so the test functions have to come from a different space of functions that don't even need to have sensible values on the boundary.)

• Thank you for your answer. How about Nitsche's method? I think computer is using $a(u,v)=f(v)$ to generate matrix $A$ and source vector $b$ by $A_{ij} = a(\psi_i,\phi_j)$ and $b_i = f(\phi_j)$ for $\psi_i$ is basis of $V_h\subset V$ and $\phi_j$ is basis of $W_h\subset V$. If test function space $W_h$ is limited to homogeneous BC, its dimension would be lower than $V_h$. Additional Dirichlet BC must be imposed to make the matrix square. If Nitsche's method is used. No extra BC need to be imposed after the matrix assemble. – CatDog Aug 3 '18 at 3:14
• For DG, you mean DG is a special kind of PG? I thought DG can be adjoint consistent because the test and trial function spaces are identical. I thought PG is of adjoint consistency problem because its test function space is discontinuous while its trial function space is continuous. – CatDog Aug 3 '18 at 3:14
• As I wrote, you do not impose any boundary conditions (homogeneous or otherwise) in Nitsche's method; see my answer scicomp.stackexchange.com/a/19913. I don't know what PG means. (If it's Petrov-Galerkin, then yes. But for these methods test and trial space are by definition not the same; that's also not what adjoint consistent means -- if you're unsure about this point, you should ask a new question.) – Christian Clason Aug 3 '18 at 6:26
• Thank you for your answer. PG = Petrov-Galerkin. I just read some thesis and paper on adjoint consistency because of its potential application and relation to optimal convergence. But it is so difficult because of the math behind it. – CatDog Aug 3 '18 at 15:25

The important point, not mentioned in the answers above, is that the weak form is a variational formulation in which we ask the question: Can we find a function $u$ that minimizes the energy functional (or is a stationary point of some other functional) in such a way that any other potential solution $u+v$ does not minimize the energy (or is not a saddle point). The $v$ here becomes the test function in the weak formulation.

If we seek solutions among the set of functions that have boundary values $u|_\Gamma=g$, then clearly the perturbed functions also have to satisfy $(u+v)|_\Gamma=g$, which can only be the case if $v|_\Gamma=0$, i.e., satisfies homogeneous Dirichlet boundary conditions.

A lengthier discussion of this issue can be found in video lecture 21.5 and following here: http://www.math.colostate.edu/~bangerth/videos.html

• Thanks! This is a very straightforward explanation if the functional exist. $v$ is just a direction in the sense of direction derivative. – CatDog Aug 3 '18 at 22:11
• That's exactly right. – Wolfgang Bangerth Aug 4 '18 at 14:03

While it is true that $V_0 \subset V$, these spaces are still both infinite dimensional, and there are no matrix sizes to speak of since no discretization has happened yet. To form matrices, you must choose finite dimensional subspaces $V_{0,h}$ and $V_h$ which "approximate" $V_0$ and $V$, in some sense. The size of the resulting matrices just depends on the dimension of $V_{0,h}$ and $V_h$. Therefore, if you use the same number of basis functions so that your finite-dimensional approximate spaces $V_{0,h}$ and $V_h$ have the same dimension, your matrices will end up square.

EDIT: Actually, the last sentence above is incorrect. You can get square matrices without having the exact same dimension for your subspaces; what matters is obviously the number of DOF you have, which also includes essential boundary conditions which are enforced in your problem. The following example should hopefully be more helpful than what I wrote above.

We seek a function $u(x)$ such that \begin{align} -u'' +u = f, \: x \in \Omega, \\ u'|_{\Gamma_N} = g_N, \: u|_{\Gamma_D} = g_D, \end{align} where $\Omega$ is some open interval whose boundary is partitioned $\partial \Omega = \Gamma_N \cup \Gamma_D$. One weak form for this problem is: find $u \in U$ such that \begin{align*} \int_{\Omega} v' u' + \int_{\Omega} vu = \int_{\Omega} vf - v|_{\Gamma_N} g_N , \: \forall v \in V, \end{align*} where \begin{align} V = \{ w \: | \: w \in H^1(\Omega), w|_{\Gamma_D} = 0 \}, \\ U = \{ w \: | \: w \in H^1(\Omega), w|_{\Gamma_D} = g_D \}. \end{align} Note that we set $w|_{\Gamma_D} = 0$ in the test function space $V$ to eliminate the $v|_{\Gamma_D} u''|_{\Gamma_D}$ term since no information about $u''$ is known on that part of the boundary. This does not affect the DOF boundary conditions which we will enforce later, since we can enforce them directly without using our test function space $V$.

Lets proceed by a Petrov-Galerkin method: we make finite-dimensional subspaces by specifying bases \begin{align*} U_m = \text{span} \{\phi_j\}_{j = 1}^m, \\ V_n = \text{span} \{\psi_i \}_{i = 1}^n. \end{align*} The weak problem reduces to problem of seeking some $u_m \in U_M$ for which the following $n+1$ equations hold: \begin{align*} \int_{\Omega} \psi_i' u_m' + \int_{\Omega} \psi_i u_m = \int_{\Omega} \psi_i f - \phi_i |_{\Gamma_N} g_N, \: i = 1, ... ,n \\ u_m|_{\Gamma_D} = g_D \text{ (enforce DOF bc directly without $V$) }. \end{align*} Expanding $u_m$ out in a basis, this becomes the linear system \begin{align*} \sum_{j = 1}^m \left( \int_{\Omega} \psi_i ' \phi_j ' + \psi_i \phi_j \right) u_j = \int_{\Omega} \psi_i f - \phi_i |_{\Gamma_N} g_N, \: i = 1, ..., n \\ \sum_{j = 1}^m u_j \phi_j = g_D. \end{align*} In matrix form (with m = 3, n = 2 for simplicity): \begin{align*} \begin{bmatrix} \int_{\Omega} (\psi_1' \phi_1' + \psi_1 \phi_1) & \int_{\Omega} (\psi_1' \phi_2' + \psi_1 \phi_2) & \int_{\Omega} (\psi_1' \phi_3' + \psi_1 \phi_3) \\ \int_{\Omega} (\psi_2' \phi_1' + \psi_2 \phi_1) & \int_{\Omega} (\psi_2' \phi_2' + \psi_2 \phi_2) & \int_{\Omega} (\psi_2' \phi_3' + \psi_2 \phi_3) \\ \phi_1 |_{\Gamma_D} & \phi_2 |_{\Gamma_D} & \phi_3 |_{\Gamma_D} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} \int_{\Omega} \psi_1 f - \phi_1|_{\Gamma_N} g_N \\ \int_{\Omega} \psi_2 f - \phi_2|_{\Gamma_N} g_N \\ g_D \end{bmatrix}. \end{align*} So, the last row to enforce essential boundary conditions extends the matrix to be square, as you mention. Once again, it was fine that we restricted our space $V$ to be zero on part of the boundary since the essential DOFs were enforced directly using the equation $\sum u_j \phi_j = g_D$, which does not involve any functions from $V$.

• Also note that the number of basis functions you use is often directly related to the number of DOFs (degrees of freedom). This is often equal to the number of nodes of your mesh if using 1st degree standard FEM spaces, but it can be greater than the number of nodes if using n-th degree standard FEM spaces, with n > 1. – amarney Aug 2 '18 at 16:42
• Homogenous conditions may be chosen for test functions to eliminate unwanted terms in the weak form. Also, a lot of the proofs in finite element theory are done for problems homogenous boundary conditions. I think not much generality is lost when doing so, as trace operators allow one to extend these results to problems with non-homgenous boundary conditions. – amarney Aug 2 '18 at 16:47
• You mean choose different $V_{0,h} \subset V$ and $V_{0,h} \subset V_h\subset V$, but the dimension of $V_{0,h}$ and $V_{h}$ are identical? I cannot imagine that. It may not be a B-Galerkin but a Petrov Galerkin. And I think linearity and sparsity of original formulation might be lost. – CatDog Aug 2 '18 at 17:01
• Sorry I cannot follow. Math is too heavy here when it comes to trace operator. I cannot even understand the definition of it. It seems like a boundary multi-value function extension but I do not know any concrete cases. – CatDog Aug 2 '18 at 17:10
• Actually, you are right in that the dimensions need not necessarily be the same. As you say, in Petrov-Galerkin you typically choose subspaces with different dimensions. I will edit my answer accordingly. – amarney Aug 2 '18 at 18:23