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?

MWE:

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,

Wouter

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
and
dolfinx.fem ā€” DOLFINx 0.6.0 documentation

2 Likes

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:
b-=Ag
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.