[petsc-users] Element connectivity of a DMPlex
Matthew Knepley
knepley at gmail.com
Fri Jun 20 16:37:24 CDT 2025
On Fri, Jun 20, 2025 at 5:14 PM Noam T. <dontbugthedevs at proton.me> wrote:
> Thank you once again, the code provides exactly what needed.
> An alternative for the VTK use was subdividing cells using corner nodes
> and integration points, such that all cells were first order. Any "better"
> alternative format/visualization software for this purpose?
>
I do have this implemented. If you give -dm_plex_high_order_view, it will
refine the grid and project into the linear space.
You can see that the code is pretty simple:
https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plex.c?ref_type=heads*L2021__;Iw!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-wmRROXs$
It uses DMPlexInterpolate() to connect the spaces.
> Possibly the last question in relation to this matter. We use first order
> meshes only, as you suggest, and let PETSc handle everything high-order
> through the approximation space. Hence, when retrieving node coordinates
> with DMGetCoordinates(Local) or DMPlexGetCellCoordinates, one gets the
> corner nodes only.
>
> Is the list of additional, high-order nodes coordinates readily available
> (stored) somewhere to be retrieved?
>
Yes. The idea is to think of coordinates as a discretized field on the
mesh, exactly as the solution field. Thus if you want higher order
coordinates, you choose a higher order coordinate space. I give the
coordinate space prefix cdm_, so you could say
-cdm_petscspace_degree 2
to get quadratic coordinates. There are some tests in Plex tests ex33.c
Thanks,
Matt
> They can be computed (e.g. using DMPlexReferenceToCoordinates knowing
> their position in the reference cell; or using corner nodes coordinates),
> but this will result in shared nodes being computed possibly several times;
> the large the mesh, the worse.
>
> E.g. in the example of the image attached before, DMPlexGetClosureIndices
> returns
> [4, 5, 6, 0, 1, 3] for the first cell
> [7, 8, 5, 1, 2, 3] for the second cell
> so that 3 nodes (5, 1, 3) will be computed twice if done naively, cell by
> cell.
>
> Thank you,
> Noam
> On Thursday, June 19th, 2025 at 12:43 AM, Matthew Knepley <
> knepley at gmail.com> wrote:
>
> On Wed, Jun 18, 2025 at 6:49 PM Noam T. <dontbugthedevs at proton.me> wrote:
>
>> 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).
>>
>
> Oh yes. VTK does this in a particularly ugly and backward way. Sigh. There
> is nothing we can do about this now, but someone should replace VTK with a
> proper interface at some point.
>
> So I understand why you want it and it is a defensible case, so here is
> how you get that (with some explanation). Those locations, I think, should
> not be understood as topological things, but rather as the locations of
> point evaluation functionals constituting a basis for the dual space (to
> your approximation space). I would call DMPlexGetClosureIndices() (
> https://urldefense.us/v3/__https://petsc.org/main/manualpages/DMPlex/DMPlexGetClosureIndices/__;!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0p5NBkz$ ) with
> a Section having the layout of P2 or Q2. This is the easy way to make that
>
> PetscSection gs;
> PetscFE fe;
> DMPolytopeType ct;
> PetscInt dim, cStart;
>
> PetscCall(DMGetDimension(dm, &dim));
> PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL));
> PetscCall(DMPlexGetCellType(dm, cStart, &ct));
> PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 2,
> PETSC_DETERMINE, &fe));
> PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
> PetscCall(PetscFEDestroy(&fe));
> PetscCall(DMCreateDS(dm));
> PetscCall(DMGetGlobalSection(dm, &gs));
>
> PetscInt *indices = NULL;
> PetscInt Nidx;
>
> PetscCall(DMPlexGetClosureIndices(dm, gs, gs, cell, PETSC_TRUE, &Nidx,
> &indices, NULL, NULL));
>
> Thanks,
>
> MAtt
>
>> 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/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY--vGRsSC$
>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$ >
>>>>
>>>>
>>>>
>>>
>>> --
>>> 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!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY--vGRsSC$
>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$ >
>>>
>>>
>>>
>>
>> --
>> 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!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY--vGRsSC$
>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$ >
>>
>>
>>
>
> --
> 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!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY--vGRsSC$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$ >
>
>
>
--
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!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY--vGRsSC$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Yrsqxk3Y_G1ua_QR9o4DIqJAMyqlrZjtniZoHJqKFNV8FYYjw0XmvDvvukBBCZwlSNnbLDWBwQlY-0bmZvjS$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250620/6e2c3d2e/attachment-0001.html>
More information about the petsc-users
mailing list