Cannot refine periodic mesh

[pefarrell]

Consider the following snippet:

from firedrake import *

mesh = PeriodicUnitSquareMesh(10, 10, direction="both")
mh = MeshHierarchy(mesh, 1)

I expected this to refine the mesh, but instead it yields

Traceback (most recent call last):
  File "bug.py", line 4, in <module>
    mh = MeshHierarchy(mesh, 1)
  File "/home/pefarrell/local/firedrake/firedrake-dev-20180914/src/firedrake/firedrake/mg/mesh.py", line 130, in MeshHierarchy
    m.init()
  File "/home/pefarrell/local/firedrake/firedrake-dev-20180914/src/firedrake/firedrake/mesh.py", line 990, in init
    self._callback(self)
  File "/home/pefarrell/local/firedrake/firedrake-dev-20180914/src/firedrake/firedrake/mesh.py", line 1356, in callback
    (self.num_vertices(), geometric_dim))
  File "firedrake/dmplex.pyx", line 851, in firedrake.dmplex.reordered_coords
ValueError: cannot reshape array of size 1200 into shape (400,2)

It looks like the coordinates on the plex are in 3D (400x3 instead of 400x2)?

[pefarrell]

After playing a lot of jiggery-pokery with the PETSc data structures, the refinement now runs:

from firedrake import *
from firedrake.petsc import PETSc

mesh = PeriodicUnitSquareMesh(100, 100, direction="both")

plex = mesh._plex
gdim = mesh.geometric_dimension()
coorddm = plex.getCoordinateDM()
x = mesh.coordinates
nodesection = x.function_space().dm.getDefaultSection()

dofsection = PETSc.Section().create(comm=mesh.mpi_comm())
dofsection.setNumFields(1)
dofsection.setFieldComponents(0, gdim)
(lo, hi) = nodesection.getChart()
dofsection.setChart(lo, hi)
for point in range(lo, hi):
    dof = nodesection.getDof(point)
    dofsection.setDof(point, dof*gdim)
    off = nodesection.getOffset(point)
    dofsection.setOffset(point, off*gdim)
coorddm.setDefaultSection(dofsection)

with x.dat.vec_ro as x_:
    plex.setCoordinates(x_)
    plex.setCoordinatesLocal(x_)

mh = MeshHierarchy(mesh, 1)

mesh = mh[1] # works with mh[0]
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)
(x, y) = SpatialCoordinate(mesh)
u_exact = cos(2*pi*x)*sin(2*pi*y)
u_ex = Function(V)
u_ex.interpolate(u_exact)
File("output/exact.pvd").write(u_ex)
f = -div(grad(u_exact))

F = inner(grad(u), grad(v))*dx - inner(f, v)*dx
nsp = VectorSpaceBasis(constant=True)
lu = {
      "snes_type": "ksponly",
      "ksp_type": "preonly",
      "pc_type": "lu",
      "pc_factor_mat_solver_type": "mumps"
     }
mg = {
      "snes_type": "ksponly",
      "ksp_type": "cg",
      "ksp_monitor_true_residual": None,
      "pc_type": "mg",
      "mg_coarse_ksp_type": "preonly",
      "mg_coarse_pc_type": "lu",
     }
solve(F == 0, u, solver_parameters=lu, nullspace=nsp)

File("output/solution.pvd").write(u)
u.assign(u - u_ex)
File("output/error.pvd").write(u)

However, the refined mesh looks screwed in paraview and solving any PDE on it doesn't work, even with LU, and with a solver that works on the coarse grid.

(but at least the arrays are all the right size, and there's no segfaults!)

[wence-]

When we turn a DM into a Mesh, we currently assume that we're getting P1 coordinates from the DM. I guess in this case we're not.

But there are bunch of things that worry me:

1. The refined plex with DG coordinates should have more dofs than fit in a P1 function space
2. I don't know if the ordering that plex applies in this case matches what we want for DG spaces.

[colinjcotter]

Patrick, did you find a workaround for this or just give up?