Merge pull request #222 from firedrakeproject/wence/fix/compile-expre…

…ssion-domain

Add optional domain argument to expression compilation

tsfc/driver.py

@@ -267,7 +267,8 @@ def name_multiindex(multiindex, name):
     return builder.construct_kernel(kernel_name, impero_c, index_names, quad_rule)
 
 
-def compile_expression_dual_evaluation(expression, to_element, coordinates, interface=None,
+def compile_expression_dual_evaluation(expression, to_element, coordinates, *,
+                                       domain=None, interface=None,
                                        parameters=None, coffee=False):
     """Compile a UFL expression to be evaluated against a compile-time known reference element's dual basis.
 
@@ -276,6 +277,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, inte
     :arg expression: UFL expression
     :arg to_element: A FIAT FiniteElement for the target space
     :arg coordinates: the coordinate function
+    :arg domain: optional UFL domain the expression is defined on (useful when expression contains no domain).
     :arg interface: backend module for the kernel interface
     :arg parameters: parameters object
     :arg coffee: compile coffee kernel instead of loopy kernel
@@ -300,7 +302,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, inte
         mapping, = set(to_element.mapping())
     except ValueError:
         raise NotImplementedError("Don't know how to interpolate onto zany spaces, sorry")
-    expression = apply_mapping(expression, mapping)
+    expression = apply_mapping(expression, mapping, domain)
 
     # Apply UFL preprocessing
     expression = ufl_utils.preprocess_expression(expression,

tsfc/ufl_utils.py

@@ -338,7 +338,7 @@ def simplify_abs(expression, complex_mode):
     return mapper(expression, False)
 
 
-def apply_mapping(expression, mapping):
+def apply_mapping(expression, mapping, domain):
     """
     This applies the appropriate transformation to the
     given expression for interpolation to a specific
@@ -387,6 +387,10 @@ def apply_mapping(expression, mapping):
     """
 
     mesh = expression.ufl_domain()
+    if mesh is None:
+        mesh = domain
+    if domain is not None and mesh != domain:
+        raise NotImplementedError("Multiple domains not supported")
     rank = len(expression.ufl_shape)
     if mapping == "affine":
         return expression
@@ -396,18 +400,15 @@ def apply_mapping(expression, mapping):
         expression = Indexed(expression, MultiIndex((*i, k)))
         return as_tensor(J.T[j, k] * expression, (*i, j))
     elif mapping == "contravariant piola":
-        mesh = expression.ufl_domain()
         K = JacobianInverse(mesh)
         detJ = JacobianDeterminant(mesh)
         *i, j, k = indices(len(expression.ufl_shape) + 1)
         expression = Indexed(expression, MultiIndex((*i, k)))
         return as_tensor(detJ * K[j, k] * expression, (*i, j))
     elif mapping == "double covariant piola" and rank == 2:
-        mesh = expression.ufl_domain()
         J = Jacobian(mesh)
         return J.T * expression * J
     elif mapping == "double contravariant piola" and rank == 2:
-        mesh = expression.ufl_domain()
         K = JacobianInverse(mesh)
         detJ = JacobianDeterminant(mesh)
         return (detJ)**2 * K * expression * K.T