DirichletBC::check_arguments is too strict for zeroing columns
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)
Comments (15)
-
-
Maybe there should just be
A->size(1)
instead ofA->size(1)
when operating on columns, right? -
- changed milestone to 2017.1
- changed component to dofmaps
-
assigned issue to
- marked as minor
-
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).
-
@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 useszero_columns
to assemble symmetric block structured systems and the function was most likely added was cbc.block in mind @kentand. -
- changed status to resolved
Fixed in 608d275.
-
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.
-
Shall we remove it? There is nothing stopping packages that build on DOLFIN from implementing their own functions to modify matrices.
-
I should add that the DOLFIN symmetric assembler was generalised ~1 year ago to handle block systems and bcs.
-
I am aware that Chris Richardson has made some work on block systems. @garth-wells is that branch what you are referring to?
-
The assembler code should be in master. @chris_richardson ?
-
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
anda10
are meant as off-diagonal blocks. The only reasonable design is that the assembler is informed about this by user (likezero_columns
decides it by (undocumented!)diag_val == 0.0
criterion). -
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. -
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.
-
@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.
- missing careful documentation (to demonstrate it, the function magically figures out whether the block is diagonal or off-diagonal depending on passed value of
- Log in to comment
@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?