interpolate only works for functionspaces with degree=1 polynomials

Question

The code below is from http://fenicsproject.org/qa/7607/working-with-boundarymesh?show=7607#q7607, but this question is in its own right:
I have a function that is defined on (a functionspace of a submesh of a boundarymesh of)
a boundary, and I wish to use it in a surface-integral of the whole domain.

from dolfin import *
import numpy as np
parameters["allow_extrapolation"] = True
degree = 2

mesh = UnitSquareMesh(20, 20)
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
leftboundary = LeftBoundary()
boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
leftboundary.mark(boundaries, 1)
ds = Measure('ds')[boundaries]

boundarymesh = BoundaryMesh(mesh, 'exterior')
left_mesh = SubMesh(boundarymesh, leftboundary)
L = FunctionSpace(left_mesh, 'CG', degree)
bf = Function(L)
# here goes some more code to find the correct bf
# but here a dummy function will do
bf.vector()[:] = np.ones(20*degree + 1, dtype=float)

V = FunctionSpace(mesh, 'CG', degree)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v)) * dx
b = interpolate(bf, V) * v * ds(1)
# b = project(bf, V) * v * ds(1)
assemble(b)

The code above works for degree=1, but when I use higher order Lagrange polynomials, I get the error

*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate function at point.
*** Reason:  No matching cell found.
*** Where:   This error was encountered inside Function.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  
*** -------------------------------------------------------------------------

Trying to use project instead of interpolate does not work either. What am I missing here?

Thanks in advance

Comments

Under Fenics 1.6 the error reads

*** Error:   Unable to create mesh entity.
*** Reason:  Mesh entity index -1 out of range [0, 20] for entity of dimension 1.
*** Where:   This error was encountered inside MeshEntity.cpp.

Answer 1

I think the problem is that bf "lives" on L (a boundary mesh) and you are trying to project it onto V, which is defined on the whole domain. I don't know why that even works for degree=1.

Imagine that you have a function $f$ on $[0,1]$ and you only know $f(0)$. How would you define the "projection" onto whatever function space on the domain $[0,1]$?. So if the code works for degree=1, then that is a coincident at best but still meaningless from a mathematical opint of view.

Unfortunately, projecting bf onto L doesn't work either (the projection works but not in combination with the ds(1) stuff)

I think you could try to "extend" bf to the domain appropriately and simply don't work with L. So you could define bf on V, initialize it as zero and then only set values at the boundary vertices. How to do this depends on how you specify your boundary data. If they are all equal to one, then you can simply define the whole function bf=1 on the whole domain. For other possibilities, I would have to see your complete code.

It is still an interesting question, because I think that your approach (together with a projection onto L instead of V) would be much closer to the mathematical definition of the problem than the "workaround" that I suggested. But it seems that such a feature is not yet supported.

Comments

I was also a bit surprised this interpolation would work. In fact it 'extends' the function on the boundary to the whole domain (in my case these functions have to be defined on the boundary, I can not initialize them on V directly).

The most natural solution would be that, if one defines a surface integral, the valid functions in the weak form can be either defined on V, or on any other domain that contains the surface to integrate over (i.e. in this case bf * v * ds(1) would be a valid input).

I noticed there is some improvement on its way:
https://bitbucket.org/fenics-project/dolfin/pull-requests/240/interpolation-improvements/diff

---

Why do you think that you cannot define a function on V. According to this (note: doc is for older fenics version), you could define a helper function on your submesh and then communicate the values via the vertex mapping parent_vertex_indices. It's not pretty but it may work...

---

I am still relatively new to FEniCS so I seem to often underestimate what is possible. If it were possible to define a function on V, the whole domain, and then set as nonzero only the dof's corresponding to L, the boundary domain, that would be avoid interpolation/projection and still give the right surface integral. That is what you mean, no?

The following post lead me to believe that for a submesh of a boundarymesh, this is not possible yet.
http://fenicsproject.org/qa/6810/vertex-mapping-from-submesh-boundarymesh-back-actual-mesh?show=6810#q6810
If you know a way (or can confirm that is it possible), I might look into it a bit more detailed.

---

Well according to the docs it seems possible. But maybe you don't even need such a mapping.

Your integral only integrates over some part of the boundary. So really, it does not matter what values your function has inside the domain. Here is an example:

Suppose you are working on the unit square with coordinate system $(x_1,x_2)$ and you want a boundaryfunction $g(x_2)=x_2^2$ only on the left edge $\Gamma_l$ of the square. Next, in analogy to your example, you define a linear form $l(v):= \int_{\Gamma_l}g(x)v\;ds$. One way to achieve this in your code would be (keeping the rest more or less as it is)

g = Expression('x[1]*x[1]') # note that g is defined on the whole domain
l = g * v * ds(1)

If $g$ cannot be expressed as an Expression, but is the solution of some other problem, then chances are that $g$ is already defined on the whole domain. in that case, it really does not matter what its values are in the interior of the domain as the integral in $l$ only integrates over the boundary anyway.

The only serious Issue I can imagine is that for some reason, $g$ actually is the solution to a lower dimensional problem (I think that such examples can occur when using e.g. domain decomposition techniques in elasticity).

So lets say that again you are working on the unit square and that you solve some 1D Problem to obtain $g$ as a 1D function on the left boundary. In that case I can think of two approaches.

Try to define your own Expression, see section 3 in this link. I think that you'd have to put something like

global g
value = g(x[1])

into the eval method. Unfortunately, I can't give much advice on this approach as the Expression class is somehow broken on my computer.

The second approach would be to solve a simple problem (such as the laplace equation) with g as boundary data to obtain a function that lives on all of $\Omega$.

Those are just very rough guesses but without further information on your problem I can't really help.

---

First of all, many thanks for the time you already took to answer my question.

I provided a small working example above, but in practice I have a complex N-D mesh on which I define several boundaries. On these boundaries, I solve a (N-1)-D eigenvalue problem, the solutions of which define a mixed bc for the eventual PDE. If the mesh were rectangular, I could define an Expression that only uses one coordinate, and is hence defined over the whole domain. For a general shape, I cannot however.

To require these eigenfunctions to be a Dirichlet boundary condition for a dummy problem and so obtain a function over the whole domain seemed like a good idea at first, but then I realised this leads me to the original problem: how to use a function that is only defined on the boundary as a boundary condition?

The (only?) solution seems to lie in defining a function in V, with an array filled with zeros everywhere, and then map the dof's of the boundary-functionspace (a submesh of a boundarymesh) to that of V. Perhaps I should have a more in depth look on how to do this (and if this is possible for polynomial degree >= 2).

---

As a small update, I looked into the last option of the reply above, as you suggested it was possible. I found a solution that I posted as an answer to this question:
http://fenicsproject.org/qa/6810/vertex-mapping-from-submesh-boundarymesh-back-actual-mesh?show=8110#a8110

It allows me to use functions defined on a boundary as boundary conditions (since as commented on the original question, the simple interpolate workaround stopped working in 1.6, which I also use).

Unfortunately, when trying this approach to Lagrange polynomials for degree 2, I got:

*** -------------------------------------------------------------------------
*** Error:   Unable to tabulate dof to vertex map.
*** Reason:  Can only tabulate dofs on vertices.
*** Where:   This error was encountered inside DofMap.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  
*** -------------------------------------------------------------------------

See also http://fenicsproject.org/qa/5917/dof_to_vertex_map-for-higher-order for this error. It seems that now not only the vertices need to be mapped (submesh --> boundarymesh --> mesh) but also other nodes, and I am not longer sure if this is already possible. I might try to figure it out later.

---

Another update (which answers the question this time for arbitrary degree): http://fenicsproject.org/qa/8522/mapping-degrees-of-freedom-between-related-meshes