[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Matthew Knepley
knepley at gmail.com
Mon May 15 10:28:31 CDT 2023
On Mon, May 15, 2023 at 9:55 AM Zongze Yang <yangzongze at gmail.com> wrote:
> On Mon, 15 May 2023 at 17:24, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Sun, May 14, 2023 at 7:23 PM Zongze Yang <yangzongze at gmail.com> wrote:
>>
>>> Could you try to project the coordinates into the continuity space by
>>> enabling the option
>>> `-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true`?
>>>
>>
>> There is a comment in the code about that:
>>
>> /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
>>
>> So what is currently done is you project into the discontinuous space
>> from the GMsh coordinates,
>> and then we get the continuous coordinates from those later. This is why
>> we get the right answer.
>>
>>
> Sorry, I'm having difficulty understanding the comment and fully
> understanding your intended meaning. Are you saying that we can only
> project the space to a discontinuous space?
>
For higher order simplices, because we do not have the mapping to the GMsh
order yet.
> Additionally, should we always set
> `dm_plex_gmsh_project_petscdualspace_lagrange_continuity` to false for
> high-order gmsh files?
>
This is done automatically if you do not override it.
> With the option set to `true`, I got the following error:
>
Yes, do not do that.
Thanks,
Matt
> ```
> $ $PETSC_DIR/$PETSC_ARCH/tests/dm/impls/plex/tests/runex33_gmsh_3d_q2.sh
> -e "-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true"
> not ok dm_impls_plex_tests-ex33_gmsh_3d_q2 # Error code: 77
> # Volume: 0.46875
> # [0]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
> # [0]PETSC ERROR: Petsc has generated inconsistent data
> # [0]PETSC ERROR: Calculated volume 0.46875 != 1. actual volume
> (error 0.53125 > 1e-06 tol)
> # [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble
> shooting.
> # [0]PETSC ERROR: Petsc Development GIT revision:
> v3.19.1-294-g9cc24bc9b93 GIT Date: 2023-05-15 12:07:10 +0000
> # [0]PETSC ERROR: ../ex33 on a arch-linux-c-debug named AMA-PC-RA18
> by yzz Mon May 15 21:53:43 2023
> # [0]PETSC ERROR: Configure options
> --CFLAGS=-I/opt/intel/oneapi/mkl/latest/include
> --CXXFLAGS=-I/opt/intel/oneapi/mkl/latest/include
> --LDFLAGS="-Wl,-rpath,/opt/intel/oneapi/mkl/latest/lib/intel64
> -L/opt/intel/oneapi/mkl/latest/lib/intel64" --download-bison
> --download-chaco --download-cmake
> --download-eigen="/home/yzz/firedrake/complex-int32-mkl-X-debug/src/eigen-3.3.3.tgz
> " --download-fftw --download-hdf5 --download-hpddm --download-hwloc
> --download-libpng --download-metis --download-mmg --download-mpich
> --download-mumps --download-netcdf --download-p4est --download-parmmg
> --download-pastix --download-pnetcdf --download-ptscotch
> --download-scalapack --download-slepc --download-suitesparse
> --download-superlu_dist --download-tetgen --download-triangle
> --with-blaslapack-dir=/opt/intel/oneapi/mkl/latest --with-c2html=0
> --with-debugging=1 --with-fortran-bindings=0
> --with-mkl_cpardiso-dir=/opt/intel/oneapi/mkl/latest
> --with-mkl_pardiso-dir=/opt/intel/oneapi/mkl/latest
> --with-scalar-type=complex --with-shared-libraries=1 --with-x=1 --with-zlib
> PETSC_ARCH=arch-linux-c-debug
> # [0]PETSC ERROR: #1 CheckVolume() at
> /home/yzz/opt/petsc/src/dm/impls/plex/tests/ex33.c:246
> # [0]PETSC ERROR: #2 main() at
> /home/yzz/opt/petsc/src/dm/impls/plex/tests/ex33.c:261
> # [0]PETSC ERROR: PETSc Option Table entries:
> # [0]PETSC ERROR: -coord_space 0 (source: command line)
> # [0]PETSC ERROR: -dm_plex_filename
> /home/yzz/opt/petsc/share/petsc/datafiles/meshes/cube_q2.msh (source:
> command line)
> # [0]PETSC ERROR: -dm_plex_gmsh_project (source: command line)
> # [0]PETSC ERROR:
> -dm_plex_gmsh_project_petscdualspace_lagrange_continuity true (source:
> command line)
> # [0]PETSC ERROR: -tol 1e-6 (source: command line)
> # [0]PETSC ERROR: -volume 1.0 (source: command line)
> # [0]PETSC ERROR: ----------------End of Error Message -------send
> entire error message to petsc-maint at mcs.anl.gov----------
> # application called MPI_Abort(MPI_COMM_SELF, 77) - process 0
> ok dm_impls_plex_tests-ex33_gmsh_3d_q2 # SKIP Command failed so no diff
> ```
>
>
> Best wishes,
> Zongze
>
> Thanks,
>>
>> Matt
>>
>>
>>> Best wishes,
>>> Zongze
>>>
>>>
>>> On Mon, 15 May 2023 at 04:24, Matthew Knepley <knepley at gmail.com> wrote:
>>>
>>>> On Sun, May 14, 2023 at 12:27 PM Zongze Yang <yangzongze at gmail.com>
>>>> wrote:
>>>>
>>>>>
>>>>>
>>>>>
>>>>> 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!
>>>>>
>>>>
>>>> I have created the MR with your tests. They are working for me:
>>>>
>>>> https://gitlab.com/petsc/petsc/-/merge_requests/6463
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>>
>>>>> 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/>
>>>>>>
>>>>>
>>>>
>>>> --
>>>> 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/4df9824d/attachment-0001.html>
More information about the petsc-users
mailing list