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