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/>
