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