[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Zongze Yang
yangzongze at gmail.com
Mon May 15 22:29:46 CDT 2023
Got it. Thank you for your explanation!
Best wishes,
Zongze
On Mon, 15 May 2023 at 23:28, Matthew Knepley <knepley at gmail.com> wrote:
> 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/20230516/b09a6a1b/attachment-0001.html>
More information about the petsc-users
mailing list