<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<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">PetscErrorCode</a> <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/MatSetValues.html#MatSetValues">MatSetValues</a>(<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Mat/Mat.html#Mat">Mat</a> mat,<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</a> m,const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</a> idxm[],<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</a> n,const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscInt.html#PetscInt">PetscInt</a> idxn[],const <a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/PetscScalar.html#PetscScalar">PetscScalar</a> v[],<a href="https://www.mcs.anl.gov/petsc/petsc-3.13/docs/manualpages/Sys/InsertMode.html#InsertMode">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;</tt><tt><br>
</tt><tt> for (int e=0; e<N-1; ++e) {</tt><tt><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+1, 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">I'd be very grateful if someone could
point out what's the issue here.<br>
</div>
<div class="moz-signature"><br>
</div>
<div class="moz-signature">Thank you,<br>
Nidish</div>
</body>
</html>