Representation fails with coefficient from different mesh

[Jan Blechta]

The following example works with quadrature but fails with uflacs

from dolfin import *
parameters['form_compiler']['representation'] = 'uflacs'

# Create mesh and submesh
mesh = UnitSquareMesh(4, 4)
cf = CellFunction('size_t', mesh)
cf[0] = 1
submesh = SubMesh(mesh, cf, 1)

V = FunctionSpace(mesh, 'Lagrange', 1)
V_loc = FunctionSpace(submesh, 'Lagrange', 4)

# Prepare functional on V_loc with coefficient from V
u = Function(V)
#R = derivative(0.5*u**2*dx(submesh), u, TestFunction(V_loc)) # works
#R = derivative(0.5*grad(u)**2*dx(submesh), u, TestFunction(V_loc)) # fails
#R = u*TestFunction(V_loc)*dx(submesh) # works
R = inner(grad(u), grad(TestFunction(V_loc)))*dx(submesh) # fails

# Solve problem on V_loc with R
u_loc = Function(V_loc)
a_loc = inner(grad(TrialFunction(V_loc)), grad(TestFunction(V_loc)))*dx(submesh)
bc_loc = DirichletBC(V_loc, 0.0, lambda x, onb: onb)
solve(a_loc == R, u_loc, bc_loc)
plot(u_loc, interactive=True)

The message is

  File "uflacs/elementtables/table_utils.py", line 125, in get_ffc_table_values
    element_table = tables[num_points][element]
KeyError: VectorElement('Lagrange', Domain(Cell('triangle', 2), label='dolfin_mesh_with_id_0', data='<data with id 0>'), 1, dim=2, quad_scheme=None)

[Jan Blechta]

Simplify MWE a bit.

[Jan Blechta]

Apperently, the problem is in mixing meshes. When one does submesh = mesh instead of submesh = SubMesh(...) the program passes.

[Jan Blechta]

Apparently, the problem is in dtables. With mass term R = u*TestFunction(V_loc)*dx(submesh) works while Laplacian R = inner(grad(u), grad(TestFunction(V_loc)))*dx(submesh) fails.

[Martin Sandve Alnæs]

This test program still fails with the new updates pushed to master today.

But I have to ask: is this supposed to work?

AFAIK, what actually happens behind the scene in the 'quadrature' case is that ffc just ignores the difference in domain/mesh, compiling as if there was just one mesh, resulting in a ufc kernel that assumes a representation of u on the same cell as the cell from the other mesh. Dolfin compensates for this by silently interpolating u onto the integration mesh, resulting in a hidden interpolation error.

[Martin Sandve Alnæs]

@garth-wells @logg @blechta @meg opinions?

[Martin Sandve Alnæs]

I guess I'll have to fix it for now, but I'm not sure of the long term status of this 'feature'. This is closely related to the silent interpolation of Expressions, which we've more or less agreed is not desireable. A more consistent approach would be to evaluate u at quadrature points instead of interpolating in the case of nonmatching meshes, and for ufl to know about the submesh relation with the planned meshview features, avoiding costly interpolation in that case.

[Marie Elisabeth Rognes]

I was not aware that we supported this, but in this example it is at least mathematically well-defined.

[Jan Blechta]

resulting in a hidden interpolation error.

Anyway even for quadrature approach you suggest, user should expect there is some kind of approximation error in the term R = inner(grad(u), grad(TestFunction(V_loc)))*dx(submesh) if mesh and submesh do not match. grad(u) is not polymomial on cells of submesh so one does not have exact quadrature (unless cut cells are introduced and integrated separately).

[Martin Sandve Alnæs]

True. Is the error the same now? Domains are gone from elements now. It might work or fail elsewhere.

[Jan Blechta]

Now it fails in C++ compilation with FEniCS master

/tmp/tmpdm66Zt2015-10-29-18-17_instant_c644cdf006a0907ddf4acf30a9bcbb85415221f4/ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d/ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d.h: In member function ‘virtual void ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d_cell_integral_0_otherwise::tabulate_tensor(double*, const double* const*, const double*, int) const’:
/tmp/tmpdm66Zt2015-10-29-18-17_instant_c644cdf006a0907ddf4acf30a9bcbb85415221f4/ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d/ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d.h:3204:18: error: redeclaration of ‘const double J_c0’
     const double J_c0 = coordinate_dofs[0] * FE0_C0_D10_Q6[0][0][0 - 0] + coordinate_dofs[2] * FE0_C0_D10_Q6[0][0][1 - 0];
                  ^
/tmp/tmpdm66Zt2015-10-29-18-17_instant_c644cdf006a0907ddf4acf30a9bcbb85415221f4/ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d/ffc_form_6833d9bd1c855c35d0be4046a3ca48c86dee447d.h:3200:18: note: ‘const double J_c0’ previously declared here
     const double J_c0 = coordinate_dofs[0] * FE0_C0_D10_Q6[0][0][0 - 0] + coordinate_dofs[2] * FE0_C0_D10_Q6[0][0][1 - 0];
                  ^

[Jan Blechta]

It's in

  void tabulate_tensor(double * A,
                       const double * const * w,
                       const double * coordinate_dofs,
                       int cell_orientation) const final override
  {
...
    // Section for piecewise constant computations
    const double J_c0 = coordinate_dofs[0] * FE0_C0_D10_Q6[0][0][0 - 0] + coordinate_dofs[2] * FE0_C0_D10_Q6[0][0][1 - 0];
...
    const double J_c0 = coordinate_dofs[0] * FE0_C0_D10_Q6[0][0][0 - 0] + coordinate_dofs[2] * FE0_C0_D10_Q6[0][0][1 - 0];

[Martin Sandve Alnæs]

Ok, that's because the Jacobian is defined in terms of the domain.

One way to "fix" this for now can be to throw away information about which domain things are defined on. That will work for legacy cases but will block new work where domain identity plays a role. At least that will make the reinterpretation explicit in the code.

I really don't like this feature... Alternatives are to ask the user to call interpolate explicitly instead, or to add an Interpolant(Function(U), V) usable in forms to make it explicit that way. I don't know what the consequences of that design choice will be in other code.

[Jan Blechta]

I really don't like this feature...

I don't have a strong opinion here. But note that an approximate behaviour is already present in code generation for forms, e.g. non-polynomial [e|E]xpressions. It's always a compromise between simplicity and precision of notation.

Interpolant(Function(U), V)

Is it necessary to interpolate to V instance of dolfin.FunctionSpace which needs dofmap to be built? Could not interpolation happen at element level in a generated code once ufl.FunctionSpace V has nodal basis so it's capable of cheap element-wise interpolation?

[Martin Sandve Alnæs]

I haven't explored concequences of design decisions, this was just a possible solution. Short term I'll try to make a quick fix to keep current behaviour.

And we've already discussed Expression and there's general agreement that it's not a good design and should change.

[Jan Blechta]

For a sake of argument, just ignore dolfin.Expression for a while. The point is that UFL expression gets approximate quadrature (at least) when

1. it's on a same mesh but not polynomial
2. it's on a different mesh hence not polynomial

It would be inconsistent if 1. would not require explicit Interpolant request (as it is now) while we would require it from 2.

EDIT: I'm not saying I'm against Interpolant. I think it could be nice having it as it would shed a more light on mathematics behind the UFL code and it's generated implementation.

[Martin Sandve Alnæs]

Following that line of thought I would say it's more consistent to evaluate a function from another mesh in the quadrature points. A non-polynomial function on the same mesh is not interpolated into a polynomial space.

[Martin Sandve Alnæs]

Argh. I'm not sure where to make the fix. The jacobian generation code in uflacs is supposed to be(come) capable of handling geometry from multiple domains, it basically just doesn't add a domain number to J yet.

I think the most consistent way will be to fix it in ufl compute_form_data next to element mappings.

Is this a rule that captures the current behaviour properly?

If an integral over integration_mesh has type cell/interor_facet/exterior_facet but not the new multimesh types:
    replace any
        Coefficient(FunctionSpace(other_mesh, element))
    with
        Coefficient(FunctionSpace(integration_mesh, element))
    and any geometry on other_mesh with the same geometry on integration_mesh

For multimesh there is no unique integration_mesh to assume the function is defined over. Is this correct @logg ?

[Jan Blechta]

@martinal, is there a good place to file somewhere a general problem with hidden interpolation (or is already filed)?

BTW, any idea why referencing @martinal does not work for you? Or at least it doesn't get rendered.

[Martin Sandve Alnæs]

@blechta that's a dolfin issue really, ufl and ffc doesn't know about that.

It's related to this: https://bitbucket.org/fenics-project/dolfin/issues/355/expression-s-interpolate-to-linear

Testing @martinal ... weird.

[Jan Blechta]

ufl and ffc doesn't know about that

The question is whether UFL should know about that. It would either raise an error or allow to generate some code if explicitly requested to do so (using ufl.NodalInterpolant).

from dolfin import *
import numpy as np

mesh1 = UnitIntervalMesh(1)
mesh2 = UnitIntervalMesh(2)
V2 = FunctionSpace(mesh2, 'Lagrange', 1)
u2 = Function(V2)
u2.vector()[1] = 1.0

# Works fine
assert np.isclose(assemble(u2*dx(mesh2)), 0.5)

# Fails; should raise an error
assert np.isclose(assemble(u2*dx(mesh1)), 0.5)

# Request nodal interpolation allowing execution of
#  GenericFunction::restrict_as_ufc_function()
#  finite_element::evaluate_dofs()
assemble(ufl.NodalInterpolant(u2, mesh1)*dx(mesh1)) # Will return zero

# Or interpolation to different element (which should not happen?!)
assemble(ufl.NodalInterpolant(u2, ufl_function_space)*dx(mesh1))
assemble(ufl_function_space.nodal_interpolant(u2)*dx(mesh1))