<div dir="ltr"><div dir="ltr">On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang <<a href="mailto:yangzongze@gmail.com">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">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></blockquote><div><br></div><div>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.</div><div>My suspicion is that the space we make when reading in the Gmsh coordinates does not match the values (wrong order).</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="ltr"><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" target="_blank">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>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><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>