Lawrence,
I understand what you want and it is a reasonable request. The problem is
that currently ISLocalToGlobalMappingCreate() when used with block vectors and
matrices is always based on blocks, that is, from the manual page, "There is
one integer value in indices per block and it represents the actual indices
bs*idx + j, where j=0,..,bs-1". If you look at ISLocalToGlobalMappingApply() it
uses the code:
out[i] = idx[in[i]/bs]*bs + (in[i] % bs);
to do the mapping; that is the idx[] mapping is stored only by block, not by
point.
But I think you may be able to get the effect you want by "managing the local
to global mapping yourself before calling MatSetValues()". That is, you create
a ISLocalToGlobal object yourself, not based on blocks, then when setting your
Dirichlet conditions you call ISLocalToGlobalMapping() apply yourself and then
call MatSetValues() using the resulting indices.
Barry
> On Dec 9, 2018, at 9:55 AM, Lawrence Mitchell via petsc-users
> <petsc-users@mcs.anl.gov> wrote:
>
> 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
