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? Additionally, should we always set `dm_plex_gmsh_project_petscdualspace_lagrange_continuity` to false for high-order gmsh files? With the option set to `true`, I got the following error: ``` $ $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/> >
