[petsc-users] Element connectivity of a DMPlex
Noam T.
dontbugthedevs at proton.me
Wed Jun 18 17:49:40 CDT 2025
See image attached.
Connectivity of the top mesh (first order triangle), can be obtained with the code shared before.
Connectivity of the bottom mesh (second order triangle) is what I would be interested in obtaining.
However, given your clarification on what the Plex and the PetscSection handle, it might not work; I am trying to get form the Plex what's only available from the PetscSection.
The purpose of this extended connectivity is plotting; in particular, using VTU files, where the "connectivity" of cells is required, and the extra nodes would be needed when using higher-order elements (e.g. VTK_QUADRATIC_TRIANGLE, VTK_QUADRATIC_QUAD, etc).
Perhaps I am over complicating things, and all this information can be obtained in a different, simpler way.
Thanks.
Noam
On Tuesday, June 17th, 2025 at 5:42 PM, Matthew Knepley <knepley at gmail.com> wrote:
> On Tue, Jun 17, 2025 at 12:43 PM Noam T. <dontbugthedevs at proton.me> wrote:
>
>> Thank you. For now, I am dealing with vertices only.
>>
>> Perhaps I did not explain myself properly, or I misunderstood your response.
>> What I meant to say is, given an element of order higher than one, the connectivity matrix I obtain this way only contains as many entries as the first order element: 3 for a triangle, 4 for a tetrahedron, etc.
>>
>> Looking at the closure of any cell in the mesh, this is also the case.However, the nodes are definitely present; e.g. from
>>
>> DMPlexGetCellCoordinates(dm, cell, NULL, nc, NULL, NULL)
>>
>> nc returns the expected value (12 for a 2nd order 6-node planar triangle, 30 for a 2nd order 10-node tetrahedron, etc).
>>
>> The question is, are the indices of these extra nodes obtainable in a similar way as with the code shared before? So that one can have e.g. [0, 1, 2, 3, 4, 5] for a second order triangle, not just [0, 1, 2].
>
> I am having a hard time understanding what you are after. I think this is because many FEM approaches confuse topology with analysis.
>
> The Plex stores topology, and you can retrieve adjacencies between any two mesh points.
>
> The PetscSection maps mesh points (cells, faces, edges , vertices) to sets of dofs. This is how higher order elements are implemented. Thus, we do not have to change topology to get different function spaces.
>
> The intended interface is for you to call DMPlexVecGetClosure() to get the closure of a cell (or face, or edge). You can also call DMPlexGetClosureIndices(), but index wrangling is what I intended to eliminate.
>
> What exactly are you looking for here?
>
> Thanks,
>
> Matt
>
>> Thank you.
>> Noam
>> On Friday, June 13th, 2025 at 3:05 PM, Matthew Knepley <knepley at gmail.com> wrote:
>>
>>> On Thu, Jun 12, 2025 at 4:26 PM Noam T. <dontbugthedevs at proton.me> wrote:
>>>
>>>> Thank you for the code; it provides exactly what I was looking for.
>>>>
>>>> Following up on this matter, does this method not work for higher order elements? For example, using an 8-node quadrilateral, exporting to a PETSC_VIEWER_HDF5_VIZ viewer provides the correct matrix of node coordinates in geometry/vertices
>>>
>>> If you wanted to include edges/faces, you could do it. First, you would need to decide how you would number things For example, would you number all points contiguously, or separately number cells, vertices, faces and edges. Second, you would check for faces/edges in the closure loop. Right now, we only check for vertices.
>>>
>>> I would say that this is what convinced me not to do FEM this way.
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>> (here a quadrilateral in [0, 10])
>>>> 5.0, 5.0
>>>> 0.0, 0.0
>>>> 10.0, 0.0
>>>> 10.0, 10.0
>>>> 0.0, 10.0
>>>> 5.0, 0.0
>>>> 10.0, 5.0
>>>> 5.0, 10.0
>>>> 0.0, 5.0
>>>>
>>>> but the connectivity in viz/topology is
>>>>
>>>> 0 1 2 3
>>>>
>>>> which are likely the corner nodes of the initial, first-order element, before adding extra nodes for the higher degree element.
>>>>
>>>> This connectivity values [0, 1, 2, 3, ...] are always the same, including for other elements, whereas the coordinates are correct
>>>>
>>>> E.g. for 3rd order triangle in [0, 1], coordinates are given left to right, bottom to top
>>>> 0, 0
>>>> 1/3, 0,
>>>> 2/3, 0,
>>>> 1, 0
>>>> 0, 1/3
>>>> 1/3, 1/3
>>>> 2/3, 1/3
>>>> 0, 2/3,
>>>> 1/3, 2/3
>>>> 0, 1
>>>>
>>>> but the connectivity (viz/topology/cells) is [0, 1, 2].
>>>>
>>>> Test meshes were created with gmsh from the python API, using
>>>> gmsh.option.setNumber("Mesh.ElementOrder", n), for n = 1, 2, 3, ...
>>>>
>>>> Thank you.
>>>> Noam
>>>> On Friday, May 23rd, 2025 at 12:56 AM, Matthew Knepley <knepley at gmail.com> wrote:
>>>>
>>>>> On Thu, May 22, 2025 at 12:25 PM Noam T. <dontbugthedevs at proton.me> wrote:
>>>>>
>>>>>> Hello,
>>>>>>
>>>>>> Thank you the various options.
>>>>>>
>>>>>> Use case here would be obtaining the exact output generated by option 1), DMView() with PETSC_VIEWER_HDF5_VIZ; in particular, the matrix generated under /viz/topology/cells.
>>>>>>
>>>>>>> There are several ways you might do this. It helps to know what you are aiming for.
>>>>>>>
>>>>>>> 1) If you just want this output, it might be easier to just DMView() with the PETSC_VIEWER_HDF5_VIZ format, since that just outputs the cell-vertex topology and coordinates
>>>>>>
>>>>>> Is it possible to get this information in memory, onto a Mat, Vec or some other Int array object directly? it would be handy to have it in order to manipulate it and/or save it to a different format/file. Saving to an HDF5 and loading it again seems redundant.
>>>>>>
>>>>>>> 2) You can call DMPlexUninterpolate() to produce a mesh with just cells and vertices, and output it in any format.
>>>>>>>
>>>>>>> 3) If you want it in memory, but still with global indices (I don't understand this use case), then you can use DMPlexCreatePointNumbering() for an overall global numbering, or DMPlexCreateCellNumbering() and DMPlexCreateVertexNumbering() for separate global numberings.
>>>>>>
>>>>>> Perhaps I missed it, but getting the connectivity matrix in /viz/topology/cells/ did not seem directly trivial to me from the list of global indices returned by DMPlexGetCell/Point/VertexNumbering() (i.e. I assume all the operations done when calling DMView()).
>>>>>
>>>>> Something like
>>>>>
>>>>> DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
>>>>> DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
>>>>> DMPlexGetVertexNumbering(dm, &globalVertexNumbers);
>>>>> ISGetIndices(globalVertexNumbers, &gv);
>>>>> for (PetscInt c = cStart; c < cEnd; ++c) {
>>>>> PetscInt *closure = NULL;
>>>>>
>>>>> DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);
>>>>> for (PetscInt cl = 0; c < Ncl * 2; cl += 2) {
>>>>> if (closure[cl] < vStart || closure[cl] >= vEnd) continue;
>>>>> const PetscInt v = gv[closure[cl]] < 0 ? -(gv[closure[cl]] + 1) : gv[closure[cl]];
>>>>>
>>>>> // Do something with v
>>>>> }
>>>>> DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);
>>>>> }
>>>>> ISRestoreIndices(globalVertexNumbers, &gv);
>>>>> ISDestroy(&globalVertexNumbers);
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>>> Thanks,
>>>>>> Noam.
>>>>>
>>>>> --
>>>>>
>>>>> 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/*(http:/*www.cse.buffalo.edu/*knepley/)__;fl0vfg!!G_uCfscf7eWS!fRqwbCoztVyyzAJCKwByl_ypV2r7m5XCOxm8fu0m5Tu0Fp8Ghj8I_m2t3XrOU9m7STyykqh6j29GaFAsb7TDQa2L3jVTMgl9$
>>>
>>> --
>>>
>>> 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/*(http:/*www.cse.buffalo.edu/*knepley/)__;fl0vfg!!G_uCfscf7eWS!fRqwbCoztVyyzAJCKwByl_ypV2r7m5XCOxm8fu0m5Tu0Fp8Ghj8I_m2t3XrOU9m7STyykqh6j29GaFAsb7TDQa2L3jVTMgl9$
>
> --
>
> 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/*(http:/*www.cse.buffalo.edu/*knepley/)__;fl0vfg!!G_uCfscf7eWS!fRqwbCoztVyyzAJCKwByl_ypV2r7m5XCOxm8fu0m5Tu0Fp8Ghj8I_m2t3XrOU9m7STyykqh6j29GaFAsb7TDQa2L3jVTMgl9$
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250618/a2921dba/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: connectivity.pdf
Type: application/pdf
Size: 12595 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250618/a2921dba/attachment-0001.pdf>
More information about the petsc-users
mailing list