[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 15:24:31 CDT 2023
On Sun, May 14, 2023 at 12:27 PM Zongze Yang <yangzongze at gmail.com> wrote:
>
>
>
> On Sun, 14 May 2023 at 23:54, Matthew Knepley <knepley at gmail.com> wrote:
>
>> 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.
>>
>>
> I appreciate your offer to handle the MR. Please go ahead and take care of
> it. Thank you!
>
I have created the MR with your tests. They are working for me:
https://gitlab.com/petsc/petsc/-/merge_requests/6463
Thanks,
Matt
> Best Wishes,
> Zongze
>
>
>> 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/>
>>
>
--
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/6674c8c7/attachment-0001.html>
More information about the petsc-users
mailing list