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!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qbM0Fge5$ ) 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!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qcjKXTR_$ >>> >>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qWpg0hOI$ >>> > >>> >>> >>> >> >> -- >> 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!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qcjKXTR_$ >> >> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qWpg0hOI$ >> > >> >> >> > > -- > 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!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qcjKXTR_$ > > <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qWpg0hOI$ > > > > > -- 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!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qcjKXTR_$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YK9T8_rEaNbJIndfx3tWHNlZlZ5GTj_KRDshR-nc5atZltZCcw73PtSNcrtvjWXZY73l5jQsAUx6qWpg0hOI$ >
