how to import a 2nd order mesh?

[APaganini]

The minimal example

from firedrake import * mesh = Mesh("disc.msh")

fails with message

Error: error code 62 [0] DMPlexCreateGmsh() line 200 in /private/tmp/pip-8t0buL-build/src/dm/impls/plex/plexgmsh.c [0] DMPlexCreateGmsh_ReadElement() line 385 in /private/tmp/pip-8t0buL-build/src/dm/impls/plex/plexgmsh.c [0] Invalid argument [0] Unsupported Gmsh element type 8

when disc.msh is generated with the command gmsh disc.geo -2 -order 2,
which creates a "2nd order" mesh (no problem arises when disc.msh is generated
with gmsh disc.geo -2 -order 1). The file disc.geo describes a disc; its content reads

Point(1) = {0, 0, 0, 1.0};
Point(2) = {-0.5, -0, 0, 1.0};
Point(3) = {0, -0.5, 0, 1.0};
Point(4) = {0, 0.5, 0, 1.0};
Point(5) = {0.5, -0, 0, 1.0};
Characteristic Length { 1,2,3,4,5 } = 0.25;

Circle(1) = {2, 1, 3}; 
Circle(2) = {3, 1, 5}; 
Circle(3) = {5, 1, 4}; 
Circle(4) = {4, 1, 2}; 
Line Loop(5) = {4, 1, 2, 3}; 
Plane Surface(6) = {5};

Physical Surface("MyDisk") = {6};
Physical Line("Bdry") = {4, 1, 2, 3}; 

[wence-]

The PETSc gmsh reader doesn't know about second order (or higher) coordinate fields.

I thought that @mlange05 had done some work on this for nektar, but it doesn't appear to have been merged in PETSc. If that was available we would just need to do a little work in firedrake to recognise these cases.

[APaganini]

Thanks for the reply. For the moment I can try a workaround by interpolating
a transformation field. However, this won't be general&automatic, as it requires
knowledge on the geometry.

To emphasize the relevance of higher order meshes, let me mention that
2nd order FEM loses half an order of convergence (wrt mesh refinement) if a
curved boundary is approximated with a polytope. This is a problem when
you're doing convergence experiments.

[APaganini]

It seems that DMPlex has been updated. Now the error message reads:

Traceback (most recent call last):
  File "importmesh.py", line 2, in <module>
    mesh = Mesh("disc.msh")
  File "<decorator-gen-262>", line 2, in Mesh
  File "/Volumes/Scratch/Users/paalbert/Documents/FIREDRAKE/firedrakeAdjoint/src/PyOP2/pyop2/profiling.py", line 61, in wrapper
    return f(*args, **kwargs)
  File "/Volumes/Scratch/Users/paalbert/Documents/FIREDRAKE/firedrakeAdjoint/src/firedrake/firedrake/mesh.py", line 1223, in Mesh
    topology = MeshTopology(plex, name=name, reorder=reorder, distribute=distribute)
  File "<decorator-gen-260>", line 2, in __init__
  File "/Volumes/Scratch/Users/paalbert/Documents/FIREDRAKE/firedrakeAdjoint/src/PyOP2/pyop2/profiling.py", line 61, in wrapper
    return f(*args, **kwargs)
  File "/Volumes/Scratch/Users/paalbert/Documents/FIREDRAKE/firedrakeAdjoint/src/firedrake/firedrake/mesh.py", line 349, in __init__
    dmplex.validate_mesh(plex)
  File "firedrake/dmplex.pyx", line 1009, in firedrake.dmplex.validate_mesh (firedrake/dmplex.c:14452)
ValueError: Provided mesh has some entities not reachable by traversing cells (maybe rogue vertices?)

[wence-]

That check is traversing the mesh by cell and then gathering all the mesh points in the closure of the cells. We require in Firedrake that union(c, closure(c)) == mesh. I.e. all mesh points can be visited by walking the closure of some cell. Sometimes when you make a gmsh mesh, that is not the case. Have you done the necessary magic with Physical Line/Point etc...?

[APaganini]

the *.geo file is in my first post (see above). note that there are no errors if the mesh is of order 1

[wence-]

In this case, it looks like somehow things that dmplex reads as "facets" are not in the closure of the cells.

[wence-]

e.g.

ValueError: Provided mesh has some entities not reachable by traversing cells (maybe rogue vertices?)

In [3]: %debug
> /data/lmitche1/src/doodles/firedrake/dmplex.pyx(1009)firedrake.dmplex.validate_mesh (firedrake/dmplex.c:14452)()

ipdb> up
> /data/lmitche1/src/firedrake/src/firedrake/firedrake/mesh.py(349)__init__()
    347         """
    348         # Do some validation of the input mesh
--> 349         dmplex.validate_mesh(plex)
    350         utils._init()
    351 

ipdb> plex.getChart()
(0, 313)
ipdb> plex.getDepthStratum(0)
(64, 209)
ipdb> plex.getDepthStratum(1)
(209, 313)
ipdb> plex.getDepthStratum(2)
(0, 64)
ipdb> import numpy
ipdb> !a = numpy.zeros(313, dtype=int)
ipdb> !for c in range(64): a[plex.getTransitiveClosure(c)[0]] = 1
ipdb> print a
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0
 0 1 1 1 0 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
ipdb> a[:64]
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
ipdb> a[64:209]
array([1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0,
       0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0])
ipdb> a[209:313]
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
ipdb> numpy.where(a == 0)
(array([ 71,  72,  73,  74,  78,  79,  80,  81,  85,  86,  87,  88,  92,
        93,  94,  95, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130,
       131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
       144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156,
       157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,
       170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182,
       183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195,
       196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208]),)
ipdb> plex.getSupport(71)
array([], dtype=int32)
ipdb> plex.getSupport(72)
array([], dtype=int32)
ipdb> plex.getCone(72)
array([], dtype=int32)
ipdb> 

[wence-]

So I think you need to fix the mesh reading first...

[APaganini]

I literally have no idea how to do that. The DMPlexCreateGmsh-code is fairly obscure to me, and I do not see where a field is generated (if it is generated at all). Also, it is not quite clear to me what

ipdb> !for c in range(64): a[plex.getTransitiveClosure(c)[0]] = 1
ipdb> print a

returns :(

[wence-]

So this is telling you that there are 313 mesh points in the mesh. But not all of them are reachable by traversing the closure of the cells. As evidenced by some of those entries being 0. In particular a bunch of the facets have no support (i.e. they are not connected to a cell). If the terminology is confusing, I recommend Matt's notes http://www.caam.rice.edu/~caam519/CSBook.pdf, particularly chapter 7.

[danshapero]

@wence- I'm running into a similar issue but your link to Matt's notes has gone dead -- do you have an updated link?

[wence-]

Updated link (Matt is now at Buffalo): https://cse.buffalo.edu/~knepley/classes/caam519/CSBook.pdf

[danshapero]

I now get a different error when trying Alberto's mesh. I can lead the order-2 mesh, but when I try to read the .coordinate field, I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [9], line 1
----> 1 mesh.coordinates

File ~/.pyenv/versions/firedrake-python3.10/src/firedrake/firedrake/mesh.py:1579, in MeshGeometry.coordinates(self)
   1576 @property
   1577 def coordinates(self):
   1578     """The :class:`.Function` containing the coordinates of this mesh."""
-> 1579     return self._coordinates_function

File ~/.pyenv/versions/firedrake-python3.10/src/firedrake/firedrake/adjoint/mesh.py:30, in MeshGeometryMixin._ad_annotate_coordinates_function.<locals>.wrapper(self, *args, **kwargs)
     27 @wraps(coordinates_function)
     28 def wrapper(self, *args, **kwargs):
     29     from .blocks import MeshInputBlock, MeshOutputBlock
---> 30     f = coordinates_function(self)
     31     f.block_class = MeshInputBlock
     32     f._ad_floating_active = True

File ~/.pyenv/versions/firedrake-python3.10/src/firedrake/firedrake/mesh.py:1569, in MeshGeometry._coordinates_function(self)
   1567     return self.ufl_cargo()._coordinates_function
   1568 else:
-> 1569     self.init()
   1570     coordinates_fs = self._coordinates.function_space()
   1571     V = functionspaceimpl.WithGeometry.create(coordinates_fs, self)

File ~/.pyenv/versions/firedrake-python3.10/src/firedrake/firedrake/mesh.py:1481, in MeshGeometry.init(self)
   1476 """Finish the initialisation of the mesh.  Most of the time
   1477 this is carried out automatically, however, in some cases (for
   1478 example accessing a property of the mesh directly after
   1479 constructing it) you need to call this manually."""
   1480 if hasattr(self, '_callback'):
-> 1481     self._callback(self)

File ~/.pyenv/versions/firedrake-python3.10/src/firedrake/firedrake/mesh.py:1504, in MeshGeometry._init_topology.<locals>.callback(self)
   1502 self.topology.init()
   1503 coordinates_fs = functionspace.FunctionSpace(self.topology, self.ufl_coordinate_element())
-> 1504 coordinates_data = dmcommon.reordered_coords(topology.topology_dm, coordinates_fs.dm.getDefaultSection(),
   1505                                              (self.num_vertices(), self.ufl_coordinate_element().cell().geometric_dimension()))
   1506 coordinates = function.CoordinatelessFunction(coordinates_fs,
   1507                                               val=coordinates_data,
   1508                                               name=_generate_default_mesh_coordinates_name(self.name))
   1509 self.__init__(coordinates)

File ~/.pyenv/versions/firedrake-python3.10/src/firedrake/firedrake/cython/dmcommon.pyx:1383, in firedrake.cython.dmcommon.reordered_coords()

ValueError: cannot reshape array of size 768 into shape (41,2)

I think PETSc can read higher-order meshes from gmsh now and the problem might just be some impedance mismatch on the Firedrake side?

[wence-]

We don't inspect the PETSc object to determine the element for the coordinate space (so VectorElement("P", mesh.ufl_cell(), 1) is hard-coded). Does it work if you hand-hack that to say 2 rather than 1? (I suspect not because the coordinate unscrambling will be wrong, but perhaps yes)

[danshapero]

I get the same error when I manually change 1 to 2 at this line in mesh.py.

That line is definitely going to be a problem because it hard-codes the element degree. This code is also going to hurt because the assumed shape is (num_vertices, geometric_dimension), but that's no longer true for higher-order coordinate fields.

Do you know how to sniff the degree of the coordinate field from the DMPlex? Can't seem to find this in the PETSc docs.

[wence-]

I get the same error when I manually change 1 to 2 at this line in mesh.py.

That line is definitely going to be a problem because it hard-codes the element degree. This code is also going to hurt because the assumed shape is (num_vertices, geometric_dimension), but that's no longer true for higher-order coordinate fields.

Do you know how to sniff the degree of the coordinate field from the DMPlex? Can't seem to find this in the PETSc docs.

Let's ask @knepley. I think you have to pull the PetscDS out of the coordinateDM and inspect that.

[knepley]

@wence- is correct, but we can just use the Field directly here. The sequence would be something like

DMGetCoordinateDM(dm, &cdm);
DMGetField(cdm, 0, &cfe);
PetscFEGetBasisSpace(cfe, &cbasis);
PetscSpaceGetDegree(cbasis, &degree, NULL);

[danshapero]

Ok, I used the Python equivalent of the code that @knepley wrote to get the degree and made some changes on this branch. When I manually create the plex for Alberto's mesh, I'm able to get out the FE, basis space, and degree. For some of the utility meshes I only get a PETSc Object from .getField, which I don't understand, but I was able to work around that with a try/except.

Something's still not right though. For piecewise parabolic functions on triangles, there should be one degree of freedom for every vertex and one for every edge. Alberto's mesh has 41 vertices and 104 edges. I'm still hitting the same error as before only now it's unable to reshape an array of size 768 into an array of size (145, 2), whereas before when it was only counting the number of vertices that was an array of size (41, 2). I don't understand how the coordinate field from the PETSc side has 768 elements if this is the number of mesh vertices and edges.

[knepley]

@danshapero The reason you get a PetscObject sometimes is that Plex allows you not to create an FE space for the coordinates, and just falls back to linear. I felt I needed to do this when FE was very immature. Now it seems like a mistake.

We should be able to see what it thinks about coordinates by printing the section. So DMGetLocalSection(cdm, &cs) and then view it.

[danshapero]

Gotcha, thanks for clarifying -- the code that I have falls back to linear when the result is PetscObject instead of an FE. So whether or not this suggests some change to PETSc it's good to know that it's at least doing the right thing.

The following code

from firedrake.petsc import PETSc
viewer = PETSc.Viewer().create()
viewer.setType("ascii")
viewer.setFileMode("r")
viewer.setFileName("disc.msh")
plex = PETSc.DMPlex().createGmsh(viewer)
print(plex.getDepthStratum(0), plex.getDepthStratum(1), plex.getDepthStratum(2), end="\n\n")
dm = plex.getCoordinateDM()
section = dm.getSection()
section.view()

gives me this:

(64, 105), (105, 209), (0, 64)

PetscSection Object: PetscSection_0x84000002_2 1 MPI process
  type not yet set
1 fields
  field 0 with 2 components
Process 0:
  (   0) dim 12 offset   0
  (   1) dim 12 offset  12
...
  (  62) dim 12 offset 744
  (  63) dim 12 offset 756
  (  64) dim  0 offset 768
  (  65) dim  0 offset 768
...
  ( 207) dim  0 offset 768
  ( 208) dim  0 offset 768
  PetscSectionSym Object: 1 MPI process
    type: label
    Label 'depth'
    Symmetry for stratum value 0 (0 dofs per point): no symmetries
    Symmetry for stratum value 1 (0 dofs per point): no symmetries
    Symmetry for stratum value 2 (12 dofs per point):
      Orientation range: [-3, 3)
    Symmetry for stratum value -1 (0 dofs per point): no symmetries

So it looks like there are 12 degrees of freedom for each of the entities in the 2nd depth stratum. All of this suggests that internally the plex is storing everything on the triangles and not on the vertices or edges (64 triangles * 6 points per triangle for parabolics * 2 DoFs per point = 768). Am I reading this right?

[knepley]

That is the DG coordinate field, so the mesh is periodic. However, Plex will not give this back for dm.getCoordinates(), you need dm.getCellCoordinates().

[danshapero]

I don't understand -- this shouldn't be a periodic mesh at all. When I generate an order-1 mesh, all the degrees of freedom for the coordinate field are clearly located at the vertices. Is the code that I wrote above just accessing a DG view of the coordinates and I have to do something different to access them with a different order?

DMGetCellCoordinates and the related functions aren't wrapped in petsc4py. I added them and tried using them on the order-2 disc mesh but I just get a seg fault whenever I try and inspect the Vec that I get out of getCellCoordinates when I use it either on the plex itself or on the coordinate DM.

Also, if this is all written down some place please let me know -- there's something I'm missing that falls between the high-level overview of DMPlex and the API reference in the PETSc docs.

[knepley]

I don't understand either. Plex only makes a DG coordinate field for a periodic mesh, and even in that case, the default coordinate field is CG and you have to ask specifically for the DG one. Maybe we need to Zoom on this. I will be back in the US on Monday.

[danshapero]

That sounds great, I'll email you to figure out a time. Meanwhile I wrote to the PETSc mailing list to see if anyone there can clear this up for me.

[jedbrown]

I replied over on petsc-users:

It's confusing, but this line makes high order simplices always read as discontinuous coordinate spaces. I would love if someone would revisit that, perhaps also using DMPlexSetIsoperiodicFaceSF(), which should simplify the code and avoid the confusing cell coordinates pattern. Sadly, I don't have time to dive in.

https://gitlab.com/petsc/petsc/-/commit/066ea43f7f75752f012be6cd06b6107ebe84cc6d#3616cad8148970af5b97293c49492ff893e25b59_1552_1724

[jedbrown]

An alternative to avoid getting hands dirty would be to do a projection of the discontinuous coordinates into a continuous space.

[knepley]

Yes, I suggested that over on petsc-users.

[danshapero]

Sorry for the cross-posting and thanks @jedbrown for clearing this up! I'll see if we can hack together a workaround on the Firedrake side for now.