On Friday, June 20th, 2025 at 9:37 PM, Matthew Knepley <knep...@gmail.com> wrote:
> 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: Thanks ,that will be handy. Perhaps this whole idea of using higher-order elements will offer no benefit for visualization and we;ll end up using linear elements only. On Friday, June 20th, 2025 at 9:37 PM, Matthew Knepley <knep...@gmail.com> wrote: > 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 The use of the flag "-..._petscspace_degree N" for higher order approximation space is indeed something we use. However, you mention this being a field, which is the part confusing me; I am not quite sure how to get coordinates values form it; I am only aware of functions such as where DMGetCoordinates / DMPlexGetCellCoordinates, which won't do the job here. The higher order space is indeed created, as shown with -dm_petscds_view --- Discrete System with 1 fields cell total dim 6 total comp 1 Field P2 FEM 1 component (implicit) (Nq 6 Nqc 1) 1-jet PetscFE Object: P2 (cdm_) 1 MPI process type: basic Basic Finite Element in 2 dimensions with 1 components PetscSpace Object: P2 (cdm_) 1 MPI process type: poly Space in 2 variables with 1 components, size 6 Polynomial space of degree 2 PetscDualSpace Object: P2 (cdm_) 1 MPI process type: lagrange Dual space with 1 components, size 6 Continuous Lagrange dual space Quadrature on a triangle of order 4 on 6 points (dim 2) Weak Form System with 1 fields --- Thanks, Noam ------- On Friday, June 20th, 2025 at 9:37 PM, Matthew Knepley <knep...@gmail.com> wrote: > 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!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPshVSYYl$ > > > 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!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPhIhe3Gt$ >>> ) 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/*(http:/*www.cse.buffalo.edu/*knepley/)__;fl0vfg!!G_uCfscf7eWS!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPnN3rfK9$ >>>>>>>>> >>>>>>> >>>>>>> -- >>>>>>> >>>>>>> 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!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPnN3rfK9$ >>>>>>> >>>>> >>>>> -- >>>>> >>>>> 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!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPnN3rfK9$ >>>>> >>> >>> -- >>> >>> 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!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPnN3rfK9$ >>> > > -- > > 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!f8aNNKGwSGN8vZeoKzWun79cKAhlzOkGLwSR4km2j_SaiJJhvUBmTsk094zqOSUvigxdI4J6WmRkIJN617sEAidRPnN3rfK9$ >
