<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
  </head>
  <body>
    <p>Update: I had made some mistakes in the assembly in the previous
      version. They're fixed now (see below). <br>
    </p>
    <div class="moz-cite-prefix">On 8/9/20 9:13 PM, Nidish wrote:<br>
    </div>
    <blockquote type="cite"
      cite="mid:4960b5fb-d8fc-ecbf-4083-618a8a8d6642@rice.edu">
      <meta http-equiv="content-type" content="text/html; charset=UTF-8">
      <p>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.</p>
      <p>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:</p>
      <blockquote>
        <pre><a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscErrorCode.html#PetscErrorCode" moz-do-not-send="true">PetscErrorCode</a> <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/MatSetValues.html#MatSetValues" moz-do-not-send="true">MatSetValues</a>(<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/Mat.html#Mat" moz-do-not-send="true">Mat</a> mat,<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> m,const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> idxm[],<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> n,const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt" moz-do-not-send="true">PetscInt</a> idxn[],const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscScalar.html#PetscScalar" moz-do-not-send="true">PetscScalar</a> v[],<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/InsertMode.html#InsertMode" moz-do-not-send="true">InsertMode</a> addv)
</pre>
      </blockquote>
      <div class="moz-signature">Here's the DM setup for my code: <tt><br>
        </tt></div>
      <blockquote>
        <div class="moz-signature"><tt>DM mshdm;</tt><tt><br>
          </tt><tt>Mater steel = {1.0, 2e11, 7800, 3.0, 400.0};</tt><tt>
            // A structure with {Area, Modulus, Density, Length, Force}
            (runtime vars)<br>
          </tt><tt>KSP ksp;</tt><tt><br>
          </tt><tt><br>
          </tt><tt>DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N,
            1, 1, NULL, &mshdm);</tt><tt> // N is a runtime variable<br>
          </tt><tt>DMSetFromOptions(mshdm);</tt><tt><br>
          </tt><tt>DMSetUp(mshdm);</tt><tt><br>
          </tt><tt>DMSetApplicationContext(mshdm, &steel);  /* User
            Context */</tt><tt><br>
          </tt><tt><br>
          </tt><tt>KSPCreate(PETSC_COMM_WORLD, &ksp);</tt><tt><br>
          </tt><tt>KSPSetDM(ksp, mshdm);</tt><tt><br>
          </tt><tt>KSPSetComputeRHS(ksp, RHSFun, &steel);</tt><tt><br>
          </tt><tt>KSPSetComputeOperators(ksp, JACFun, &steel);</tt><tt><br>
          </tt><tt>KSPSetFromOptions(ksp);</tt><br>
        </div>
      </blockquote>
      <div class="moz-signature">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:</div>
      <blockquote>
        <div class="moz-signature"><tt>PetscErrorCode JACFun(KSP ksp,
            Mat J, Mat jac, void* ctx)</tt><tt><br>
          </tt><tt>{</tt><tt><br>
          </tt><tt>  Mater* pblm = (Mater*)ctx;</tt><tt><br>
          </tt><tt>  PetscInt N, ix, n;</tt><tt><br>
          </tt><tt>  PetscScalar Le;</tt><tt><br>
          </tt><tt>  PetscErrorCode ierr;</tt><tt><br>
          </tt><tt>  MatNullSpace nullspc;</tt><tt><br>
          </tt><tt><br>
          </tt><tt>  PetscScalar idx[2], idy[2], elstiff[2][2];</tt><tt> 
            // idx, idy are ids for MatSetValues()<br>
          </tt><tt>  DM mshdm;</tt><tt><br>
          </tt><tt><br>
          </tt><tt>  ierr = KSPGetDM(ksp, &mshdm);</tt><tt><br>
          </tt><tt>  ierr = DMDAGetInfo(mshdm, NULL, &N, NULL, NULL,
            &n, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);</tt><tt><br>
          </tt><tt>  DMDAGetCorners(mshdm, &ix, NULL, NULL, &n,
            NULL, NULL);</tt><tt><br>
          </tt><tt><br>
          </tt><tt>  // Construct Element Stiffness</tt></div>
        <div class="moz-signature"><tt>  Le = pblm->L/N;</tt><tt><br>
          </tt><tt>  elstiff[0][0] = pblm->A*pblm->E/Le;</tt><tt><br>
          </tt><tt>  elstiff[0][1] = -pblm->A*pblm->E/Le;</tt><tt><br>
          </tt><tt>  elstiff[1][0] = -pblm->A*pblm->E/Le;  </tt><tt><br>
          </tt><tt>  elstiff[1][1] = pblm->A*pblm->E/Le;<br>
          </tt><br>
             <tt>PetscInt is, ie;</tt><tt><br>
          </tt><tt>
              MatGetOwnershipRange(jac, &is, &ie);</tt><tt><br>
          </tt><tt>  for (int e=is; e<((ie<N-1)?ie:(N-1)); ++e) {</tt><tt> 
            // Only looping over local elements<br>
          </tt><tt>    idx[0] = e; idx[1] = e+1;</tt><tt><br>
          </tt><tt>    idy[0] = e; idy[1] = e+1;    </tt><tt><br>
          </tt><tt><br>
          </tt><tt>    // This is the version I want to get working, but
            throws an argument out of range error<br>
          </tt></div>
        <div class="moz-signature"><tt>    MatSetValuesBlocked(jac, 2,
            (const PetscInt*)idx, 2, (const PetscInt*)idy,</tt><tt><br>
          </tt><tt>                (const PetscScalar*)elstiff,
            ADD_VALUES);</tt></div>
        <div class="moz-signature"><tt><br>
          </tt><tt>    // This version seemed to not construct the
            correct matrix (overlapped elements were added)<br>
          </tt><tt>    /* MatSetValue(jac, e, e, elstiff[0][0],
            ADD_VALUES); */</tt><tt><br>
          </tt><tt>    /* MatSetValue(jac, e, e+1, elstiff[0][1],
            ADD_VALUES); */</tt><tt><br>
          </tt><tt>    /* MatSetValue(jac, e+1, e, elstiff[1][0],
            ADD_VALUES); */</tt><tt><br>
          </tt><tt>    /* MatSetValue(jac, e+1, e+1, elstiff[1][1],
            ADD_VALUES); */</tt><tt><br>
          </tt><tt>  }</tt><tt><br>
          </tt><tt>  MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);</tt><tt><br>
          </tt><tt>  ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);</tt><tt><br>
          </tt><tt><br>
          </tt><tt>  // I'm doing a "free-free" simulation, so
            constraining out the nullspace</tt></div>
        <div class="moz-signature"><tt> 
            MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL,
            &nullspc);</tt><tt><br>
          </tt><tt>  MatSetNullSpace(J, nullspc);</tt><tt><br>
          </tt><tt>  MatNullSpaceDestroy(&nullspc);</tt><tt><br>
          </tt><tt>  PetscFunctionReturn(0);</tt><tt><br>
          </tt><tt>}</tt><br>
        </div>
      </blockquote>
      <div class="moz-signature"><br>
      </div>
      <div class="moz-signature">Thank you,<br>
        Nidish</div>
    </blockquote>
    <div class="moz-signature">-- <br>
      Nidish</div>
  </body>
</html>