[petsc-users] Dropping single entries in matrix insertion for block matrices

[Lawrence Mitchell]

Dear petsc-users,

I have a matrix with a block size > 1. I would sometimes like to insert into it, applying Dirichlet conditions to one component of the block. Now, for normal block (or when the block size is 1) dirichlet conditions, I do this by swapping out the local to global map for one with negative entries. That is, I do:

MatCreate(..., &mat);

MatSetSizes(mat, ...);

MatSetUp(mat);

MatSetLocalToGlobalMapping(mat, rmap, cmap);
// normal insertion
MatSetValuesLocal(mat, ...);
// or MatSetValuesBlockedLocal(mat, ...);

When I want to drop some entries I do:

ISLocalToGlobalMappingCreate(..., &rmap_with_bcs);

MatSetLocalToGlobalMapping(mat, rmap_with_bcs, cmap_with_bcs);

MatSetValuesLocal(mat, ...);

And those entries for which rmap_with_bcs maps to negative rows get dropped.

The problem arises when I want to do this for matrices with blocksize > 1.

Now the ISLocalToGlobalMapping must have the same block size as the matrix. But I can't then use it to drop single components.

In petsc4py speak I would like to be able to do:

from petsc4py import PETSc
mat = PETSc.Mat().create()
mat.setSizes((3, 3))
mat.setBlockSize(3)
mat.setUp()
lgmap = PETSc.LGMap().create([0, 1, -1], bsize=1)
mat.setLGMap(lgmap, lgmap)

This results in an error:

[0] MatSetLocalToGlobalMapping() line 2009 in matrix.c
[0] PetscLayoutSetISLocalToGlobalMapping() line 247 in pmap.c
[0] Petsc has generated inconsistent data
[0] Blocksize of layout 3 must match that of mapping 1

I suppose this is to catch me from shooting myself in the foot when subsequently doing MatSetValuesBlockedLocal? But I can arrange that this doesn't occur.

Thoughts?

If I could do this, it would make my code a lot simpler, but I am aware it might make error checking more difficult.

Lawrence