[petsc-users] Question regarding

[Matthew Knepley]

On Fri, Mar 20, 2026 at 3:09 PM Aldo Bonfiglioli <aldo.bonfiglioli at unibas.it>
wrote:

> Hi there,
>
> I am (successfully) using the following piece of code to insert elemental
> matrices into the "global" matrix "jac":
>
> !
>         if( matrix_assembly )then
> !
>            i = 0
>            PetscCall(DMPlexGetTransitiveClosure(dm, cell, useCone,
> PETSC_NULL_INTEGER, p_points, ierr))
> !
> !     loop over the vertices of the cell
> !
>            do j = vtxbgn, vtxend, 2  ! this is ndim-dependent
>               i = i + 1
>               point = p_points(j)
>               PetscCall(PetscSectionGetDof(section, point, ndofs, ierr))
>               PetscCall(PetscSectionGetOffset(section, point, offset,
> ierr))
>               row(i) = offset/ndofs
>               col(i) = offset/ndofs
>            end do
> !
>            do i = 1, nofvert
>               do j = 1, nofvert
>                  dummy = transpose(eltmat(:, :, i, j))
>                  ir = row(i)
>                  ic = col(j)
> !
> ! p_dummy is a 1d view of dummy
> !
>                  PetscCall(MatSetValuesBlockedLocal(jac, 2, [ir], 1,
> [ic], p_dummy, ADD_VALUES, ierr))
>               end do
>            end do
> !
>            PetscCall(DMPlexRestoreTransitiveClosure(dm, cell, useCone,
> PETSC_NULL_INTEGER, p_points, ierr))
> !
>         end if ! matrix_assembly
>
> I would however be happy to replace all of the above with a single call to:
>
> PetscCall(DMPlexMatSetClosure(dm, PETSC_NULL_SECTION, PETSC_NULL_SECTION,
> jac, cell, values, ADD_VALUES, ierr ))
>
> which I actually succeeded to do as long as ndofs = 1.
>
> The question is:
>
> 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 ?
>
> 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

  I = {v0,c0 v0,c1 ... v0,cN v1,c0 ,,, v1,cN ... vk,c0 ... vk,cN}

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.

  Does that make sense?

      Matt

> Thanks,
>
> Aldo
>
> --
> 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: https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!aau-IpT4Lq8cW7t_107COba4WHFpEfDnaq7pja4oDx-G4HY6csFSX8xjXHW4TLmAjFh6g9ucgfMj1MCZfTgg$  <https://urldefense.us/v3/__http://docenti.unibas.it/site/home/docente.html?m=002423__;!!G_uCfscf7eWS!ZusjOf_nCidORTOchtNsX_xGJZNjiAXMKKmn4rDhnR-3Y6mFzMF9Z7f7iJFeGj6tgcBzT2L5gtfxO56AMvPr3pgXHaCPziR-QL4$>
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aau-IpT4Lq8cW7t_107COba4WHFpEfDnaq7pja4oDx-G4HY6csFSX8xjXHW4TLmAjFh6g9ucgfMj1JExicVp$  <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aau-IpT4Lq8cW7t_107COba4WHFpEfDnaq7pja4oDx-G4HY6csFSX8xjXHW4TLmAjFh6g9ucgfMj1FEGLbyU$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20260320/b3fad6dd/attachment-0001.html>