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
(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!cZvIn6-WAARgI-p1R_s8_n7yI_bYSvNOSQA6c46k1dLOGF9AxVvPAYS3bARKuXY6jcvg_zko2D-CUgfpOxYtN35083cDfLZY$
>
