diff --git a/tests/test_dual_evaluation.py b/tests/test_dual_evaluation.py new file mode 100644 index 00000000..8cabba70 --- /dev/null +++ b/tests/test_dual_evaluation.py @@ -0,0 +1,61 @@ +import pytest +import ufl +from tsfc.finatinterface import create_element +from tsfc import compile_expression_dual_evaluation + + +def test_ufl_only_simple(): + mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1)) + V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2)) + v = ufl.Coefficient(V) + expr = ufl.inner(v, v) + W = V + to_element = create_element(W.ufl_element()) + ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False) + assert first_coeff_fake_coords is False + + +def test_ufl_only_spatialcoordinate(): + mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1)) + V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2)) + x, y = ufl.SpatialCoordinate(mesh) + expr = x*y - y**2 + x + W = V + to_element = create_element(W.ufl_element()) + ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False) + assert first_coeff_fake_coords is True + + +def test_ufl_only_from_contravariant_piola(): + mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1)) + V = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1)) + v = ufl.Coefficient(V) + expr = ufl.inner(v, v) + W = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2)) + to_element = create_element(W.ufl_element()) + ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False) + assert first_coeff_fake_coords is True + + +def test_ufl_only_to_contravariant_piola(): + mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1)) + V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2)) + v = ufl.Coefficient(V) + expr = ufl.as_vector([v, v]) + W = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1)) + to_element = create_element(W.ufl_element()) + ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False) + assert first_coeff_fake_coords is True + + +def test_ufl_only_shape_mismatch(): + mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1)) + V = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1)) + v = ufl.Coefficient(V) + expr = ufl.inner(v, v) + assert expr.ufl_shape == () + W = V + to_element = create_element(W.ufl_element()) + assert to_element.value_shape == (2,) + with pytest.raises(ValueError): + ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False) diff --git a/tsfc/driver.py b/tsfc/driver.py index b8160a85..d6db7120 100644 --- a/tsfc/driver.py +++ b/tsfc/driver.py @@ -267,7 +267,7 @@ 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, *, +def compile_expression_dual_evaluation(expression, to_element, *, 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,8 +276,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *, :arg expression: UFL expression :arg to_element: A FInAT element 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 domain: optional UFL domain the expression is defined on (required 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 @@ -328,17 +327,21 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *, argument_multiindices = tuple(builder.create_element(arg.ufl_element()).get_indices() for arg in arguments) - # Replace coordinates (if any) - domain = expression.ufl_domain() - if domain: - assert coordinates.ufl_domain() == domain - builder.domain_coordinate[domain] = coordinates - builder.set_cell_sizes(domain) + # Replace coordinates (if any) unless otherwise specified by kwarg + if domain is None: + domain = expression.ufl_domain() + assert domain is not None # Collect required coefficients + first_coefficient_fake_coords = False coefficients = extract_coefficients(expression) if has_type(expression, GeometricQuantity) or any(fem.needs_coordinate_mapping(c.ufl_element()) for c in coefficients): - coefficients = [coordinates] + coefficients + # Create a fake coordinate coefficient for a domain. + coords_coefficient = ufl.Coefficient(ufl.FunctionSpace(domain, domain.ufl_coordinate_element())) + builder.domain_coordinate[domain] = coords_coefficient + builder.set_cell_sizes(domain) + coefficients = [coords_coefficient] + coefficients + first_coefficient_fake_coords = True builder.set_coefficients(coefficients) # Split mixed coefficients @@ -346,7 +349,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *, # Translate to GEM kernel_cfg = dict(interface=builder, - ufl_cell=coordinates.ufl_domain().ufl_cell(), + ufl_cell=domain.ufl_cell(), argument_multiindices=argument_multiindices, index_cache={}, scalar_type=parameters["scalar_type"]) @@ -431,7 +434,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *, # Handle kernel interface requirements builder.register_requirements([ir]) # Build kernel tuple - return builder.construct_kernel(return_arg, impero_c, index_names) + return builder.construct_kernel(return_arg, impero_c, index_names, first_coefficient_fake_coords) def lower_integral_type(fiat_cell, integral_type): diff --git a/tsfc/kernel_interface/firedrake_loopy.py b/tsfc/kernel_interface/firedrake_loopy.py index ee5cafa4..eb90cf72 100644 --- a/tsfc/kernel_interface/firedrake_loopy.py +++ b/tsfc/kernel_interface/firedrake_loopy.py @@ -17,7 +17,7 @@ # Expression kernel description type -ExpressionKernel = namedtuple('ExpressionKernel', ['ast', 'oriented', 'needs_cell_sizes', 'coefficients', 'tabulations']) +ExpressionKernel = namedtuple('ExpressionKernel', ['ast', 'oriented', 'needs_cell_sizes', 'coefficients', 'first_coefficient_fake_coords', 'tabulations']) def make_builder(*args, **kwargs): @@ -153,12 +153,14 @@ def register_requirements(self, ir): provided by the kernel interface.""" self.oriented, self.cell_sizes, self.tabulations = check_requirements(ir) - def construct_kernel(self, return_arg, impero_c, index_names): + def construct_kernel(self, return_arg, impero_c, index_names, first_coefficient_fake_coords): """Constructs an :class:`ExpressionKernel`. :arg return_arg: loopy.GlobalArg for the return value :arg impero_c: gem.ImperoC object that represents the kernel :arg index_names: pre-assigned index names + :arg first_coefficient_fake_coords: If true, the kernel's first + coefficient is a constructed UFL coordinate field :returns: :class:`ExpressionKernel` object """ args = [return_arg] @@ -173,7 +175,8 @@ def construct_kernel(self, return_arg, impero_c, index_names): loopy_kernel = generate_loopy(impero_c, args, self.scalar_type, "expression_kernel", index_names) return ExpressionKernel(loopy_kernel, self.oriented, self.cell_sizes, - self.coefficients, self.tabulations) + self.coefficients, first_coefficient_fake_coords, + self.tabulations) class KernelBuilder(KernelBuilderBase):