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 <knep...@gmail.com> wrote: > On Sun, May 14, 2023 at 12:27 PM Zongze Yang <yangzon...@gmail.com> wrote: > >> >> >> >> On Sun, 14 May 2023 at 23:54, Matthew Knepley <knep...@gmail.com> wrote: >> >>> On Sun, May 14, 2023 at 9:21 AM Zongze Yang <yangzon...@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 <knep...@gmail.com> >>>> wrote: >>>> >>>>> On Sat, May 13, 2023 at 6:08 AM Zongze Yang <yangzon...@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 <yangzon...@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 <knep...@gmail.com> 于2022年6月18日周六 20:02写道: >>>>>>> >>>>>>>> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang <yangzon...@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 <yangzon...@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 <knep...@gmail.com> 于2022年6月17日周五 01:11写道: >>>>>>>>>> >>>>>>>>>>> On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang < >>>>>>>>>>> yangzon...@gmail.com> wrote: >>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> 在 2022年6月16日,23:22,Matthew Knepley <knep...@gmail.com> 写道: >>>>>>>>>>>> >>>>>>>>>>>> >>>>>>>>>>>> On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang < >>>>>>>>>>>> yangzon...@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/> >