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].
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!auwV8wTSwLcsJ8sDHyW_i-VTpb-2nu5SMxYOBxNVljOhAIQEeIy8yRyp5hfOho00flBooHSX92kCWgO6wrTmYwQKaJd1QT1x$
>>>
>
> --
>
> 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!auwV8wTSwLcsJ8sDHyW_i-VTpb-2nu5SMxYOBxNVljOhAIQEeIy8yRyp5hfOho00flBooHSX92kCWgO6wrTmYwQKaJd1QT1x$
>
