I am trying to implement a FEM solver for the following initial value problem \begin{align} \frac{\partial u}{\partial t} - \nabla^2 u &= f\quad \text{ in } \Omega\times (0,T)\\ u &= g\quad \text{ on }\partial\Omega\times(0,T)\\ u(x,y,0) &= u_0(x,y) \end{align} where $\Omega\subset\mathbb{R}^2$ is just a rectangular domain. The functions $f$ and $g$ can vary with time. I have discretized this IVP into $$\mathbf{B}\frac{\partial}{\partial t}\mathbf{u} + \mathbf{A}\mathbf{u} = \mathbf{l} $$ where \begin{align} \text{B}_{ij} &= \int_\Omega \phi_i\phi_j\\ \text{A}_{ij} &= \int_\Omega \nabla\phi_i\nabla\phi_j\\ l_i &= \int_\Omega \phi_i f \end{align} and I can assemble these matrices and vector. I have also discretized time to get the time-stepping scheme $$(\mathbf{B}+k\theta\mathbf{A})\mathbf{u}^{n+1} = (\mathbf{B}-k(1-\theta)\mathbf{A})\mathbf{u}^n + k(\theta \mathbf{l}^n + (1-\theta)\mathbf{l}^{n+1}), $$ where $\theta\in[0,1]$ and $k$ is the time-step.

Here comes my problem. How exactly do I impose the boundary and initial condition, implementation-wise? What changes should I make to the matrices $\mathbf{A},\mathbf{B}$, etc? Should I do anything at each time-step? I feel like it shouldn't be too hard, however when the time is introduced it suddenly confuses me.



2 Answers 2


There are several equivalent ways of implementing Dirichlet boundary conditions with the finite element method. I'll give a brief overview, but you'll probably need more details, which can be found in the book Understanding and Implementing the Finite Element Method.

It may be easier to analyze this for the Poisson problem $-\nabla^2u = f$ than for the heat equation; the techniques are the same in both cases, and if you understand them for the steady-state problem then it shouldn't be a big leap to adapt them to the time-dependent one. For the initial value, that amounts to setting $u^0_i = \mathbf{u}_0(x_i)$. You'll then use this as part of the right-hand side in the linear system you've written down when you go to solve for the next time-step $\mathbf{u}^1$. The boundary conditions on $\partial\Omega$ are trickier.

One way to implement Dirichlet boundary conditions is to alter the matrices $A$ and $B$ and the load vector in order to guarantee that the boundary conditions are satisfied exactly. Suppose that $i$ is some unknown of the linear system corresponding to a Dirichlet boundary point; the modification you have to make is to zero out row $i$ of $B$, set row $i$ of $A$ to be the corresponding row in the identity matrix and set $\mathbf{l}_i$ equal to $g(x_i)$.

To illustrate this approach, suppose you had a $3 \times 3$ chunk of the linear system that looked like this:

$\left[\begin{matrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \\ u_3\end{matrix}\right] = \left[\begin{matrix} l_1 \\ l_2 \\ l_3\end{matrix}\right]$

and that $u_2$ is a degree of freedom corresponding to a boundary node. The modified system would look like this:

$\left[\begin{matrix} a_{11} & a_{12} & a_{13} \\ 0 & 1 & 0 \\ a_{31} & a_{32} & a_{33} \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \\ u_3\end{matrix}\right] = \left[\begin{matrix} l_1 \\ g_2 \\ l_3\end{matrix}\right]$

You can see how the modified system isn't symmetric anymore. Nonetheless, $u_1$ and $u_3$ are still coupled to $u_2$. The asymmetry of the system isn't usually an issue in practice provided that your initial guess for $u$ exactly satisfies the Dirichlet boundary conditions; try throwing the problem into a naive CG solver and see what happens. (Can you come up with a reason why this might be so?)

Another approach is to reduce the system size. Instead of solving for all of the coefficients $u_i$ in the expansion $\mathbf{u} = \sum_iu_i\phi_i$, you'll only solve for the coefficients corresponding to interior nodes. This means extracting sub-matrices out of $A$ and $B$ for those unknowns, and adding an extra component to the load vector $\mathbf{l}$ consisting of the effective forcing on the interior from boundary conditions. Referring back to the example matrices from before, we're left with the following reduced system for $u_1$ and $u_3$:

$\left[\begin{matrix} a_{11} & a_{13} \\ a_{31} & a_{33} \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_3\end{matrix}\right] = \left[\begin{matrix} l_1 - a_{12}g_2 \\ l_3 - a_{32}g_2\end{matrix}\right]$

Changing the linear system or deflating it are, from what I've observed, the most common approaches. You can also add on boundary conditions later as a constraint, or you can use penalty methods, but these are less common. Of the two, I tend to favor the first approach.

  • $\begingroup$ (+1) Thanks a lot for the answer. I will use the first type of approach, by altering the system. When I tried to do the changes myself, I encountered exactly the problem that the system no longer was postive definite, and I thought I did something wrong. I would love if you expanded on this. $\endgroup$
    – Eff
    Jan 22, 2015 at 22:50
  • $\begingroup$ Also, whenever I impose the BC for a certain boundary vertex, shouldn't I eliminate this unknown from the other rows of the system. Example, say $u_2$ is a boundary vertex. Something like this $$\begin{bmatrix}a_{1,1} & 0 & a_{1,3}\\ 0 & 1 & 0\\ a_{3,1} & 0 & a_{3,3} \end{bmatrix} \begin{bmatrix}u_1\\u_2\\u_3\end{bmatrix} = \begin{bmatrix}l_1-a_{1,2}g_2\\g_2\\l_3-a_{3,2}g_2\end{bmatrix}$$ $\endgroup$
    – Eff
    Jan 22, 2015 at 23:18
  • $\begingroup$ You can eliminate $u_2$ from the system that you've written down; you're sort of half-way between altering the matrix and reducing the system. I've edited my answer; let me know if it's clearer now. $\endgroup$ Jan 23, 2015 at 0:30
  • $\begingroup$ I will have to read it in detail tomorrow, it is getting late where I live. I guess om just still in a slight state of confusion on which matrices/vectors we do changes in when imposing boundary conditions. If you consider the first linear system in my question, I have seen in some literature that they also in each row subtract the factor $b_{i,2}\frac{\partial g}{\partial t}$ from the load vector, since my function $g$ depends on time. $\endgroup$
    – Eff
    Jan 23, 2015 at 0:51

I guess you are also missing a discretized weak form for the boundary condition $u = g$, i.e. $\mathbf{G} \mathbf{u} = \mathbf{m}$, where $\mathbf{G}$ is assembled from the boundary elements (the basis function set for elements of $\partial \Omega$). I found the treatment in Volakis, et al., "Finite Element Method for Electromagnetics" to be quite helpful regarding this basis set.

Anyway, if both $f$ and $g$ change in time, you need to recompute the respective right-hand sides before the time-stepping iteration.

Just a guess: what about unifying the volume balance equation and the boundary flux equation in one? As in

$\mathbf{B}_{\Omega} \, \partial_t \mathbf{u}_{\Omega} + \mathbf{A}_\Omega \mathbf{u}_\Omega + \mathbf{G}_{\partial \Omega} \mathbf{u}_{\partial \Omega} - \mathbf{l}_{\Omega} -\mathbf{m}_{\partial \Omega} = 0$,

where the $\partial \Omega$ and $\Omega$ subscripts should indicate the restrictions to the respective subdomains (the first method in Daniel's comment above)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.