On Fri, Jun 20, 2025 at 5:14 PM Noam T. <dontbugthed...@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 < > knep...@gmail.com> wrote: > > On Wed, Jun 18, 2025 at 6:49 PM Noam T. <dontbugthed...@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 < >> 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/__;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$ >
