[petsc-users] MatSetValues vs MatSetValuesBlocked

[Nidish]

Hello,

I've written a basic FE code of an Euler-Bernoulli Beam (4th order 
spatial deriv interpolated using Hermite elements).

While assembling matrices I noticed something peculiar - if I conducted 
assembly using MatSetValues it works, while if I did so using 
MatSetValuesBlocked it throws me an argument out of range error. You can 
observe this in the attached file in lines 94-95 where I have one 
version commented out and the other version active.

I.O.W., using
    MatSetValues(jac, 4, (const PetscInt*)idx, 4, (const 
PetscInt*)idx,(const PetscScalar*)elstiff, ADD_VALUES);
works, while
    MatSetValuesBlocked(jac, 4, (const PetscInt*)idx, 4, (const 
PetscInt*)idx,(const PetscScalar*)elstiff, ADD_VALUES);
does NOT work.

Also from the documentation I could not really understand what the 
difference was between these two.

Thank you,
Nidish
static char help[] = "1D Euler-Bernoulli Beam.\n";

#include <petscdm.h>
#include <petscdmda.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscmat.h>

PetscErrorCode RHSFun(KSP, Vec, void*);
PetscErrorCode StiffMat(KSP, Mat, Mat, void*);

typedef struct {
  PetscScalar A, E, I2, rho, L, F;
}Mater;

int main(int nargs, char *sargs[])
{
  PetscErrorCode ierr;
  PetscInt N=10;
  PetscMPIInt rank, size;
  ierr = PetscInitialize(&nargs, &sargs, (char*)0, help);
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  MPI_Comm_size(PETSC_COMM_WORLD, &size);

  PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);

  DM mshdm;
  Mater steelbm = {1.0, 2e11, 1.0/12.0, 7800, 3.0, 1e6};
  KSP ksp;
  
  DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 2, 1, NULL, &mshdm);
  DMSetFromOptions(mshdm);
  DMSetUp(mshdm);
  DMSetApplicationContext(mshdm, &steelbm);  /* User Context */

  KSPCreate(PETSC_COMM_WORLD, &ksp);
  KSPSetDM(ksp, mshdm);
  KSPSetComputeRHS(ksp, RHSFun, &steelbm);
  KSPSetComputeOperators(ksp, StiffMat, &steelbm);
  KSPSetFromOptions(ksp);
  
  KSPSolve(ksp, NULL, NULL);

  DMDestroy(&mshdm);
  KSPDestroy(&ksp);
  PetscFinalize();
  return 0;
}

PetscErrorCode RHSFun(KSP ksp, Vec b, void* ctx)
{
  Mater* steelbm = (Mater*)ctx;
  PetscInt N;
  PetscErrorCode ierr;  

  VecGetSize(b, &N);
  for (int i=0; i<N; ++i)
    VecSetValue(b, i, 0.0, INSERT_VALUES);
  VecSetValue(b, N-2, steelbm->F, INSERT_VALUES);
  VecAssemblyBegin(b);
  ierr = VecAssemblyEnd(b);
  
  PetscFunctionReturn(0);
}

PetscErrorCode StiffMat(KSP ksp, Mat J, Mat jac, void* ctx)
{
  Mater* steelbm = (Mater*)ctx;
  PetscInt N, ix, n;
  PetscScalar Le, EIbLec;
  PetscErrorCode ierr;
  DM mshdm;

  ierr = KSPGetDM(ksp, &mshdm);
  ierr = DMDAGetInfo(mshdm, NULL, &N, NULL, NULL, &n, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
  DMDAGetCorners(mshdm, &ix, NULL, NULL, &n, NULL, NULL);

  Le = steelbm->L/(N-1);
  EIbLec = steelbm->E*steelbm->I2/pow(Le, 3);
  PetscInt idx[4];
  PetscScalar elstiff[4][4] = {{EIbLec*12.0, EIbLec*6*Le, -EIbLec*12.0, EIbLec*6*Le},
  			       {EIbLec*6*Le, EIbLec*4*Le*Le, -EIbLec*6*Le, EIbLec*2*Le*Le},
  			       {-EIbLec*12.0, -EIbLec*6.0*Le, EIbLec*12.0, -EIbLec*6*Le},
  			       {EIbLec*6*Le, EIbLec*2*Le*Le, -EIbLec*6.0*Le, EIbLec*4*Le*Le}};

  PetscInt is, ie;
  is = ix;  ie = ix+n;
  PetscMPIInt rank;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  
  for (int e=is; e<((ie<N-1)?ie:(N-1)); ++e) {
    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); */
    MatSetValues(jac, 4, (const PetscInt*)idx, 4, (const PetscInt*)idx,
    		 (const PetscScalar*)elstiff, ADD_VALUES);
  }
  MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
  ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

  PetscInt rows[] = {0, 1};
  MatZeroRowsColumns(jac, 2, rows, EIbLec*12.0, NULL, NULL);

  PetscFunctionReturn(0);
}