Got it. Thank you for your explanation! Best wishes, Zongze
On Mon, 15 May 2023 at 23:28, Matthew Knepley <[email protected]> wrote: > On Mon, May 15, 2023 at 9:55 AM Zongze Yang <[email protected]> wrote: > >> On Mon, 15 May 2023 at 17:24, Matthew Knepley <[email protected]> wrote: >> >>> On Sun, May 14, 2023 at 7:23 PM Zongze Yang <[email protected]> >>> wrote: >>> >>>> Could you try to project the coordinates into the continuity space by >>>> enabling the option >>>> `-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true`? >>>> >>> >>> There is a comment in the code about that: >>> >>> /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */ >>> >>> So what is currently done is you project into the discontinuous space >>> from the GMsh coordinates, >>> and then we get the continuous coordinates from those later. This is why >>> we get the right answer. >>> >>> >> Sorry, I'm having difficulty understanding the comment and fully >> understanding your intended meaning. Are you saying that we can only >> project the space to a discontinuous space? >> > > For higher order simplices, because we do not have the mapping to the GMsh > order yet. > > >> Additionally, should we always set >> `dm_plex_gmsh_project_petscdualspace_lagrange_continuity` to false for >> high-order gmsh files? >> > > This is done automatically if you do not override it. > > >> With the option set to `true`, I got the following error: >> > > Yes, do not do that. > > Thanks, > > Matt > > >> ``` >> $ $PETSC_DIR/$PETSC_ARCH/tests/dm/impls/plex/tests/runex33_gmsh_3d_q2.sh >> -e "-dm_plex_gmsh_project_petscdualspace_lagrange_continuity true" >> not ok dm_impls_plex_tests-ex33_gmsh_3d_q2 # Error code: 77 >> # Volume: 0.46875 >> # [0]PETSC ERROR: --------------------- Error Message >> -------------------------------------------------------------- >> # [0]PETSC ERROR: Petsc has generated inconsistent data >> # [0]PETSC ERROR: Calculated volume 0.46875 != 1. actual volume >> (error 0.53125 > 1e-06 tol) >> # [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble >> shooting. >> # [0]PETSC ERROR: Petsc Development GIT revision: >> v3.19.1-294-g9cc24bc9b93 GIT Date: 2023-05-15 12:07:10 +0000 >> # [0]PETSC ERROR: ../ex33 on a arch-linux-c-debug named AMA-PC-RA18 >> by yzz Mon May 15 21:53:43 2023 >> # [0]PETSC ERROR: Configure options >> --CFLAGS=-I/opt/intel/oneapi/mkl/latest/include >> --CXXFLAGS=-I/opt/intel/oneapi/mkl/latest/include >> --LDFLAGS="-Wl,-rpath,/opt/intel/oneapi/mkl/latest/lib/intel64 >> -L/opt/intel/oneapi/mkl/latest/lib/intel64" --download-bison >> --download-chaco --download-cmake >> --download-eigen="/home/yzz/firedrake/complex-int32-mkl-X-debug/src/eigen-3.3.3.tgz >> " --download-fftw --download-hdf5 --download-hpddm --download-hwloc >> --download-libpng --download-metis --download-mmg --download-mpich >> --download-mumps --download-netcdf --download-p4est --download-parmmg >> --download-pastix --download-pnetcdf --download-ptscotch >> --download-scalapack --download-slepc --download-suitesparse >> --download-superlu_dist --download-tetgen --download-triangle >> --with-blaslapack-dir=/opt/intel/oneapi/mkl/latest --with-c2html=0 >> --with-debugging=1 --with-fortran-bindings=0 >> --with-mkl_cpardiso-dir=/opt/intel/oneapi/mkl/latest >> --with-mkl_pardiso-dir=/opt/intel/oneapi/mkl/latest >> --with-scalar-type=complex --with-shared-libraries=1 --with-x=1 --with-zlib >> PETSC_ARCH=arch-linux-c-debug >> # [0]PETSC ERROR: #1 CheckVolume() at >> /home/yzz/opt/petsc/src/dm/impls/plex/tests/ex33.c:246 >> # [0]PETSC ERROR: #2 main() at >> /home/yzz/opt/petsc/src/dm/impls/plex/tests/ex33.c:261 >> # [0]PETSC ERROR: PETSc Option Table entries: >> # [0]PETSC ERROR: -coord_space 0 (source: command line) >> # [0]PETSC ERROR: -dm_plex_filename >> /home/yzz/opt/petsc/share/petsc/datafiles/meshes/cube_q2.msh (source: >> command line) >> # [0]PETSC ERROR: -dm_plex_gmsh_project (source: command line) >> # [0]PETSC ERROR: >> -dm_plex_gmsh_project_petscdualspace_lagrange_continuity true (source: >> command line) >> # [0]PETSC ERROR: -tol 1e-6 (source: command line) >> # [0]PETSC ERROR: -volume 1.0 (source: command line) >> # [0]PETSC ERROR: ----------------End of Error Message -------send >> entire error message to [email protected] >> # application called MPI_Abort(MPI_COMM_SELF, 77) - process 0 >> ok dm_impls_plex_tests-ex33_gmsh_3d_q2 # SKIP Command failed so no diff >> ``` >> >> >> Best wishes, >> Zongze >> >> Thanks, >>> >>> Matt >>> >>> >>>> Best wishes, >>>> Zongze >>>> >>>> >>>> On Mon, 15 May 2023 at 04:24, Matthew Knepley <[email protected]> >>>> wrote: >>>> >>>>> On Sun, May 14, 2023 at 12:27 PM Zongze Yang <[email protected]> >>>>> wrote: >>>>> >>>>>> >>>>>> >>>>>> >>>>>> On Sun, 14 May 2023 at 23:54, Matthew Knepley <[email protected]> >>>>>> wrote: >>>>>> >>>>>>> On Sun, May 14, 2023 at 9:21 AM Zongze Yang <[email protected]> >>>>>>> wrote: >>>>>>> >>>>>>>> Hi, Matt, >>>>>>>> >>>>>>>> The issue has been resolved while testing on the latest version of >>>>>>>> PETSc. It seems that the problem has been fixed in the following merge >>>>>>>> request: https://gitlab.com/petsc/petsc/-/merge_requests/5970 >>>>>>>> >>>>>>> >>>>>>> No problem. Glad it is working. >>>>>>> >>>>>>> >>>>>>>> I sincerely apologize for any inconvenience caused by my previous >>>>>>>> message. However, I would like to provide you with additional >>>>>>>> information >>>>>>>> regarding the test files. Attached to this email, you will find two >>>>>>>> Gmsh >>>>>>>> files: "square_2rd.msh" and "square_3rd.msh." These files contain >>>>>>>> high-order triangulated mesh data for the unit square. >>>>>>>> >>>>>>>> ``` >>>>>>>> $ ./ex33 -coord_space 0 -dm_plex_filename square_2rd.msh >>>>>>>> -dm_plex_gmsh_project >>>>>>>> -dm_plex_gmsh_project_petscdualspace_lagrange_continuity true >>>>>>>> -dm_plex_gmsh_project_fe_view -volume 1 >>>>>>>> PetscFE Object: P2 1 MPI process >>>>>>>> type: basic >>>>>>>> Basic Finite Element in 2 dimensions with 2 components >>>>>>>> PetscSpace Object: P2 1 MPI process >>>>>>>> type: sum >>>>>>>> Space in 2 variables with 2 components, size 12 >>>>>>>> Sum space of 2 concatenated subspaces (all identical) >>>>>>>> PetscSpace Object: sum component (sumcomp_) 1 MPI process >>>>>>>> type: poly >>>>>>>> Space in 2 variables with 1 components, size 6 >>>>>>>> Polynomial space of degree 2 >>>>>>>> PetscDualSpace Object: P2 1 MPI process >>>>>>>> type: lagrange >>>>>>>> Dual space with 2 components, size 12 >>>>>>>> Continuous Lagrange dual space >>>>>>>> Quadrature on a triangle of order 5 on 9 points (dim 2) >>>>>>>> Volume: 1. >>>>>>>> $ ./ex33 -coord_space 0 -dm_plex_filename square_3rd.msh >>>>>>>> -dm_plex_gmsh_project >>>>>>>> -dm_plex_gmsh_project_petscdualspace_lagrange_continuity true >>>>>>>> -dm_plex_gmsh_project_fe_view -volume 1 >>>>>>>> PetscFE Object: P3 1 MPI process >>>>>>>> type: basic >>>>>>>> Basic Finite Element in 2 dimensions with 2 components >>>>>>>> PetscSpace Object: P3 1 MPI process >>>>>>>> type: sum >>>>>>>> Space in 2 variables with 2 components, size 20 >>>>>>>> Sum space of 2 concatenated subspaces (all identical) >>>>>>>> PetscSpace Object: sum component (sumcomp_) 1 MPI process >>>>>>>> type: poly >>>>>>>> Space in 2 variables with 1 components, size 10 >>>>>>>> Polynomial space of degree 3 >>>>>>>> PetscDualSpace Object: P3 1 MPI process >>>>>>>> type: lagrange >>>>>>>> Dual space with 2 components, size 20 >>>>>>>> Continuous Lagrange dual space >>>>>>>> Quadrature on a triangle of order 7 on 16 points (dim 2) >>>>>>>> Volume: 1. >>>>>>>> ``` >>>>>>>> >>>>>>>> Thank you for your attention and understanding. I apologize once >>>>>>>> again for my previous oversight. >>>>>>>> >>>>>>> >>>>>>> Great! If you make an MR for this, you will be included on the next >>>>>>> list of PETSc contributors. Otherwise, I can do it. >>>>>>> >>>>>>> >>>>>> I appreciate your offer to handle the MR. Please go ahead and take >>>>>> care of it. Thank you! >>>>>> >>>>> >>>>> I have created the MR with your tests. They are working for me: >>>>> >>>>> https://gitlab.com/petsc/petsc/-/merge_requests/6463 >>>>> >>>>> Thanks, >>>>> >>>>> Matt >>>>> >>>>> >>>>>> Best Wishes, >>>>>> Zongze >>>>>> >>>>>> >>>>>>> Thanks, >>>>>>> >>>>>>> Matt >>>>>>> >>>>>>> >>>>>>>> Best wishes, >>>>>>>> Zongze >>>>>>>> >>>>>>>> >>>>>>>> On Sun, 14 May 2023 at 16:44, Matthew Knepley <[email protected]> >>>>>>>> wrote: >>>>>>>> >>>>>>>>> On Sat, May 13, 2023 at 6:08 AM Zongze Yang <[email protected]> >>>>>>>>> wrote: >>>>>>>>> >>>>>>>>>> Hi, Matt, >>>>>>>>>> >>>>>>>>>> There seem to be ongoing issues with projecting high-order >>>>>>>>>> coordinates from a gmsh file to other spaces. I would like to inquire >>>>>>>>>> whether there are any plans to resolve this problem. >>>>>>>>>> >>>>>>>>>> Thank you for your attention to this matter. >>>>>>>>>> >>>>>>>>> >>>>>>>>> Yes, I will look at it. The important thing is to have a good >>>>>>>>> test. Here are the higher order geometry tests >>>>>>>>> >>>>>>>>> >>>>>>>>> https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/tests/ex33.c >>>>>>>>> >>>>>>>>> I take shapes with known volume, mesh them with higher order >>>>>>>>> geometry, and look at the convergence to the true volume. Could you >>>>>>>>> add a >>>>>>>>> GMsh test, meaning the .msh file and known volume, and I will fix it? >>>>>>>>> >>>>>>>>> Thanks, >>>>>>>>> >>>>>>>>> Matt >>>>>>>>> >>>>>>>>> >>>>>>>>>> >>>>>>>>>> Best wishes, >>>>>>>>>> Zongze >>>>>>>>>> >>>>>>>>>> >>>>>>>>>> On Sat, 18 Jun 2022 at 20:31, Zongze Yang <[email protected]> >>>>>>>>>> wrote: >>>>>>>>>> >>>>>>>>>>> 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/> >>>>>>>>>>>> >>>>>>>>>>> >>>>>>>>> >>>>>>>>> -- >>>>>>>>> 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/> >>>>> >>>> >>> >>> -- >>> 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/> >
