subject:"\[petsc\-users\] How to find the map between the high order coordinates of DMPlex and vertex numbering\?"

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-15 Thread Zongze Yang

Got it. Thank you for your explanation!

Best wishes,
Zongze


On Mon, 15 May 2023 at 23:28, Matthew Knepley  wrote:

> On Mon, May 15, 2023 at 9:55 AM Zongze Yang  wrote:
>
>> On Mon, 15 May 2023 at 17:24, Matthew Knepley  wrote:
>>
>>> On Sun, May 14, 2023 at 7:23 PM Zongze Yang 
>>> 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 +
>> #   [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-ma...@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 
 wrote:

> On Sun, May 14, 2023 at 12:27 PM Zongze Yang 
> wrote:
>
>>
>>
>>
>> On Sun, 14 May 2023 at 23:54, Matthew Knepley 
>> wrote:
>>
>>> On Sun, May 14, 2023 at 9:21 AM Zongze Yang 
>>> wrote:
>>>
 Hi, Matt,

 The issue has been resolved

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-15 Thread Matthew Knepley

On Mon, May 15, 2023 at 9:55 AM Zongze Yang  wrote:

> On Mon, 15 May 2023 at 17:24, Matthew Knepley  wrote:
>
>> On Sun, May 14, 2023 at 7:23 PM Zongze Yang  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 +
> #   [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-ma...@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  wrote:
>>>
 On Sun, May 14, 2023 at 12:27 PM Zongze Yang 
 wrote:

>
>
>
> On Sun, 14 May 2023 at 23:54, Matthew Knepley 
> wrote:
>
>> On Sun, May 14, 2023 at 9:21 AM Zongze Yang 
>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-15 Thread Zongze Yang

On Mon, 15 May 2023 at 17:24, Matthew Knepley  wrote:

> On Sun, May 14, 2023 at 7:23 PM Zongze Yang  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?

Additionally, should we always set
`dm_plex_gmsh_project_petscdualspace_lagrange_continuity` to false for
high-order gmsh files?

With the option set to `true`, I got the following error:
```
$ $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 +
#   [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-ma...@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  wrote:
>>
>>> On Sun, May 14, 2023 at 12:27 PM Zongze Yang 
>>> wrote:
>>>



 On Sun, 14 May 2023 at 23:54, Matthew Knepley 
 wrote:

> On Sun, May 14, 2023 at 9:21 AM Zongze Yang 
> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-15 Thread Matthew Knepley

On Sun, May 14, 2023 at 7:23 PM Zongze Yang  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.

  Thanks,

 Matt


> Best wishes,
> Zongze
>
>
> On Mon, 15 May 2023 at 04:24, Matthew Knepley  wrote:
>
>> On Sun, May 14, 2023 at 12:27 PM Zongze Yang 
>> wrote:
>>
>>>
>>>
>>>
>>> On Sun, 14 May 2023 at 23:54, Matthew Knepley  wrote:
>>>
 On Sun, May 14, 2023 at 9:21 AM Zongze Yang 
 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 
> wrote:
>
>> On Sat, May 13, 2023 at 6:08 AM Zongze Yang 
>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Zongze Yang

Could you try to project the coordinates into the continuity space by
enabling the option
`-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true`?

Best wishes,
Zongze


On Mon, 15 May 2023 at 04:24, Matthew Knepley  wrote:

> On Sun, May 14, 2023 at 12:27 PM Zongze Yang  wrote:
>
>>
>>
>>
>> On Sun, 14 May 2023 at 23:54, Matthew Knepley  wrote:
>>
>>> On Sun, May 14, 2023 at 9:21 AM Zongze Yang 
>>> 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 
 wrote:

> On Sat, May 13, 2023 at 6:08 AM Zongze Yang 
> 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 
>> 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  于2022年6月18日周六 20:02写道：
>>>
 On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
 wrote:

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Matthew Knepley

On Sun, May 14, 2023 at 12:27 PM Zongze Yang  wrote:

>
>
>
> On Sun, 14 May 2023 at 23:54, Matthew Knepley  wrote:
>
>> On Sun, May 14, 2023 at 9:21 AM Zongze Yang  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  wrote:
>>>
 On Sat, May 13, 2023 at 6:08 AM Zongze Yang 
 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 
> 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  于2022年6月18日周六 20:02写道：
>>
>>> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
>>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Zongze Yang

On Sun, 14 May 2023 at 23:54, Matthew Knepley  wrote:

> On Sun, May 14, 2023 at 9:21 AM Zongze Yang  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!

Best Wishes,
Zongze


>   Thanks,
>
>  Matt
>
>
>> Best wishes,
>> Zongze
>>
>>
>> On Sun, 14 May 2023 at 16:44, Matthew Knepley  wrote:
>>
>>> On Sat, May 13, 2023 at 6:08 AM Zongze Yang 
>>> 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  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  于2022年6月18日周六 20:02写道：
>
>> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Matthew Knepley

On Sun, May 14, 2023 at 9:21 AM Zongze Yang  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.

  Thanks,

 Matt


> Best wishes,
> Zongze
>
>
> On Sun, 14 May 2023 at 16:44, Matthew Knepley  wrote:
>
>> On Sat, May 13, 2023 at 6:08 AM Zongze Yang  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  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  于2022年6月18日周六 20:02写道：

> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Zongze Yang

Yes, you are correct. I have conducted tests using high-order 3D meshes of
a unit cube, and regrettably, the tests have failed. I have attached the
files for your reference.

Kindly review the output provided below: ( The volume should be 1)

```
$ mpiexec -n 3 ./ex33 -coord_space 0 -dm_plex_filename cube_2rd.msh
-dm_plex_gmsh_project
-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true
-dm_plex_gmsh_project_fe_view
PetscFE Object: P2 3 MPI processes
  type: basic
  Basic Finite Element in 3 dimensions with 3 components
  PetscSpace Object: P2 3 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_) 3 MPI processes
type: poly
Space in 3 variables with 1 components, size 10
Polynomial space of degree 2
  PetscDualSpace Object: P2 3 MPI processes
type: lagrange
Dual space with 3 components, size 30
Continuous Lagrange dual space
Quadrature on a tetrahedron of order 5 on 27 points (dim 3)
Volume: 0.46875
$ mpiexec -n 3 ./ex33 -coord_space 0 -dm_plex_filename cube_3rd.msh
-dm_plex_gmsh_project
-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true
-dm_plex_gmsh_project_fe_view
PetscFE Object: P3 3 MPI processes
  type: basic
  Basic Finite Element in 3 dimensions with 3 components
  PetscSpace Object: P3 3 MPI processes
type: sum
Space in 3 variables with 3 components, size 60
Sum space of 3 concatenated subspaces (all identical)
  PetscSpace Object: sum component (sumcomp_) 3 MPI processes
type: poly
Space in 3 variables with 1 components, size 20
Polynomial space of degree 3
  PetscDualSpace Object: P3 3 MPI processes
type: lagrange
Dual space with 3 components, size 60
Continuous Lagrange dual space
Quadrature on a tetrahedron of order 7 on 64 points (dim 3)
Volume: 0.536855
```

Best wishes,
Zongze


On Sun, 14 May 2023 at 22:36, Jed Brown  wrote:

> Good to hear this works for you. I believe there is still a problem with
> high order tetrahedral elements (we've been coping with it for months and
> someone asked last week) and plan to look at it as soon as possible now
> that my semester finished.
>
> Zongze Yang  writes:
>
> > 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
> >
> > 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.
> >
> > Best wishes,
> > Zongze
> >
> >
> > On Sun, 14 May 2023 at 16:44, Matthew Knepley  wrote:
> >
> >> On Sat, May 13, 2023 at 6:08

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Jed Brown

Good to hear this works for you. I believe there is still a problem with high 
order tetrahedral elements (we've been coping with it for months and someone 
asked last week) and plan to look at it as soon as possible now that my 
semester finished.

Zongze Yang  writes:

> 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
>
> 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.
>
> Best wishes,
> Zongze
>
>
> On Sun, 14 May 2023 at 16:44, Matthew Knepley  wrote:
>
>> On Sat, May 13, 2023 at 6:08 AM Zongze Yang  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  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  于2022年6月18日周六 20:02写道：

> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Zongze Yang

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

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.

Best wishes,
Zongze


On Sun, 14 May 2023 at 16:44, Matthew Knepley  wrote:

> On Sat, May 13, 2023 at 6:08 AM Zongze Yang  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  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  于2022年6月18日周六 20:02写道：
>>>
 On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-14 Thread Matthew Knepley

On Sat, May 13, 2023 at 6:08 AM Zongze Yang  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  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  于2022年6月18日周六 20:02写道：
>>
>>> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang 
>>> 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  于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, ) )
> 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 @@

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2023-05-13 Thread Zongze Yang

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.

Best wishes,
Zongze


On Sat, 18 Jun 2022 at 20:31, Zongze Yang  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  于2022年6月18日周六 20:02写道：
>
>> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang  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  于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, ) )
 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 =

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2022-06-18 Thread Zongze Yang

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  于2022年6月18日周六 20:02写道：

> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang  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  于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, ) )
>>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2022-06-18 Thread Matthew Knepley

On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang  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  于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, ) )
>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2022-06-17 Thread Zongze Yang

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, ) )
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  于2022年6月17日周五 01:11写道：

> On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang  wrote:
>
>>
>>
>> 在 2022年6月16日，23:22，Matthew Knepley  写道：
>>
>> 
>> On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang 
>> 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
>>  seems 
>> return
>> the coordinates in lex order.
>>
>> Some code in reading gmsh file reads that
>>
>>
>> 1756: if (isSimplex) continuity = PETSC_FALSE
>> ; /* XXX FIXME
>> Requires DMPlexSetClosurePermutationLexicographic() */
>>
>>
>> 1758: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim,
>> coordDim, order, )
>>
>>
>> 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

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2022-06-16 Thread Matthew Knepley

On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang  wrote:

>
>
> 在 2022年6月16日，23:22，Matthew Knepley  写道：
>
> 
> On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang  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
>  seems 
> return
> the coordinates in lex order.
>
> Some code in reading gmsh file reads that
>
>
> 1756: if (isSimplex) continuity = PETSC_FALSE
> ; /* XXX FIXME
> Requires DMPlexSetClosurePermutationLexicographic() */
>
>
> 1758: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim,
> coordDim, order, )
>
>
> 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/
> 
>
>

-- 
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/

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2022-06-16 Thread Zongze Yang



> 在 2022年6月16日，23:22，Matthew Knepley  写道：
> 
> 
>> On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang  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 seems return the coordinates in lex order.

Some code in reading gmsh file reads that


1756: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires 
DMPlexSetClosurePermutationLexicographic() */

1758: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, 
coordDim, order, )

The continuity is set to false for simplex.

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/

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

2022-06-16 Thread Matthew Knepley

On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang  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.

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/

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?

19 matches

Site Navigation

Mail list logo

Footer information