[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 11:27:16 CDT 2023
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!
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/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230515/d079d9b6/attachment-0001.html>
More information about the petsc-users
mailing list