Ah I get it now, MatSetBlocked has to be set node-wise. I tried this and it works, thank you.

The other question I had was why are the arguments for MatSetValues() and MatSetValuesBlocked() set to const PetscInt* and const PetscScalar*  instead of just PetscInt* and PetscScalar* ? I have the typecast there so my flycheck doesn't keep throwing me warnings on emacs ;)

Thank You,
Nidish

On 8/10/20 5:16 PM, Jed Brown wrote:
Nidish <[email protected]> writes:

It's a 1D model with displacements and rotations as DoFs at each node.

I couldn't find much in the manual on MatSetBlockSize - could you
provide some more information on its use?

I thought since I've setup the system using DMDACreate1d (I've given
2dofs per node and a stencil width of 1 there), the matrix object should
have the nonzero elements preallocated. Here's the call to DMDACreate1d:

    DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 2, 1, NULL, &mshdm);
Ah, that will set the block size, but then it'll expect elstiff to be an 8x8 
matrix where you've only passed 4x4.

     idx[0] = 2*e; idx[1] = 2*e+1; idx[2] = 2*e+2; idx[3] = 2*e+3;

     MatSetValuesBlocked(jac, 4, (const PetscInt*)idx, 4, (const PetscInt*)idx,
                        (const PetscScalar*)elstiff, ADD_VALUES);

You don't need the casts in either case, BTW.  You probably want something like 
this.

     idx[0] = e; idx[1] = e + 1;

     MatSetValuesBlocked(jac, 2, idx, 2, idx, elstiff, ADD_VALUES);

Also, it might be more convenient to call MatSetValuesBlockedStencil(), 
especially if you move to a multi-dimensional problem at some point.
--
Nidish

Reply via email to