[petsc-users] Using MatSetValuesBlocked to build Finite Element Matrices with variable coefficients

[Nidish]

I'm trying to write a simple 1D finite element test case (linear 
elasticity). It was given in the manual that using an appropriate DM and 
calling MatSetValuesBlocked(). I'm however unable to do this for cases 
where I want to take runtime inputs for physical constants in the matrices.

One thing that puzzles me is the use of "const PetscInt*" and "const 
PetscScalar*" in the declaration of MatSetValuesBlocked(). I can't think 
of why this is done (if I am to loop through the elements and call this 
function, I don't see why the indices and/or values must be constant). 
The declaration I found was:

    PetscErrorCode  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode>  MatSetValues  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/MatSetValues.html#MatSetValues>(Mat  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/Mat.html#Mat>  mat,PetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  m,constPetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  idxm[],PetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  n,constPetscInt  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt>  idxn[],constPetscScalar  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscScalar.html#PetscScalar>  v[],InsertMode  <https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/InsertMode.html#InsertMode>  addv)

Here's the DM setup for my code:

    DM mshdm;
    Mater steel = {1.0, 2e11, 7800, 3.0, 400.0};// A structure with
    {Area, Modulus, Density, Length, Force} (runtime vars)
    KSP ksp;

    DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1, 1, NULL,
    &mshdm);// N is a runtime variable
    DMSetFromOptions(mshdm);
    DMSetUp(mshdm);
    DMSetApplicationContext(mshdm, &steel);  /* User Context */

    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetDM(ksp, mshdm);
    KSPSetComputeRHS(ksp, RHSFun, &steel);
    KSPSetComputeOperators(ksp, JACFun, &steel);
    KSPSetFromOptions(ksp);

I have no problems with the ComputeRHS function that I wrote, it seems 
to work fine. The ComputeOperators function on the other hand is giving 
me an "Argument out of range error" when MatSetValues() is called the 
second time (for e=1, discovered by debugging). Here's the function 
definition:

    PetscErrorCode JACFun(KSP ksp, Mat J, Mat jac, void* ctx)
    {
       Mater* pblm = (Mater*)ctx;
       PetscInt N, ix, n;
       PetscScalar Le;
       PetscErrorCode ierr;
       MatNullSpace nullspc;

       PetscScalar idx[2], idy[2], elstiff[2][2];// idx, idy are ids for
    MatSetValues()
       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);

       // Construct Element Stiffness
       Le = pblm->L/N;
       elstiff[0][0] = pblm->A*pblm->E/Le;
       elstiff[0][1] = -pblm->A*pblm->E/Le;
       elstiff[1][0] = -pblm->A*pblm->E/Le;
       elstiff[1][1] = pblm->A*pblm->E/Le;
       for (int e=0; e<N-1; ++e) {
         idx[0] = e; idx[1] = e+1;
         idy[0] = e; idy[1] = e+1;

         // This is the version I want to get working, but throws an
    argument out of range error
         MatSetValuesBlocked(jac, 2, (const PetscInt*)idx, 2, (const
    PetscInt*)idy,
                     (const PetscScalar*)elstiff, ADD_VALUES);

         // This version seemed to not construct the correct matrix
    (overlapped elements were added)
         /* MatSetValue(jac, e, e, elstiff[0][0], ADD_VALUES); */
         /* MatSetValue(jac, e, e+1, elstiff[0][1], ADD_VALUES); */
         /* MatSetValue(jac, e+1, e+1, elstiff[1][0], ADD_VALUES); */
         /* MatSetValue(jac, e+1, e+1, elstiff[1][1], ADD_VALUES); */
       }
       MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
       ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

       // I'm doing a "free-free" simulation, so constraining out the
    nullspace
    MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &nullspc);
       MatSetNullSpace(J, nullspc);
       MatNullSpaceDestroy(&nullspc);
       PetscFunctionReturn(0);
    }

I'd be very grateful if someone could point out what's the issue here.

Thank you,
Nidish
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200809/33c0c8ef/attachment.html>