Good catch Stefano. Matt
On Thu, Oct 25, 2018 at 9:36 AM Stefano Zampini <stefano.zamp...@gmail.com> 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 < > stefano.zamp...@gmail.com> 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 < >> knep...@gmail.com> ha scritto: >> >>> On Wed, Oct 24, 2018 at 11:36 AM Maximilian Hartig < >>> imilian.har...@gmail.com> wrote: >>> >>>> >>>> >>>> On 24. Oct 2018, at 12:49, Matthew Knepley <knep...@gmail.com> wrote: >>>> >>>> On Wed, Oct 24, 2018 at 6:29 AM Lawrence Mitchell <we...@gmx.li> 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 < >>>>> imilian.har...@gmail.com> 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/>