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