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