[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Matthew Knepley
knepley at gmail.com
Sun May 14 10:54:37 CDT 2023
On Sun, May 14, 2023 at 9:21 AM Zongze Yang <yangzongze at gmail.com> wrote:
> 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
>
No problem. Glad it is working.
> 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.
>
Great! If you make an MR for this, you will be included on the next list of
PETSc contributors. Otherwise, I can do it.
Thanks,
Matt
> 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/>
>>
>
--
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/1bc63197/attachment-0001.html>
More information about the petsc-users
mailing list