[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Zongze Yang
yangzongze at gmail.com
Sat Jun 18 07:31:10 CDT 2022
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/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220618/0dc5c7a2/attachment.html>
More information about the petsc-users
mailing list