[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 03:44:03 CDT 2023
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230514/485efecb/attachment-0001.html>
More information about the petsc-users
mailing list