Define a FunctionSpace containing only Functions living on a coarser mesh?

Question

Suppose we are given two nested triangulations created via

mesh2 = refine(mesh1)

is it possible to create a piecewise linear space associated with mesh2, but containing only those functions present on the coarser level mesh1?

Let me add that I would willingly trade some unnecessary DOFs for simplicity of a potential solution, since I only want to do some prototyping right now. Thus, it is perfectly fine if a solution is proposed which only constrains the fine-level DOFs but leaves them in the system.

The question is related to the discussion 'Mixed elements with two different meshes' which did not come up with a solution yet.

A way which doesn't work: Suppose we have the nested function spaces

Q1 = FunctionSpace(mesh1, 'CG', 1)
Q2 = FunctionSpace(mesh2, 'CG', 1)

These spaces have the problem that they live on different meshes. Particularly, I cannot create a mixed space

W = Q1*Q2

and get away with it. From my experiments I found out that it is sort of allowed by the interface that you create and use mixed spaces using components defined on different meshes, but (understandably) things will break whenever you try to solve some problem defined using such a mixed space.

Some background: If it is not yet clear what I mean, let me try to give some more detail. I want to to test a projection-based stabilization technique for equal-order elements for Stokes flow. Therefore, I need to be able to incorporate a term of the form

$$ \int_\Omega (p-Pp)\cdot(q-Pq) dx $$

where P projects pressures from the fine-level to the coarse-level (e.g., by injection).

To my knowledge, something like this cannot be assembled out of the box. The worst case would be to implement the assembly of the matrix-block associated with the stabilization term by hand. I just want to make sure that there is no easier way before starting that, so I thought about possible workarounds

Alternatively to assembling these terms, everything would be fine if I could define a FunctionSpace containing the functions of Q1, but defined on the finer mesh2.

In that case, I would add an auxiliary pressure variable Pp and define these terms as

(p-Pp)*(q-Pq)*dx

with test functions q and Pq from the respective pressure function spaces. This would then ensure that Pp is the projection of p to the coarser level (at the price of additional DOFs for Pp, but I don't care about such issues at the moment) and the discrete problem should be equivalent to having a realization of operator P.

Any help is appreciated!

Answer 1

Okay, after some tedious fiddling, I found part of the solution, which is the placement of velocities and pressures on different levels. The rest should be not a problem anymore. What the code listed below does is to solve the Stokes equations using a Taylor-Hood(ish) element, but instead of the usual $P_2-P_1$-combination it interpolates the velocity using piecewise linears on a finer mesh (sometimes called $P_1-iso-P_2$ element).

What I did to put the pressure on the coarser level is simply to constrain the unnecessary pressures at the midpoint nodes. Be aware that I made some assumptions regarding DOF-numbering etc. and maybe I forgot about some corner pressures when other meshes are used. I consider this a hack... However, I am pretty happy with it, since it is not awfully slow. If somebody has a more elegant/faster/shorter suggestion to do the things below, please let me know!

Example code:

from dolfin import *
import numpy, sys, math

dolfin.parameters["reorder_dofs_serial"] = False
maxlevel = 5

def print_rates(h, E):
  print   'err(u,L2)  rate  err(u,H1)  rate  err(p,L2)  rate'
  print   '%.4e -     %.4e -     %.4e -     ' \
    %(E[0][0], E[0][1], E[0][2])
  for i in range(1, len(E)):
    b = h[i-1]/h[i]
    err0 = math.log(E[i-1][0]/E[i][0], b)
    err1 = math.log(E[i-1][1]/E[i][1], b)
    err2 = math.log(E[i-1][2]/E[i][2], b)
    print '%.4e %.2f  %.4e %.2f  %.4e %.2f' \
      %(E[i][0], err0, E[i][1], err1, E[i][2], err2)
  print '\n'
  return

mesh = UnitSquareMesh(1, 1, 'crossed')
h, err = [], []

for i in xrange(maxlevel):

  nv0 = mesh.size(0) # remember for later
  mesh = refine(mesh)

  V = VectorFunctionSpace(mesh, 'CG', 1)
  Q = FunctionSpace(mesh, 'CG', 1)
  R = FunctionSpace(mesh, 'R', 0)
  W = MixedFunctionSpace([V, Q, R])

  # variational problem
  u, p, m = TrialFunctions(W)
  v, q, r = TestFunctions(W)

  # exact solution and rhs
  U = Expression(('x[0]+x[0]*x[0]-2*x[0]*x[1]+x[0]*x[0]*x[0]-3*x[0]*x[1]*x[1]+x[0]*x[0]*x[1]',\
                  '-x[1]-2*x[0]*x[1]+x[1]*x[1]-3*x[0]*x[0]*x[1]+x[1]*x[1]*x[1]-x[0]*x[1]*x[1]'), degree=4)
  P = Expression('x[0]*x[1]+x[0]+x[1]+x[0]*x[0]*x[0]*x[1]*x[1]-4./3', degree=5)
  F = Expression(('3*x[0]*x[0]*x[1]*x[1]-x[1]-1.0', '2*x[0]*x[0]*x[0]*x[1]+3*x[0]-1.0'), degree=4)

  bc = DirichletBC(W.sub(0), U, 'on_boundary')

  stokes = inner(grad(u),grad(v))*dx - div(v)*p*dx - div(u)*q*dx - dot(F,v)*dx + p*r*dx + m*q*dx
  A, b = assemble_system(lhs(stokes), rhs(stokes), bc, finalize_tensor = False)

  # constrain some dofs
  dm = W.sub(1).dofmap().collapse(mesh)[1]
  map = {}
  for v in vertices(mesh):
    if v.index() >= nv0 or len([c for c in cells(v)]) == 1:
      vdof = dm[v.index()]
      map[vdof] = []
      for e in edges(v):
        map[vdof] += [dm[v2.index()] for v2 in vertices(e) if v2.index() < nv0]
  rows = numpy.asarray([m[0] for m in map.items()], dtype='intc')
  for r in rows:
    A.setrow(r, numpy.asarray([r] + map[r], dtype='uintp'), numpy.asarray([1.0,-0.5,-0.5]))
  A.apply('insert')

  # solve the problem
  x = Function(W)
  solve(A, x.vector(), b)
  uh, ph, mh = x.split()
  plot(ph)

  h.append(mesh.hmin())
  err.append((errornorm(U, uh), errornorm(U, uh, 'H10'), errornorm(P, ph)))

print_rates(h,err)

Output:

err(u,L2)  rate  err(u,H1)  rate  err(p,L2)  rate
2.5905e-01 -     2.1610e+00 -     2.2103e+00 -     
6.7218e-02 1.95  9.5323e-01 1.18  8.6163e-01 1.36
1.6748e-02 2.00  4.4739e-01 1.09  2.7860e-01 1.63
4.2041e-03 1.99  2.1694e-01 1.04  8.4477e-02 1.72
1.0568e-03 1.99  1.0713e-01 1.02  2.4851e-02 1.77