Expunge ufl_domain() on non Expr

docs/source/interpolation.rst

@@ -251,7 +251,7 @@ proceeds as follows:
 .. code-block:: python3
 
    # First, grab the mesh.
-   m = V.ufl_domain()
+   m = V.mesh()
 
    # Now make the VectorFunctionSpace corresponding to V.
    W = VectorFunctionSpace(m, V.ufl_element())

firedrake/assemble.py

@@ -1338,7 +1338,7 @@ def _as_global_kernel_arg_cell_facet(_, self):
 @_as_global_kernel_arg.register(kernel_args.CellOrientationsKernelArg)
 def _as_global_kernel_arg_cell_orientations(_, self):
     # this mirrors firedrake.mesh.MeshGeometry.init_cell_orientations
-    ufl_element = finat.ufl.FiniteElement("DG", cell=self._form.ufl_domain().ufl_cell(), degree=0)
+    ufl_element = finat.ufl.FiniteElement("DG", cell=self._mesh.ufl_cell(), degree=0)
     finat_element = create_element(ufl_element)
     return self._make_dat_global_kernel_arg(finat_element)
 

firedrake/functionspaceimpl.py

@@ -118,7 +118,7 @@ def ufl_function_space(self):
 
     def ufl_cell(self):
         r"""The :class:`~ufl.classes.Cell` this FunctionSpace is defined on."""
-        return self.ufl_domain().ufl_cell()
+        return self.mesh().ufl_cell()
 
     @PETSc.Log.EventDecorator()
     def split(self):

firedrake/interpolation.py

@@ -852,7 +852,7 @@ def _interpolator(V, tensor, expr, subset, arguments, access, bcs=None):
                            % (expr.ufl_shape, V.ufl_element().value_shape))
 
     # NOTE: The par_loop is always over the target mesh cells.
-    target_mesh = V.ufl_domain()
+    target_mesh = as_domain(V)
     source_mesh = extract_unique_domain(expr) or target_mesh
 
     if target_mesh is not source_mesh:

firedrake/output.py

@@ -73,7 +73,7 @@ def is_cg(V):
 
     :arg V: A FunctionSpace.
     """
-    nvertex = V.ufl_domain().ufl_cell().num_vertices()
+    nvertex = V.mesh().ufl_cell().num_vertices()
     entity_dofs = V.finat_element.entity_dofs()
     # If there are as many dofs on vertices as there are vertices,
     # assume a continuous space.
@@ -96,7 +96,7 @@ def is_linear(V):
 
     :arg V: A FunctionSpace.
     """
-    nvertex = V.ufl_domain().ufl_cell().num_vertices()
+    nvertex = V.mesh().ufl_cell().num_vertices()
     return V.finat_element.space_dimension() == nvertex
 
 
@@ -135,7 +135,7 @@ def get_topology(coordinates):
     V = coordinates.function_space()
 
     nonLinear = not is_linear(V)
-    mesh = V.ufl_domain().topology
+    mesh = V.mesh().topology
     cell = mesh.ufl_cell()
     values = V.cell_node_map().values
     value_shape = values.shape

firedrake/plot.py

@@ -1000,7 +1000,7 @@ def _pgfplot_create_patches(f, coords, complex_component):
     fiat_cell = V.finat_element.cell
     degree = elem.degree()
     coordV = coords.function_space()
-    mesh = V.ufl_domain()
+    mesh = V.mesh()
     cdata = coords.dat.data_ro.real
     fdata = f.dat.data_ro.real if complex_component == 'real' else f.dat.data_ro.imag
     map_facet_dofs, patch_type = _pgfplot_make_perms(fiat_cell, degree)
@@ -1068,7 +1068,7 @@ def pgfplot(f, filename, degree=1, complex_component='real', print_latex_example
         raise NotImplementedError(f"complex_component must be {'real', 'imag'}: got {complex_component}")
     V = f.function_space()
     elem = V.ufl_element()
-    mesh = V.ufl_domain()
+    mesh = V.mesh()
     dim = mesh.geometric_dimension()
     if dim not in (2, 3):
         raise NotImplementedError(f"Not yet implemented for functions in spatial dimension {dim}")

firedrake/preconditioners/patch.py

@@ -192,11 +192,11 @@ def matrix_funptr(form, state):
         arg = mesh.coordinates.dat(op2.READ, get_map(mesh.coordinates))
         args.append(arg)
         if kinfo.oriented:
-            c = form.ufl_domain().cell_orientations()
+            c = mesh.cell_orientations()
             arg = c.dat(op2.READ, get_map(c))
             args.append(arg)
         if kinfo.needs_cell_sizes:
-            c = form.ufl_domain().cell_sizes
+            c = mesh.cell_sizes
             arg = c.dat(op2.READ, get_map(c))
             args.append(arg)
         for n, indices in kinfo.coefficient_numbers:
@@ -217,7 +217,7 @@ def matrix_funptr(form, state):
             args.append(all_constants[constant_index].dat(op2.READ))
 
         if kinfo.integral_type == "interior_facet":
-            arg = test.ufl_domain().interior_facets.local_facet_dat(op2.READ)
+            arg = mesh.interior_facets.local_facet_dat(op2.READ)
             args.append(arg)
         iterset = op2.Subset(iterset, [])
 
@@ -286,11 +286,11 @@ def residual_funptr(form, state):
         args.append(arg)
 
         if kinfo.oriented:
-            c = form.ufl_domain().cell_orientations()
+            c = mesh.cell_orientations()
             arg = c.dat(op2.READ, get_map(c))
             args.append(arg)
         if kinfo.needs_cell_sizes:
-            c = form.ufl_domain().cell_sizes
+            c = mesh.cell_sizes
             arg = c.dat(op2.READ, get_map(c))
             args.append(arg)
         for n, indices in kinfo.coefficient_numbers:
@@ -514,11 +514,12 @@ def load_c_function(code, name, comm):
 
 def make_c_arguments(form, kernel, state, get_map, require_state=False,
                      require_facet_number=False):
-    coeffs = [form.ufl_domain().coordinates]
+    mesh = form.ufl_domains()[kernel.kinfo.domain_number]
+    coeffs = [mesh.coordinates]
     if kernel.kinfo.oriented:
-        coeffs.append(form.ufl_domain().cell_orientations())
+        coeffs.append(mesh.cell_orientations())
     if kernel.kinfo.needs_cell_sizes:
-        coeffs.append(form.ufl_domain().cell_sizes)
+        coeffs.append(mesh.cell_sizes)
     for n, indices in kernel.kinfo.coefficient_numbers:
         c = form.coefficients()[n]
         if c is state:
@@ -550,7 +551,7 @@ def make_c_arguments(form, kernel, state, get_map, require_state=False,
         data_args.extend(all_constants[constant_index].dat._kernel_args_)
 
     if require_facet_number:
-        data_args.extend(form.ufl_domain().interior_facets.local_facet_dat._kernel_args_)
+        data_args.extend(mesh.interior_facets.local_facet_dat._kernel_args_)
     return data_args, map_args
 
 
@@ -778,7 +779,8 @@ def initialize(self, obj):
             J = ctx.Jp or ctx.J
             bcs = ctx._problem.bcs
 
-        mesh = J.ufl_domain()
+        V, _ = map(operator.methodcaller("function_space"), J.arguments())
+        mesh = V.mesh()
         self.plex = mesh.topology_dm
         # We need to attach the mesh and appctx to the plex, so that
         # PlaneSmoothers (and any other user-customised patch
@@ -810,8 +812,6 @@ def initialize(self, obj):
             Jstate = None
             is_snes = False
 
-        V, _ = map(operator.methodcaller("function_space"), J.arguments())
-
         if len(bcs) > 0:
             ghost_bc_nodes = numpy.unique(numpy.concatenate([bcdofs(bc, ghost=True)
                                                              for bc in bcs]))
@@ -838,7 +838,7 @@ def initialize(self, obj):
                                                                        require_facet_number=True)
             code, Struct = make_jacobian_wrapper(facet_Jop_data_args, facet_Jop_map_args)
             facet_Jop_function = load_c_function(code, "ComputeJacobian", obj.comm)
-            point2facet = J.ufl_domain().interior_facets.point2facetnumber.ctypes.data
+            point2facet = mesh.interior_facets.point2facetnumber.ctypes.data
             facet_Jop_struct = make_c_struct(facet_Jop_data_args, facet_Jop_map_args,
                                              Jint_facet_kernel.funptr, Struct,
                                              point2facet=point2facet)

tests/extrusion/test_real_tensorproduct.py

@@ -28,7 +28,7 @@ def variant(request):
 
 @pytest.fixture
 def expr(variant, V, fs_kind):
-    x, y, z = SpatialCoordinate(V.ufl_domain())
+    x, y, z = SpatialCoordinate(V.mesh())
     val = {"linear": z, "sin": sin(pi*z)}[variant]
     ret = {
         "scalar": val,

tests/multigrid/test_non_nested.py

@@ -25,7 +25,7 @@ def coarsen(expr, self, coefficient_mapping=None):
     @coarsen.register(functionspaceimpl.FunctionSpace)
     @coarsen.register(functionspaceimpl.WithGeometry)
     def coarsen_fs(V, self, coefficient_mapping=None):
-        mesh = self(V.ufl_domain(), self)
+        mesh = self(V.mesh(), self)
         return FunctionSpace(mesh, "CG", 1)
 
     uh = Function(V)

tests/multigrid/test_transfer_manager.py

@@ -49,7 +49,7 @@ def test_transfer_manager_inside_coarsen(sub, mesh):
 
     bc, = cctx._problem.bcs
     V = bc.function_space()
-    mesh = V.ufl_domain()
+    mesh = V.mesh()
     x, y = SpatialCoordinate(mesh)
     expect = project(as_vector([-y, x]), V)
     assert numpy.allclose(bc.function_arg.dat.data_ro, expect.dat.data_ro)

tests/regression/test_auxiliary_dm.py

@@ -28,7 +28,7 @@ def mesh(self):
         return mesh
 
     def bcs(self, Z):
-        bcs = [DirichletBC(Z.sub(0), self.analytical_solution(Z.ufl_domain()), "on_boundary"),
+        bcs = [DirichletBC(Z.sub(0), self.analytical_solution(Z.mesh()), "on_boundary"),
                DirichletBC(Z.sub(1), zero(), "on_boundary")]
         return bcs
 

tests/regression/test_helmholtz_bernstein.py

@@ -44,7 +44,7 @@ def helmholtz(V):
     u = TrialFunction(V)
     v = TestFunction(V)
     f = Function(V)
-    x = SpatialCoordinate(V.ufl_domain())
+    x = SpatialCoordinate(V.mesh())
     f.project(np.prod([cos(2*pi*xi) for xi in x]))
     a = (inner(grad(u), grad(v)) + inner(u, v)) * dx
     L = inner(f, v) * dx

tests/regression/test_interpolate_vs_project.py

@@ -29,7 +29,7 @@ def V(request, mesh):
 
 
 def test_interpolate_vs_project(V):
-    mesh = V.ufl_domain()
+    mesh = V.mesh()
     dim = mesh.geometric_dimension()
     if dim == 2:
         x, y = SpatialCoordinate(mesh)

tests/regression/test_interpolate_zany.py

@@ -41,7 +41,7 @@ def tolerance(element, which):
 
 @pytest.fixture
 def expect(V, which):
-    x, y = SpatialCoordinate(V.ufl_domain())
+    x, y = SpatialCoordinate(V.mesh())
     expr = (x + y)**(V.ufl_element().degree())
     if which == "coefficient":
         return expr

tests/regression/test_interpolation_nodes.py

@@ -39,7 +39,7 @@ def V(request, mesh, degree):
 
 
 def test_div_curl_preserving(V):
-    mesh = V.ufl_domain()
+    mesh = V.mesh()
     dim = mesh.geometric_dimension()
     if dim == 3 and V.ufl_element().degree() == 3 and "Nedelec" not in V.ufl_element().family():
         pytest.skip("N2div interpolation kernel with exact quadrature creates tensors which risk stack overflow")