[petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Zongze Yang
yangzongze at gmail.com
Fri Jun 17 01:54:03 CDT 2022
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 <knepley at gmail.com> 于2022年6月17日周五 01:11写道：
> On Thu, Jun 16, 2022 at 12:06 PM Zongze Yang <yangzongze at gmail.com> wrote:
>
>>
>>
>> 在 2022年6月16日，23:22，Matthew Knepley <knepley at gmail.com> 写道：
>>
>>
>> On Thu, Jun 16, 2022 at 11:11 AM Zongze Yang <yangzongze at gmail.com>
>> 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/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220617/8d6b6b11/attachment-0001.html>
More information about the petsc-users
mailing list