Hi, I want to replace the “dolfin” with the latest “dolfinx” in my old project, but I met a problem: the result from “fenics.assemble_system” is different compared with “dolfinx.fem.assemble_matrix” in dolfinx. I find that the result from “dolfinx.fem.assemble_matrix” is equal to the result of “fenics.assemble_system” without boundary conditions. So how can I apply the boundary conditions to the left-hand side matrix?

I am wondering what is difference between A.zeroRowsLocal and cpp.fem.petsc.insert_diagonal? It seem to result in a more correct result with A.zeroRowsLocal.

I see you do the following to manually apply boundary conditions for the LHS matrix A, where zero the selected row in local order of matrix and leave diagnoal 1:

for bc in bcs:
dofs, _ = bc._cpp_object.dof_indices()
A.zeroRowsLocal(dofs, diag=1)

I was directed to see the following insert_diagonal when calling dolfinx.fem.petsc.assemble_matrix. I guess it is also used to apply bc with changing entities and diagonal:

if a.function_spaces[0] is a.function_spaces[1]:
A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
_cpp.fem.petsc.insert_diagonal(A, a.function_spaces[0], _bcs, diagonal)

These operations do different things.
ZeroRowsLocal zeros a set of rows given by indices local to the process, and optionally sets a diagonal value.
Insert_diagonal changes the diagonal, But doesn’t Zero of diagonal entries.
Im not dure what you mean by one giving a more correct result than the other.

That does assembly with BCs as input, so rows and columns are already zeroed, as illustrated by the lifting/symmetric application of bcs

My problem is that I did not know how to manually apply bc for A correctly. In my applications, I need to first get the matrix A and then modify A with boundary conditions applied. So my code looks like

# 1. Get A and b
A = assemble_matrix(a, bcs)
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, self.bcs)
# 2. Code to modify A and b, e.g:
A += delta_A
b += delta_b
# 3. Apply bc for A and b again
_bcs = [bc._cpp_object for bc in bcs]
if a.function_spaces[0] is a.function_spaces[1]:
_cpp.fem.petsc.insert_diagonal(A, a.function_spaces[0], _bcs, diagonal=1.0)
A.assemble()
set_bc(b, bcs)

It does not give the desired values on the boundary, for instance, a non-zero velocity boundary value with no-slip bc. However, if I replace the above Part 3 by the following code, the result shows a possibly correct value on boundary:

# 3. Apply bc for A and b again
for bc in self.bcs:
dofs, _ = bc._cpp_object.dof_indices()
A.zeroRowsLocal(dofs, diag=1)
A.assemble()
set_bc(b, self.bcs)