[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Matthew Knepley
knepley at gmail.com
Sat Jun 18 07:02:44 CDT 2022
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/b6566f6c/attachment-0001.html>
More information about the petsc-users
mailing list