[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 08:20:58 CDT 2023
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/8e87f93e/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: square_3rd.msh
Type: application/octet-stream
Size: 1670 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230514/8e87f93e/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: square_2rd.msh
Type: application/octet-stream
Size: 1072 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230514/8e87f93e/attachment-0003.obj>
More information about the petsc-users
mailing list