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?
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?
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.