Opened the PR https://bitbucket.org/petsc/petsc/pull-requests/1203/fix-dump-vtk-field-with-periodic-meshes/diff
Il giorno gio 25 ott 2018 alle ore 17:44 Matthew Knepley <[email protected]> ha scritto: > Good catch Stefano. > > Matt > > On Thu, Oct 25, 2018 at 9:36 AM Stefano Zampini <[email protected]> > wrote: > >> Maybe this is a fix >> >> diff --git a/src/dm/impls/plex/plexvtu.c b/src/dm/impls/plex/plexvtu.c >> index acdea12c2f..1a8bbada6a 100644 >> --- a/src/dm/impls/plex/plexvtu.c >> +++ b/src/dm/impls/plex/plexvtu.c >> @@ -465,10 +465,11 @@ PetscErrorCode DMPlexVTKWriteAll_VTU(DM >> dm,PetscViewer viewer) >> if ((closure[v] >= vStart) && (closure[v] < vEnd)) { >> PetscScalar *xpoint; >> >> - ierr = >> DMPlexPointLocalRead(dm,v,x,&xpoint);CHKERRQ(ierr); >> + ierr = >> DMPlexPointLocalRead(dm,closure[v],x,&xpoint);CHKERRQ(ierr); >> y[cnt + off++] = xpoint[i]; >> } >> } >> + cnt += off; >> ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, >> &closureSize, &closure);CHKERRQ(ierr); >> } >> } >> >> Max, does this fix your problem? If you confirm, I'll fix this in the >> maint branch >> >> If I run the below command line with the patch and with snes tutorials >> ex12 I get the nice picture attached >> $ ./ex12 -quiet -run_type test -interpolate 1 -bc_type dirichlet >> -petscspace_degree 1 -vec_view vtk:test.vtu:vtk_vtu -f >> ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh >> -dm_plex_gmsh_periodic -x_periodicity periodic -y_periodicity periodic >> -dm_refine 4 >> >> Il giorno gio 25 ott 2018 alle ore 15:11 Stefano Zampini < >> [email protected]> ha scritto: >> >>> Matt, >>> >>> you can reproduce it via >>> >>> $ valgrind ./ex12 -quiet -run_type test -interpolate 1 -bc_type >>> dirichlet -petscspace_degree 1 -vec_view vtk:test.vtu:vtk_vtu -f >>> ${PETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh >>> -dm_plex_gmsh_periodic >>> >>> Long time ago I added support for viewing meshes with periodic vertices >>> in the VTK_VTU viewer, but I did not fix the part that writes fields >>> >>> >>> Il giorno mer 24 ott 2018 alle ore 21:04 Matthew Knepley < >>> [email protected]> ha scritto: >>> >>>> On Wed, Oct 24, 2018 at 11:36 AM Maximilian Hartig < >>>> [email protected]> wrote: >>>> >>>>> >>>>> >>>>> On 24. Oct 2018, at 12:49, Matthew Knepley <[email protected]> wrote: >>>>> >>>>> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell <[email protected]> >>>>> wrote: >>>>> >>>>>> Hi Max, >>>>>> >>>>>> (I'm cc'ing in the petsc-users mailing list which may have more >>>>>> advice, if you are using PETSc you should definitely subscribe! >>>>>> >>>>>> > On 24 Oct 2018, at 09:27, Maximilian Hartig < >>>>>> [email protected]> wrote: >>>>>> > >>>>>> > Hello Lawrence, >>>>>> > >>>>>> > sorry to message you out of the blue. My name is Max and I found >>>>>> your post on GitHub ( >>>>>> https://github.com/firedrakeproject/firedrake/issues/1246 ) on >>>>>> DMPlex being able to read periodic gmsh files. I am currently trying to >>>>>> do >>>>>> just that (creating a periodic DMPlex mesh with gmsh) in the context of >>>>>> my >>>>>> PhD work. So far I haven’ t found any documentation on the periodic BC’s >>>>>> with DMPlex and gmsh in the official petsc documentation. >>>>>> > I was wondering whether you’d be so kind as to point me in a >>>>>> general direction concerning how to achieve this. You seem experienced in >>>>>> using petsc and I would greatly appreciate your help. >>>>>> >>>>>> >>>>>> I think the answer is "it depends". If you're just using DMPlex >>>>>> directly and all the of the functionality with PetscDS, then I /think/ >>>>>> that >>>>>> reading periodic meshes via gmsh (assuming you're using the appropriate >>>>>> gmsh mesh format [v2]) "just works". >>>>>> >>>>> >>>>> There are two phases here: topological and geometric. DMPlex >>>>> represents the periodic topological entity directly. For example, a >>>>> circle >>>>> is just a segment with one end hooked to the other. Vertices are not >>>>> duplicated, or mapped to each other. This makes topology simple and easy >>>>> to >>>>> implement. However, then geometry is more complicated. What Plex does is >>>>> allow coordinates to be represented by a discontinuous field taking values >>>>> on cells, in addition to vertices. In our circle example, each cells near >>>>> the cut will have 2 coordinates, one for each vertex, but they will not >>>>> agree across the cut. If you define a periodic domain, then Plex can >>>>> construct this coordinate field automatically using DMPlexLocalize(). >>>>> These >>>>> DG coordinates are then used by the integration routines. >>>>> >>>>> >>>>> Ok, I think I understand the concept. DMPlex reads information about >>>>> both topology and coordinates from the .msh file. Creating a periodic mesh >>>>> in gmsh then should allow DMPlex to identify the periodic boundaries as >>>>> the >>>>> “cut” and build the mesh topology accordingly. Coordinate information is >>>>> handled separately. >>>>> That means, as Lawrence suggested, building periodic meshes in gmsh >>>>> and reading them in to petsc’s DMPlex should indeed “just work”. (From >>>>> the >>>>> user perspective). The only extra step is to call DMLocalizeCoordinates() >>>>> after DMPlexReadFromFile(). Sorry to reiterate, I am just trying to make >>>>> sense of this. >>>>> >>>>> >>>>> >>>>>> From my side, the issue is to do with mapping that coordinate field >>>>>> into one that I understand (within Firedrake). You may not have this >>>>>> problem. >>>>>> >>>>> >>>>> Firedrake uses its own coordinate mapping and integration routines, so >>>>> they must manage the second part independently. I hope to get change this >>>>> slightly soon by making the Firedrake representation a DMField, so that it >>>>> looks the same to Plex. >>>>> >>>>> Thanks, >>>>> >>>>> Matt >>>>> >>>>> >>>>>> Thanks, >>>>>> >>>>>> Lawrence >>>>> >>>>> >>>>> >>>>> -- >>>>> 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://www.cse.buffalo.edu/~knepley/ >>>>> >>>>> >>>>> To read periodic meshes from GMSH, you need to use the option >>>>> -dm_plex_gmsh_periodic and DMPlexCreateFromFile >>>>> >>>>> >>>>> Ahh, thanks. I was missing the option " -dm_plex_gmsh_periodic “. But >>>>> using this option I now generate a segmentation fault error when calling >>>>> VecView() on the solution vector with vtk and hdf5 viewers. Any >>>>> suggestions? >>>>> >>>> >>>> Small example? VTK is deprecated. HDF5 should work, although it will >>>> require you to have proper coordinates I think. We have to >>>> think about what you mean. If its for a checkpoint, there is no >>>> problem, but for viz, those programs do not understand periodicity. Thus I >>>> embed it in a higher dimensional space. >>>> >>>> Matt >>>> >>>>> See src/dm/impls/plex/examples/tests/ex1.c. An example runs >>>>> >>>>> $ ./ex1 -filename >>>>> ${PETSC_DIR}/share/petsc/datafiles/meshes/cube_periodic_bin.msh >>>>> -dm_plex_gmsh_periodic -dm_view ::ascii_info_detail -interpolate >>>>> -test_shape >>>>> >>>>> generating periodic meshes in gmsh may be tricky, Lisandro for sure >>>>> may advice. >>>>> >>>>> >>>>> Thanks, >>>>> Max >>>>> >>>>> >>>> >>>> -- >>>> 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://www.cse.buffalo.edu/~knepley/ >>>> <http://www.cse.buffalo.edu/~knepley/> >>>> >>> >>> >>> -- >>> Stefano >>> >> >> >> -- >> Stefano >> > > > -- > 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://www.cse.buffalo.edu/~knepley/ > <http://www.cse.buffalo.edu/~knepley/> > -- Stefano
