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 <knepley@gmail.com> 于2022年6月18日周六 20:02写道：
On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang <yangzongze@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 <yangzongze@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 <knepley@gmail.com> 于2022年6月17日周五 01:11写道：
On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang <yangzongze@gmail.com> wrote:


在 2022年6月16日，23:22，Matthew Knepley <knepley@gmail.com> 写道：


On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang <yangzongze@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 seems return the coordinates in lex order.

Some code in reading gmsh file reads that


1756:     if (isSimplex) continuity = 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


  Thanks,

     Matt
 
Thanks,
Zongze
You can ask for the coordinates of a vertex or an edge directly using


by giving the vertex or edge point. You can get all the coordinates on a cell, in the closure order, using

  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



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



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