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