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 <knep...@gmail.com>
wrote:
> On Tue, Jun 17, 2025 at 12:43 PM Noam T. <dontbugthed...@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 <knep...@gmail.com>
>> wrote:
>>
>>> On Thu, Jun 12, 2025 at 4:26 PM Noam T. <dontbugthed...@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 <knep...@gmail.com>
>>>> wrote:
>>>>
>>>>> On Thu, May 22, 2025 at 12:25 PM Noam T. <dontbugthed...@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$
>
connectivity.pdf
Description: Adobe PDF document
