solving on submesh: meaning of "allow_extrapolation"?

I don't think the problem here is that you are using a sub mesh. It's
just that for that particular mesh, you happen to get round-off errors
which prevent DOLFIN (or rather CGAL) to find some specific vertex or
quadrature point inside the mesh (it happens to be located eps > 0
outside of the domain).

Setting allow_extrapolation = true will allow evaluation at those
points even if they seem to be just outside of the domain. It should
be perfectly safe to set that option.

--
Anders

On Sat, Oct 06, 2012 at 11:41:09AM -0000, Nico Schlömer wrote:
> New question #210508 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/210508
>
> When solving the Poisson problem on a subdomain, i.e.,
>
> mesh = Mesh("coil.xml")
> subdomains = MeshFunction('uint', mesh, 'coil_physical_region.xml')
> submesh = SubMesh(mesh, subdomains, 3)
>
> all runs fine if the RHS function is defined on the entire mesh,
>
> rhs = project(Constant(1.0), FunctionSpace(mesh, "Lagrange", 1))
>
> -- Nice.
> The situation appears to be different when the right hand side is a function of Nedelec elements:
>
> T = project(Constant([0,0,1]), FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1))
> rhs = inner(T, T)
>
> gives
>
> Traceback (most recent call last):
> File "poisson.py", line 41, in <module>
> solve(a == L, u, bc)
> [...]
> *** Error: Unable to evaluate function at point.
> *** Reason: The point is not inside the domain. Consider setting "allow_extrapolation" to allow extrapolation.
> *** Where: This error was encountered inside Function.cpp.
> [...]
>
> The same error message appears when trying to interpolate T to the submesh,
>
> T = project(Constant([0,0,1]), FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1))
> T_submesh = interpolate(T, FunctionSpace(submesh, 'Nedelec 1st kind H(curl)', 1))
>
> Setting
>
> parameters['allow_extrapolation'] = True
>
> indeed works around the problem, but of course I'm asking myself what this actually does and why it is necessary.
>
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Oh, I was thinking with a MeshFunction provided in an XML file, no more float/double comparisons would be necessary. After all, the XML file specifies exactly what the regions should be.
Is that not correct?

On Sat, Oct 06, 2012 at 02:41:15PM -0000, Nico Schlömer wrote:
> Question #210508 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210508
>
> Nico Schlömer posted a new comment:
> Oh, I was thinking with a MeshFunction provided in an XML file, no more float/double comparisons would be necessary. After all, the XML file specifies exactly what the regions should be.
> Is that not correct?

It depends on what you do. If you project or integrate some expresion,
it will need to be interpolated into a finite element space which
leads to point evaluation at nodal points.

--
Anders

Oh I see. And since the Nedelec elements are defined on the edges, the values get interpolated onto the Gauss points to quadrature. The calculations of the Gauss points may then be slightly off, resulting in function evaluation outside of the actual domain (but only by machine precision).
Is this correct?

Perfect, thanks!

I got another instance of the same error message, cf.

=============== *snip* ===============
from dolfin import *

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5 + DOLFIN_EPS
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return 0.5 - DOLFIN_EPS < x[0]

left = Left()
right = Right()

mesh = UnitSquareMesh(20, 20)

domains = CellFunction('size_t', mesh)
domains.set_all(0)
left.mark(domains, 1)
right.mark(domains, 2)

f = Constant(1.0)

V = FunctionSpace(mesh, 'CG', 1)

bcs = DirichletBC(V, 0.0, 'on_boundary')

dx = Measure('dx')[domains]

K = {0: 1.0, 1: 0.1, 2: 1.0}

u = TrialFunction(V)
v = TestFunction(V)

a = zero() * dx(0)
L = zero() * dx(0)
for i in [0,1,2]:
    a += inner(K[i]*grad(u), grad(v)) * dx(i)
    L += f*v*dx(i)

# Set a convection in the left subdomain (index 0).
submesh = SubMesh(mesh, domains, 1)
VL = FunctionSpace(submesh, 'CG', 1)
c = project(Constant((1.0,2.0)), VL*VL)

# This works.
#c = Constant((1.0,2.0))

a += v * inner(grad(u), c) * dx(1)

assemble(a, bcs=bcs)
=============== *snap* ===============

This time it's certainly related to subdomains, more precisely the convection term `c` that's only defined on a subdomain. I'm guessing that something wants to evaluate `c` outside of `submesh`, but really I can't see why that would be.

This small example is harmless, but still this error message is making me really uncomfortable. In the larger and much more involved application code, allow_extrapolation might allow extrapolation at places where it isn't desired, e.g., where there's actually a bug.