[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 01:16:08 CDT 2022
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.
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/>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220618/4a41551c/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cube-p2.msh
Type: application/octet-stream
Size: 5210 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220618/4a41551c/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_gmsh_load_2rd.c
Type: application/octet-stream
Size: 2297 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220618/4a41551c/attachment-0003.obj>
More information about the petsc-users
mailing list