<div dir="ltr"><div dir="ltr">On Fri, Mar 20, 2026 at 3:09 PM Aldo Bonfiglioli <<a href="mailto:aldo.bonfiglioli@unibas.it">aldo.bonfiglioli@unibas.it</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><u></u>
<div>
<p>Hi there,</p>
<p>I am (successfully) using the following piece of code to insert
elemental matrices into the "global" matrix "jac":</p>
<p>
</p><blockquote type="cite">
<span style="font-family:monospace"><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">if(
matrix_assembly )then</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">i = 0</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">PetscCall(DMPlexGetTransitiveClosure(dm,
cell, useCone, PETSC_NULL_INTEGER, p_points, ierr))</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!
loop over the vertices of the cell</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">do j =
vtxbgn, vtxend, 2 ! this is ndim-dependent</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
i = i + 1</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
point = p_points(j)</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">PetscCall(PetscSectionGetDof(section,
point, ndofs, ierr))</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
PetscCall(PetscSectionGetOffset(section,
point, offset, ierr))</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
row(i) = offset/ndofs</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
col(i) = offset/ndofs</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
end do</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">do i = 1,
nofvert</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">do j = 1,
nofvert</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">dummy =
transpose(eltmat(:, :, i, j))</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">ir = row(i)</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">ic = col(j)</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">p_dummy</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">is</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">a</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">1d</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">view</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">of</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)"> </span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">dummy</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">PetscCall(MatSetValuesBlockedLocal(jac,
2, [ir], 1, [ic], p_dummy, ADD_VALUES, ierr))</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">end do</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">end do</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">PetscCall(DMPlexRestoreTransitiveClosure(dm,
cell, useCone, PETSC_NULL_INTEGER, p_points, ierr))</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">!</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">end if !
matrix_assembly</span><br>
<span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">
</span><br>
</span>
</blockquote>
I would however be happy to replace all of the above with a single
call to:<p></p>
<p>
</p><blockquote type="cite">
<span style="font-family:monospace"><span style="color:rgb(0,0,0);background-color:rgb(255,255,255)">PetscCall(DMPlexMatSetClosure(dm,
PETSC_NULL_SECTION, PETSC_NULL_SECTION, jac, cell, values,
ADD_VALUES, ierr ))</span><br>
</span>
</blockquote>
which I actually succeeded to do as long as ndofs = 1.<p></p>
<p>The question is:</p>
<p>how (in which order) should I "transfer" the entries of
eltmat(1:ndofs,1:ndofs,1:nofvert,1:nofvert) into the values array
required by DMPlexMatSetClosure ?</p>
<p></p></div></blockquote><div>For a given field, the dofs on the closure are ordered first by mesh point, and then by field component. So all field components are contiguous on a mesh point. So in your example, the dof order would be</div><div><br></div><div> I = {v0,c0 v0,c1 ... v0,cN v1,c0 ,,, v1,cN ... vk,c0 ... vk,cN}</div><div><br></div><div>The matrix "values" is in row-major order, meaning that all column values for the first row are given, and then the next row begins.</div><div><br></div><div> Does that make sense?</div><div><br></div><div> Matt</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><p>Thanks,</p>
<p>Aldo</p>
<pre cols="72">--
Dr. Aldo Bonfiglioli
Associate professor of Fluid Mechanics
Dipartimento di Ingegneria
Universita' della Basilicata
V.le dell'Ateneo Lucano, 10 85100 Potenza ITALY
tel:+39.0971.205203 fax:+39.0971.205215
web: <a href="https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!ZusjOf_nCidORTOchtNsX_xGJZNjiAXMKKmn4rDhnR-3Y6mFzMF9Z7f7iJFeGj6tgcBzT2L5gtfxO56AMvPr3pgXHaCPziR-QL4$" target="_blank">http://docenti.unibas.it/site/home/docente.html?m=002423</a></pre>
</div>
</blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aau-IpT4Lq8cW7t_107COba4WHFpEfDnaq7pja4oDx-G4HY6csFSX8xjXHW4TLmAjFh6g9ucgfMj1FEGLbyU$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>