Reapplying boundary conditions to the matrix

Hi all,

To solve some stationary problem, use dolfinx to assemble a matrix A which I modify slightly to obtain the matrix problem I need. However, now I need to apply the boundary conditions (again), to ensure my problem still satisfies the boundary conditions. However, I cannot find how to do this, can anybody give me a hint where to look?


from dolfinx import mesh
from dolfinx.fem import FunctionSpace, locate_dofs_topological, dirichletbc, form
from dolfinx.fem.function import Function
from dolfinx.fem.petsc import assemble_matrix, set_bc
from mpi4py import MPI
from ufl import grad, dx, inner, TrialFunction, TestFunction

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("Lagrange", 1))

u = TrialFunction(V)
v = TestFunction(V)

# Boundary stuff.
tdim = domain.topology.dim
fdim = tdim - 1
uD = Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = locate_dofs_topological(V, 1, boundary_facets)
bc = dirichletbc(uD, boundary_dofs)

A = assemble_matrix(form(inner(grad(u), grad(v)) * dx))

A.setValue(5, 5, 49)  # Some modification to the A matrix

set_bc(A, [bc])  # Such function does not exist for the matrix, only for the vector?

Many thanks in advance,


Dolfinx uses lifting. This means that in the matrix every row and colum with a bc dof is set to the identity row, and then the corresponding Ag is subtracted from the RHs. See Why the coupling blocks in the stiffness matrix are zero? - #4 by dokken
How to set BCs for PETSc Matrices and PETSc Vectors in dolfinx - #2 by dokken
dolfinx.fem ā€” DOLFINx 0.6.0 documentation


Dear Dokken, thank you for your answer.

If I understand this correctly we put the rows and columns of the boundary dofs to zero, with the identity at the diagonal (in these rows/columns). Next, we multiply this with a vector \vec{g} and substract the result from the rhs vector b. However, If I look at a toy example where the indices of the boundary dofs are 1,2,3, I arrive at the following result for Ag:

Ag =\begin{bmatrix}1&0&0&\\0&1&0&\cdots\\0&0&1&\\& \vdots&&B \end{bmatrix} \cdot \begin{bmatrix}g_1\\g_2\\g_3\\0\\\vdots\\0\end{bmatrix} = \vec{g}

It seems to me that this is not archiving anything, as the result of this matrix-vector product just equals \vec{g}. Hence, I think Iā€™m missing something, can you give me a hint what Iā€™m missing?

  • Wouter

This is not what we do.
Imagine you have a system A x = b where no boundary conditions have been applied yet.
We take the matrix A, and multiply it by the g you specified, creating a vector Ag
We take the rhs of the equation, and subtract this vector:
On the lhs we create another matrix A_sub described by the matrix A from your post.
For the dofs that are not constrained by Dirichlet boundary conditions We are now solving B x_nonbc = (L - Ag)_nonbc
We finally use the set_bc command to set rhs for the dofs with Dirihclet bc to the correct solution.

Se for instance 3.1.1 of: A review on some discrete variational techniques for the approximation of essential boundary conditions - Archive ouverte HAL

1 Like

Thank you very much for the answer and paper! This has resolved my questions thus far.