<div dir="ltr">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.<div>If the code is correct, there must be something wrong with `DMProjectCoordinates` or `DMPlexCreateGmshFromFile` for high-order mesh.</div><div><br></div><div>The command and the output are listed below: (Obviously the bounding box is changed.)<div>```</div><div>$ ./test_gmsh_load_2rd -filename cube-p2.msh -old_fe_view -new_fe_view<br>Old Bounding Box:<br>  0: lo = 0. hi = 1.<br>  1: lo = 0. hi = 1.<br>  2: lo = 0. hi = 1.<br>PetscFE Object: OldCoordinatesFE 1 MPI processes<br>  type: basic<br>  Basic Finite Element in 3 dimensions with 3 components<br>  PetscSpace Object: P2 1 MPI processes<br>    type: sum<br>    Space in 3 variables with 3 components, size 30<br>    Sum space of 3 concatenated subspaces (all identical)<br>      PetscSpace Object: sum component (sumcomp_) 1 MPI processes<br>        type: poly<br>        Space in 3 variables with 1 components, size 10<br>        Polynomial space of degree 2<br>  PetscDualSpace Object: P2 1 MPI processes<br>    type: lagrange<br>    Dual space with 3 components, size 30<br>    Discontinuous Lagrange dual space<br>    Quadrature of order 5 on 27 points (dim 3)<br>PetscFE Object: NewCoordinatesFE 1 MPI processes<br>  type: basic<br>  Basic Finite Element in 3 dimensions with 3 components<br>  PetscSpace Object: P2 1 MPI processes<br>    type: sum<br>    Space in 3 variables with 3 components, size 30<br>    Sum space of 3 concatenated subspaces (all identical)<br>      PetscSpace Object: sum component (sumcomp_) 1 MPI processes<br>        type: poly<br>        Space in 3 variables with 1 components, size 10<br>        Polynomial space of degree 2<br>  PetscDualSpace Object: P2 1 MPI processes<br>    type: lagrange<br>    Dual space with 3 components, size 30<br>    Continuous Lagrange dual space<br>    Quadrature of order 5 on 27 points (dim 3)<br>New Bounding Box:<br>  0: lo = 2.5624e-17 hi = 8.<br>  1: lo = -9.23372e-17 hi = 7.<br>  2: lo = 2.72091e-17 hi = 8.5<br><div>``` </div></div></div><div><br></div><div>Thanks,</div><div>Zongze</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Zongze Yang <<a href="mailto:yangzongze@gmail.com">yangzongze@gmail.com</a>> 于2022年6月17日周五 14:54写道：<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">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.<div><br></div><div>First, I patch the petsc4py by adding `DMProjectCoordinates`:</div><div>```</div><div>diff --git a/src/binding/petsc4py/src/PETSc/DM.pyx b/src/binding/petsc4py/src/PETSc/DM.pyx<br>index d8a58d183a..dbcdb280f1 100644<br>--- a/src/binding/petsc4py/src/PETSc/DM.pyx<br>+++ b/src/binding/petsc4py/src/PETSc/DM.pyx<br>@@ -307,6 +307,12 @@ cdef class DM(Object):<br>         PetscINCREF(c.obj)<br>         return c<br><br>+    def projectCoordinates(self, FE fe=None):<br>+        if fe is None:<br>+            CHKERR( DMProjectCoordinates(<a href="http://self.dm" target="_blank">self.dm</a>, NULL) )<br>+        else:<br>+            CHKERR( DMProjectCoordinates(<a href="http://self.dm" target="_blank">self.dm</a>, fe.fe) )<br>+<br>     def getBoundingBox(self):<br>         cdef PetscInt i,dim=0<br>         CHKERR( DMGetCoordinateDim(<a href="http://self.dm" target="_blank">self.dm</a>, &dim) )<br>diff --git a/src/binding/petsc4py/src/PETSc/petscdm.pxi b/src/binding/petsc4py/src/PETSc/petscdm.pxi<br>index 514b6fa472..c778e39884 100644<br>--- a/src/binding/petsc4py/src/PETSc/petscdm.pxi<br>+++ b/src/binding/petsc4py/src/PETSc/petscdm.pxi<br>@@ -90,6 +90,7 @@ cdef extern from * nogil:<br>     int DMGetCoordinateDim(PetscDM,PetscInt*)<br>     int DMSetCoordinateDim(PetscDM,PetscInt)<br>     int DMLocalizeCoordinates(PetscDM)<br>+    int DMProjectCoordinates(PetscDM, PetscFE)<br><br>     int DMCreateInterpolation(PetscDM,PetscDM,PetscMat*,PetscVec*)<br>     int DMCreateInjection(PetscDM,PetscDM,PetscMat*)<br></div><div>```</div><div><br></div><div>Then in python, I load a mesh and project the coordinates to P2:</div><div><font face="arial, sans-serif">```</font></div><div><font face="arial, sans-serif">import firedrake as fd<br>from firedrake.petsc import PETSc<br><br># plex = fd.mesh._from_gmsh('test-fd-load-p2.msh')<br>plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh')<br>print('old bbox:', plex.getBoundingBox())<br><br>dim = plex.getDimension()<br>#                                 (dim,  nc, isSimplex, k,          qorder, comm=None)<br>fe_new = PETSc.FE().createLagrange(dim, dim,      True, 2, PETSc.DETERMINE)<br>plex.projectCoordinates(fe_new)<br>fe_new.view()<br><br>print('new bbox:', plex.getBoundingBox())<br></font></div><div><font face="arial, sans-serif">```</font></div><div><font face="arial, sans-serif"><br></font></div><div><font face="arial, sans-serif">The output is (The bounding box is changed!)</font></div><div><font face="arial, sans-serif">```</font></div><div><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif">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))</font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif">```</font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif"><br></font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif">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?</font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif"><br></font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif">Thanks!</font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif"><br></font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><font face="arial, sans-serif">  Zongze</font></pre><pre style="box-sizing:unset;border:none;margin-top:0px;margin-bottom:0px;padding:0px;overflow:auto;word-break:break-all;white-space:pre-wrap"><br></pre></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> 于2022年6月17日周五 01:11写道：<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang <<a href="mailto:yangzongze@gmail.com" target="_blank">yangzongze@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto"><div dir="ltr"></div><div dir="ltr"><br></div><div dir="ltr"><br><blockquote type="cite">在 2022年6月16日，23:22，Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> 写道：<br><br></blockquote></div><blockquote type="cite"><div dir="ltr"><div dir="ltr"><div dir="ltr">On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang <<a href="mailto:yangzongze@gmail.com" target="_blank">yangzongze@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>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? </div></div></blockquote><div><br></div><div>By default, they are stored as P2, not DG.</div></div></div></div></blockquote><div><br></div><div>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.</div><div>Then the function <a href="https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexVecGetClosure/" target="_blank">DMPlexVecGetClosure</a> seems return the coordinates in lex order.</div><div><br></div><div>Some code in reading gmsh file reads that</div><div><br></div><div><br></div><div><p style="margin:0px;font-stretch:normal;font-size:medium;line-height:normal;color:rgb(0,0,0)"><span style="font-size:16px">1756:    <span> </span>if (isSimplex) continuity = <a href="https://petsc.org/main/docs/manualpages/Sys/PETSC_FALSE/" target="_blank"><span>PETSC_FALSE</span></a>; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */</span></p><p style="margin:0px;font-stretch:normal;font-size:medium;line-height:normal;min-height:21px;color:rgb(0,0,0)"><span style="font-size:16px"></span><br></p><p style="margin:0px;font-stretch:normal;font-size:medium;line-height:normal;color:rgb(0,0,0)"><span style="font-size:16px">1758:    <span> </span>GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe)</span></p><p style="margin:0px;font-stretch:normal;font-size:medium;line-height:normal;color:rgb(0,0,0)"><span style="font-size:16px"><br></span></p><p style="margin:0px;font-stretch:normal;font-size:medium;line-height:normal;color:rgb(0,0,0)"><span style="font-size:16px">The continuity is set to false for simplex.</span></p></div></div></blockquote><div><br></div><div>Oh, yes. That needs to be fixed. For now, you can just project it to P2 if you want using</div><div><br></div><div>  <a href="https://petsc.org/main/docs/manualpages/DM/DMProjectCoordinates/" target="_blank">https://petsc.org/main/docs/manualpages/DM/DMProjectCoordinates/</a></div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto"><div>Thanks,</div><div>Zongze</div><blockquote type="cite"><div dir="ltr"><div dir="ltr"><div class="gmail_quote"><div>You can ask for the coordinates of a vertex or an edge directly using</div><div><br></div><div>  <a href="https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexPointLocalRead/" target="_blank">https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexPointLocalRead/</a></div><div><br></div><div>by giving the vertex or edge point. You can get all the coordinates on a cell, in the closure order, using</div><div><br></div><div>  <a href="https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexVecGetClosure/" target="_blank">https://petsc.org/main/docs/manualpages/DMPLEX/DMPlexVecGetClosure/</a></div></div></div></div></blockquote><blockquote type="cite"><div dir="ltr"><div dir="ltr"><div class="gmail_quote"><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Below is some code load the gmsh file, I want to know the relation between `cl` and `cell_coords`.</div><div><br></div>```<div>import firedrake as fd<br>import numpy as np<br><br># Load gmsh file (2rd)<br>plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh')<br><br>cs, ce = plex.getHeightStratum(0)<br><br>cdm = plex.getCoordinateDM()<br>csec = dm.getCoordinateSection()<br>coords_gvec = dm.getCoordinates()<br><br>for i in range(cs, ce):<br>    cell_coords = cdm.getVecClosure(csec, coords_gvec, i)<br>    print(f'coordinates for cell {i} :\n{cell_coords.reshape([-1, 3])}')<br>    cl = dm.getTransitiveClosure(i)<br>    print('closure:', cl)<br>    break</div><div>```<br clear="all"><div><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><br></div><div dir="ltr">Best wishes,</div></div><div>Zongze</div></div></div></div></div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</div></blockquote></div></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div>
</blockquote></div>