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