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

  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/

--

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/