DirichletBC::check_arguments is too strict for zeroing columns

[Miroslav Kuchta]

Hi, due to the check the following snippet fails, while it should work (and did until here). The condition is not appropriate for zeroing columns.

from dolfin import *

mesh = UnitSquareMesh(20, 20)

P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)

V = FunctionSpace(mesh, P2)
Q = FunctionSpace(mesh, P1)

u = TrialFunction(V)
q = TestFunction(Q)

a = inner(div(u), q)*dx
L = inner(Constant(0), q)*dx

A = assemble(a)
b = assemble(L)

bc = DirichletBC(V, Constant((0, 0)), 'on_boundary')
bc.zero_columns(A, b)

[Jan Blechta]

@mirok-w-simula so the checks now assume an application to row space which is not appropriate for zero_columns. Do you know how to modify the check appropriately?

[Jan Blechta]

Maybe there should just be A->size(1) instead of A->size(1) when operating on columns, right?

[Prof Garth Wells]

Zeroing columns is highly undesirable with CSR storage. Much better to not add non-zero entries in the first place. For this point of view, I would like to see this function removed (I don't know why it was added in the first place).

[Miro Kuchta]

@blechta correct, A->size(1) should be used for column operations. I don't have a good solution for how to introduce this test. @garth-wells cbc.block uses zero_columns to assemble symmetric block structured systems and the function was most likely added was cbc.block in mind @kentand.

[Jan Blechta]

Fixed in 608d275.

[Jan Blechta]

I'm also inclined to remove it. It has essentially no documentation - one needs to dig into the implementation to understand what it does -, no coverage (I have added two tests in 3e35408) and probably does not work in parallel, although there is no check for it.

[Prof Garth Wells]

Shall we remove it? There is nothing stopping packages that build on DOLFIN from implementing their own functions to modify matrices.

[Prof Garth Wells]

I should add that the DOLFIN symmetric assembler was generalised ~1 year ago to handle block systems and bcs.

[Miro Kuchta]

I am aware that Chris Richardson has made some work on block systems. @garth-wells is that branch what you are referring to?

[Prof Garth Wells]

The assembler code should be in master. @chris_richardson ?

[Jan Blechta]

Garth means rather that system assembler accepts any rectangular matrix, decides whether it is diagonal or off-diagonal block, and applies bcs accordingly. The problem with this is that the logic might not work everytime, for example

u, v = TrialFunction(P1), TestFunction(P1)
a00 = inner(grad(u), grad(v))*dx
a01 = inner(grad(u), grad(v))*dx
a10 = inner(grad(u), grad(v))*dx
a11 = inner(grad(u), grad(v))*dx

The assembler won't recognize that a01 and a10 are meant as off-diagonal blocks. The only reasonable design is that the assembler is informed about this by user (like zero_columns decides it by (undocumented!) diag_val == 0.0 criterion).

[Chris Richardson]

With the system assembler, we should be supplying a set of matrices [a00, a01, a10, a11] , and appropriate BCs [bc0, bc1] to assemble together. It is still in the chris/petsc-matnest branch. As @garth says, this should be merged.

[kentand]

So, any news on chris/pets-matnest? Symmetric block matrices are convenient and I hate to see this functionality removed just because it is slightly impractical.

[Jan Blechta]

@kentand, as of I know it's not removed yet. We only proposed that for removing and it's for the discussion.

The slight impracticality would not matter. The real problem is that there is

• missing careful documentation (to demonstrate it, the function magically figures out whether the block is diagonal or off-diagonal depending on passed value of diag_val which is nowhere documented)
• missing any coverage (well I added a bit recently)
• undefined parallel behaviour (raise error or make it working)

If the issues were resolved then there could be a possibility to keep the function, I think.