[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Zongze Yang
yangzongze at gmail.com
Sun May 14 10:09:24 CDT 2023
Yes, you are correct. I have conducted tests using high-order 3D meshes of
a unit cube, and regrettably, the tests have failed. I have attached the
files for your reference.
Kindly review the output provided below: ( The volume should be 1)
```
$ mpiexec -n 3 ./ex33 -coord_space 0 -dm_plex_filename cube_2rd.msh
-dm_plex_gmsh_project
-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true
-dm_plex_gmsh_project_fe_view
PetscFE Object: P2 3 MPI processes
type: basic
Basic Finite Element in 3 dimensions with 3 components
PetscSpace Object: P2 3 MPI processes
type: sum
Space in 3 variables with 3 components, size 30
Sum space of 3 concatenated subspaces (all identical)
PetscSpace Object: sum component (sumcomp_) 3 MPI processes
type: poly
Space in 3 variables with 1 components, size 10
Polynomial space of degree 2
PetscDualSpace Object: P2 3 MPI processes
type: lagrange
Dual space with 3 components, size 30
Continuous Lagrange dual space
Quadrature on a tetrahedron of order 5 on 27 points (dim 3)
Volume: 0.46875
$ mpiexec -n 3 ./ex33 -coord_space 0 -dm_plex_filename cube_3rd.msh
-dm_plex_gmsh_project
-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true
-dm_plex_gmsh_project_fe_view
PetscFE Object: P3 3 MPI processes
type: basic
Basic Finite Element in 3 dimensions with 3 components
PetscSpace Object: P3 3 MPI processes
type: sum
Space in 3 variables with 3 components, size 60
Sum space of 3 concatenated subspaces (all identical)
PetscSpace Object: sum component (sumcomp_) 3 MPI processes
type: poly
Space in 3 variables with 1 components, size 20
Polynomial space of degree 3
PetscDualSpace Object: P3 3 MPI processes
type: lagrange
Dual space with 3 components, size 60
Continuous Lagrange dual space
Quadrature on a tetrahedron of order 7 on 64 points (dim 3)
Volume: 0.536855
```
Best wishes,
Zongze
On Sun, 14 May 2023 at 22:36, Jed Brown <jed at jedbrown.org> wrote:
> Good to hear this works for you. I believe there is still a problem with
> high order tetrahedral elements (we've been coping with it for months and
> someone asked last week) and plan to look at it as soon as possible now
> that my semester finished.
>
> Zongze Yang <yangzongze at gmail.com> writes:
>
> > Hi, Matt,
> >
> > The issue has been resolved while testing on the latest version of PETSc.
> > It seems that the problem has been fixed in the following merge request:
> > https://gitlab.com/petsc/petsc/-/merge_requests/5970
> >
> > I sincerely apologize for any inconvenience caused by my previous
> message.
> > However, I would like to provide you with additional information
> regarding
> > the test files. Attached to this email, you will find two Gmsh files:
> > "square_2rd.msh" and "square_3rd.msh." These files contain high-order
> > triangulated mesh data for the unit square.
> >
> > ```
> > $ ./ex33 -coord_space 0 -dm_plex_filename square_2rd.msh
> > -dm_plex_gmsh_project
> > -dm_plex_gmsh_project_petscdualspace_lagrange_continuity true
> > -dm_plex_gmsh_project_fe_view -volume 1
> > PetscFE Object: P2 1 MPI process
> > type: basic
> > Basic Finite Element in 2 dimensions with 2 components
> > PetscSpace Object: P2 1 MPI process
> > type: sum
> > Space in 2 variables with 2 components, size 12
> > Sum space of 2 concatenated subspaces (all identical)
> > PetscSpace Object: sum component (sumcomp_) 1 MPI process
> > type: poly
> > Space in 2 variables with 1 components, size 6
> > Polynomial space of degree 2
> > PetscDualSpace Object: P2 1 MPI process
> > type: lagrange
> > Dual space with 2 components, size 12
> > Continuous Lagrange dual space
> > Quadrature on a triangle of order 5 on 9 points (dim 2)
> > Volume: 1.
> > $ ./ex33 -coord_space 0 -dm_plex_filename square_3rd.msh
> > -dm_plex_gmsh_project
> > -dm_plex_gmsh_project_petscdualspace_lagrange_continuity true
> > -dm_plex_gmsh_project_fe_view -volume 1
> > PetscFE Object: P3 1 MPI process
> > type: basic
> > Basic Finite Element in 2 dimensions with 2 components
> > PetscSpace Object: P3 1 MPI process
> > type: sum
> > Space in 2 variables with 2 components, size 20
> > Sum space of 2 concatenated subspaces (all identical)
> > PetscSpace Object: sum component (sumcomp_) 1 MPI process
> > type: poly
> > Space in 2 variables with 1 components, size 10
> > Polynomial space of degree 3
> > PetscDualSpace Object: P3 1 MPI process
> > type: lagrange
> > Dual space with 2 components, size 20
> > Continuous Lagrange dual space
> > Quadrature on a triangle of order 7 on 16 points (dim 2)
> > Volume: 1.
> > ```
> >
> > Thank you for your attention and understanding. I apologize once again
> for
> > my previous oversight.
> >
> > Best wishes,
> > Zongze
> >
> >
> > On Sun, 14 May 2023 at 16:44, Matthew Knepley <knepley at gmail.com> wrote:
> >
> >> On Sat, May 13, 2023 at 6:08 AM Zongze Yang <yangzongze at gmail.com>
> wrote:
> >>
> >>> Hi, Matt,
> >>>
> >>> There seem to be ongoing issues with projecting high-order coordinates
> >>> from a gmsh file to other spaces. I would like to inquire whether
> there are
> >>> any plans to resolve this problem.
> >>>
> >>> Thank you for your attention to this matter.
> >>>
> >>
> >> Yes, I will look at it. The important thing is to have a good test. Here
> >> are the higher order geometry tests
> >>
> >>
> >>
> https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/tests/ex33.c
> >>
> >> I take shapes with known volume, mesh them with higher order geometry,
> and
> >> look at the convergence to the true volume. Could you add a GMsh test,
> >> meaning the .msh file and known volume, and I will fix it?
> >>
> >> Thanks,
> >>
> >> Matt
> >>
> >>
> >>> Best wishes,
> >>> Zongze
> >>>
> >>>
> >>> On Sat, 18 Jun 2022 at 20:31, Zongze Yang <yangzongze at gmail.com>
> wrote:
> >>>
> >>>> Thank you for your reply. May I ask for some references on the order
> of
> >>>> the dofs on PETSc's FE Space (especially high order elements)?
> >>>>
> >>>> Thanks,
> >>>>
> >>>> Zongze
> >>>>
> >>>> Matthew Knepley <knepley at gmail.com> 于2022年6月18日周六 20:02写道：
> >>>>
> >>>>> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang <yangzongze at gmail.com>
> >>>>> wrote:
> >>>>>
> >>>>>> In order to check if I made mistakes in the python code, I try to
> use
> >>>>>> c code to show the issue on DMProjectCoordinates. The code and mesh
> file is
> >>>>>> attached.
> >>>>>> If the code is correct, there must be something wrong with
> >>>>>> `DMProjectCoordinates` or `DMPlexCreateGmshFromFile` for high-order
> mesh.
> >>>>>>
> >>>>>
> >>>>> Something is definitely wrong with high order, periodic simplices
> from
> >>>>> Gmsh. We had not tested that case. I am at a conference and cannot
> look at
> >>>>> it for a week.
> >>>>> My suspicion is that the space we make when reading in the Gmsh
> >>>>> coordinates does not match the values (wrong order).
> >>>>>
> >>>>> Thanks,
> >>>>>
> >>>>> Matt
> >>>>>
> >>>>>
> >>>>>> The command and the output are listed below: (Obviously the bounding
> >>>>>> box is changed.)
> >>>>>> ```
> >>>>>> $ ./test_gmsh_load_2rd -filename cube-p2.msh -old_fe_view
> -new_fe_view
> >>>>>> Old Bounding Box:
> >>>>>> 0: lo = 0. hi = 1.
> >>>>>> 1: lo = 0. hi = 1.
> >>>>>> 2: lo = 0. hi = 1.
> >>>>>> PetscFE Object: OldCoordinatesFE 1 MPI processes
> >>>>>> type: basic
> >>>>>> Basic Finite Element in 3 dimensions with 3 components
> >>>>>> PetscSpace Object: P2 1 MPI processes
> >>>>>> type: sum
> >>>>>> Space in 3 variables with 3 components, size 30
> >>>>>> Sum space of 3 concatenated subspaces (all identical)
> >>>>>> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
> >>>>>> type: poly
> >>>>>> Space in 3 variables with 1 components, size 10
> >>>>>> Polynomial space of degree 2
> >>>>>> PetscDualSpace Object: P2 1 MPI processes
> >>>>>> type: lagrange
> >>>>>> Dual space with 3 components, size 30
> >>>>>> Discontinuous Lagrange dual space
> >>>>>> Quadrature of order 5 on 27 points (dim 3)
> >>>>>> PetscFE Object: NewCoordinatesFE 1 MPI processes
> >>>>>> type: basic
> >>>>>> Basic Finite Element in 3 dimensions with 3 components
> >>>>>> PetscSpace Object: P2 1 MPI processes
> >>>>>> type: sum
> >>>>>> Space in 3 variables with 3 components, size 30
> >>>>>> Sum space of 3 concatenated subspaces (all identical)
> >>>>>> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
> >>>>>> type: poly
> >>>>>> Space in 3 variables with 1 components, size 10
> >>>>>> Polynomial space of degree 2
> >>>>>> PetscDualSpace Object: P2 1 MPI processes
> >>>>>> type: lagrange
> >>>>>> Dual space with 3 components, size 30
> >>>>>> Continuous Lagrange dual space
> >>>>>> Quadrature of order 5 on 27 points (dim 3)
> >>>>>> New Bounding Box:
> >>>>>> 0: lo = 2.5624e-17 hi = 8.
> >>>>>> 1: lo = -9.23372e-17 hi = 7.
> >>>>>> 2: lo = 2.72091e-17 hi = 8.5
> >>>>>> ```
> >>>>>>
> >>>>>> Thanks,
> >>>>>> Zongze
> >>>>>>
> >>>>>> Zongze Yang <yangzongze at gmail.com> 于2022年6月17日周五 14:54写道：
> >>>>>>
> >>>>>>> I tried the projection operation. However, it seems that the
> >>>>>>> projection gives the wrong solution. After projection, the
> bounding box is
> >>>>>>> changed! See logs below.
> >>>>>>>
> >>>>>>> First, I patch the petsc4py by adding `DMProjectCoordinates`:
> >>>>>>> ```
> >>>>>>> diff --git a/src/binding/petsc4py/src/PETSc/DM.pyx
> >>>>>>> b/src/binding/petsc4py/src/PETSc/DM.pyx
> >>>>>>> index d8a58d183a..dbcdb280f1 100644
> >>>>>>> --- a/src/binding/petsc4py/src/PETSc/DM.pyx
> >>>>>>> +++ b/src/binding/petsc4py/src/PETSc/DM.pyx
> >>>>>>> @@ -307,6 +307,12 @@ cdef class DM(Object):
> >>>>>>> PetscINCREF(c.obj)
> >>>>>>> return c
> >>>>>>>
> >>>>>>> + def projectCoordinates(self, FE fe=None):
> >>>>>>> + if fe is None:
> >>>>>>> + CHKERR( DMProjectCoordinates(self.dm, NULL) )
> >>>>>>> + else:
> >>>>>>> + CHKERR( DMProjectCoordinates(self.dm, fe.fe) )
> >>>>>>> +
> >>>>>>> def getBoundingBox(self):
> >>>>>>> cdef PetscInt i,dim=0
> >>>>>>> CHKERR( DMGetCoordinateDim(self.dm, &dim) )
> >>>>>>> diff --git a/src/binding/petsc4py/src/PETSc/petscdm.pxi
> >>>>>>> b/src/binding/petsc4py/src/PETSc/petscdm.pxi
> >>>>>>> index 514b6fa472..c778e39884 100644
> >>>>>>> --- a/src/binding/petsc4py/src/PETSc/petscdm.pxi
> >>>>>>> +++ b/src/binding/petsc4py/src/PETSc/petscdm.pxi
> >>>>>>> @@ -90,6 +90,7 @@ cdef extern from * nogil:
> >>>>>>> int DMGetCoordinateDim(PetscDM,PetscInt*)
> >>>>>>> int DMSetCoordinateDim(PetscDM,PetscInt)
> >>>>>>> int DMLocalizeCoordinates(PetscDM)
> >>>>>>> + int DMProjectCoordinates(PetscDM, PetscFE)
> >>>>>>>
> >>>>>>> int DMCreateInterpolation(PetscDM,PetscDM,PetscMat*,PetscVec*)
> >>>>>>> int DMCreateInjection(PetscDM,PetscDM,PetscMat*)
> >>>>>>> ```
> >>>>>>>
> >>>>>>> Then in python, I load a mesh and project the coordinates to P2:
> >>>>>>> ```
> >>>>>>> import firedrake as fd
> >>>>>>> from firedrake.petsc import PETSc
> >>>>>>>
> >>>>>>> # plex = fd.mesh._from_gmsh('test-fd-load-p2.msh')
> >>>>>>> plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh')
> >>>>>>> print('old bbox:', plex.getBoundingBox())
> >>>>>>>
> >>>>>>> dim = plex.getDimension()
> >>>>>>> # (dim, nc, isSimplex, k,
> >>>>>>> qorder, comm=None)
> >>>>>>> fe_new = PETSc.FE().createLagrange(dim, dim, True, 2,
> >>>>>>> PETSc.DETERMINE)
> >>>>>>> plex.projectCoordinates(fe_new)
> >>>>>>> fe_new.view()
> >>>>>>>
> >>>>>>> print('new bbox:', plex.getBoundingBox())
> >>>>>>> ```
> >>>>>>>
> >>>>>>> The output is (The bounding box is changed!)
> >>>>>>> ```
> >>>>>>>
> >>>>>>> old bbox: ((0.0, 1.0), (0.0, 1.0), (0.0, 1.0))
> >>>>>>> PetscFE Object: P2 1 MPI processes
> >>>>>>> type: basic
> >>>>>>> Basic Finite Element in 3 dimensions with 3 components
> >>>>>>> PetscSpace Object: P2 1 MPI processes
> >>>>>>> type: sum
> >>>>>>> Space in 3 variables with 3 components, size 30
> >>>>>>> Sum space of 3 concatenated subspaces (all identical)
> >>>>>>> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
> >>>>>>> type: poly
> >>>>>>> Space in 3 variables with 1 components, size 10
> >>>>>>> Polynomial space of degree 2
> >>>>>>> PetscDualSpace Object: P2 1 MPI processes
> >>>>>>> type: lagrange
> >>>>>>> Dual space with 3 components, size 30
> >>>>>>> Continuous Lagrange dual space
> >>>>>>> Quadrature of order 5 on 27 points (dim 3)
> >>>>>>> new bbox: ((-6.530133708576188e-17, 36.30670832662781),
> (-3.899962995254311e-17, 36.2406171632539), (-8.8036464152166e-17,
> 36.111577025012224))
> >>>>>>>
> >>>>>>> ```
> >>>>>>>
> >>>>>>>
> >>>>>>> By the way, for the original DG coordinates, where can I find the
> relation of the closure and the order of the dofs for the cell?
> >>>>>>>
> >>>>>>>
> >>>>>>> Thanks!
> >>>>>>>
> >>>>>>>
> >>>>>>> Zongze
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>> Matthew Knepley <knepley at gmail.com> 于2022年6月17日周五 01:11写道：
> >>>>>>>
> >>>>>>>> On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang <
> yangzongze at gmail.com>
> >>>>>>>> wrote:
> >>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> 在 2022年6月16日，23:22，Matthew Knepley <knepley at gmail.com> 写道：
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang <
> yangzongze at gmail.com>
> >>>>>>>>> wrote:
> >>>>>>>>>
> >>>>>>>>>> Hi, if I load a `gmsh` file with second-order elements, the
> >>>>>>>>>> coordinates will be stored in a DG-P2 space. After obtaining the
> >>>>>>>>>> coordinates of a cell, how can I map the coordinates to vertex
> and edge?
> >>>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> By default, they are stored as P2, not DG.
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> I checked the coordinates vector, and found the dogs only defined
> >>>>>>>>> on cell other than vertex and edge, so I said they are stored as
> DG.
> >>>>>>>>> Then the function DMPlexVecGetClosure
> >>>>>>>>> <
> https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexVecGetClosure/>
> seems return
> >>>>>>>>> the coordinates in lex order.
> >>>>>>>>>
> >>>>>>>>> Some code in reading gmsh file reads that
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> 1756: if (isSimplex) continuity = PETSC_FALSE
> >>>>>>>>> <https://petsc.org/main/docs/manualpages/Sys/PETSC_FALSE/>; /*
> XXX
> >>>>>>>>> FIXME Requires DMPlexSetClosurePermutationLexicographic() */
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> 1758: GmshCreateFE(comm, NULL, isSimplex, continuity,
> >>>>>>>>> nodeType, dim, coordDim, order, &fe)
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> The continuity is set to false for simplex.
> >>>>>>>>>
> >>>>>>>>
> >>>>>>>> Oh, yes. That needs to be fixed. For now, you can just project it
> to
> >>>>>>>> P2 if you want using
> >>>>>>>>
> >>>>>>>>
> https://petsc.org/main/docs/manualpages/DM/DMProjectCoordinates/
> >>>>>>>>
> >>>>>>>> Thanks,
> >>>>>>>>
> >>>>>>>> Matt
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>> Thanks,
> >>>>>>>>> Zongze
> >>>>>>>>>
> >>>>>>>>> You can ask for the coordinates of a vertex or an edge directly
> >>>>>>>>> using
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexPointLocalRead/
> >>>>>>>>>
> >>>>>>>>> by giving the vertex or edge point. You can get all the
> coordinates
> >>>>>>>>> on a cell, in the closure order, using
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexVecGetClosure/
> >>>>>>>>>
> >>>>>>>>> Thanks,
> >>>>>>>>>
> >>>>>>>>> Matt
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>> Below is some code load the gmsh file, I want to know the
> relation
> >>>>>>>>>> between `cl` and `cell_coords`.
> >>>>>>>>>>
> >>>>>>>>>> ```
> >>>>>>>>>> import firedrake as fd
> >>>>>>>>>> import numpy as np
> >>>>>>>>>>
> >>>>>>>>>> # Load gmsh file (2rd)
> >>>>>>>>>> plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh')
> >>>>>>>>>>
> >>>>>>>>>> cs, ce = plex.getHeightStratum(0)
> >>>>>>>>>>
> >>>>>>>>>> cdm = plex.getCoordinateDM()
> >>>>>>>>>> csec = dm.getCoordinateSection()
> >>>>>>>>>> coords_gvec = dm.getCoordinates()
> >>>>>>>>>>
> >>>>>>>>>> for i in range(cs, ce):
> >>>>>>>>>> cell_coords = cdm.getVecClosure(csec, coords_gvec, i)
> >>>>>>>>>> print(f'coordinates for cell {i}
> :\n{cell_coords.reshape([-1,
> >>>>>>>>>> 3])}')
> >>>>>>>>>> cl = dm.getTransitiveClosure(i)
> >>>>>>>>>> print('closure:', cl)
> >>>>>>>>>> break
> >>>>>>>>>> ```
> >>>>>>>>>>
> >>>>>>>>>> Best wishes,
> >>>>>>>>>> Zongze
> >>>>>>>>>>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>> --
> >>>>>>>>> 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/>
> >>>>>>>>>
> >>>>>>>>>
> >>>>>>>>
> >>>>>>>> --
> >>>>>>>> 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/>
> >>>>>>>>
> >>>>>>>
> >>>>>
> >>>>> --
> >>>>> 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/>
> >>>>>
> >>>>
> >>
> >> --
> >> 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/>
> >>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230514/75008011/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cube_3rd.msh
Type: application/octet-stream
Size: 12223 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230514/75008011/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cube_2rd.msh
Type: application/octet-stream
Size: 4926 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230514/75008011/attachment-0003.obj>
More information about the petsc-users
mailing list