From bb5a9bd020d03907ceb82033359cfeebe4d2b6f9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 13 Nov 2024 15:56:58 +0000 Subject: [PATCH] fix tests --- firedrake/formmanipulation.py | 2 +- firedrake/function.py | 2 +- firedrake/interpolation.py | 2 +- firedrake/mesh.py | 2 +- firedrake/pointquery_utils.py | 40 ++++++++++++++--------- tests/output/test_io_backward_compat.py | 5 ++- tests/regression/test_2dcohomology.py | 8 ++--- tests/submesh/test_submesh_interpolate.py | 6 ++-- 8 files changed, 38 insertions(+), 29 deletions(-) diff --git a/firedrake/formmanipulation.py b/firedrake/formmanipulation.py index 394fdcca95..35a6789107 100644 --- a/firedrake/formmanipulation.py +++ b/firedrake/formmanipulation.py @@ -124,7 +124,7 @@ def argument(self, o): args += [a_[j] for j in numpy.ndindex(a_.ufl_shape)] else: args += [Zero() - for j in numpy.ndindex(V_is[i].block_size)] + for j in numpy.ndindex(V_is[i].value_shape)] return self._arg_cache.setdefault(o, as_vector(args)) diff --git a/firedrake/function.py b/firedrake/function.py index ddedebad10..83136407ea 100644 --- a/firedrake/function.py +++ b/firedrake/function.py @@ -633,7 +633,7 @@ def at(self, arg, *args, **kwargs): raise NotImplementedError("Point evaluation not implemented for variable layers") # Validate geometric dimension - gdim = mesh.ufl_cell().geometric_dimension() + gdim = mesh.geometric_dimension() if arg.shape[-1] == gdim: pass elif len(arg.shape) == 1 and gdim == 1: diff --git a/firedrake/interpolation.py b/firedrake/interpolation.py index 2d1c6ef051..32d80be3a3 100644 --- a/firedrake/interpolation.py +++ b/firedrake/interpolation.py @@ -613,7 +613,7 @@ def __init__( # I first point evaluate my expression at these locations, giving a # P0DG function on the VOM. As described in the manual, this is an # interpolation operation. - shape = V_dest.value_shape + shape = V_dest.ufl_function_space().value_shape if len(shape) == 0: fs_type = firedrake.FunctionSpace elif len(shape) == 1: diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 8624240421..fd37e7b575 100644 --- a/firedrake/mesh.py +++ b/firedrake/mesh.py @@ -2920,7 +2920,7 @@ def make_vom_from_vom_topology(topology, name, tolerance=0.5): gdim = topology.topology_dm.getCoordinateDim() cell = topology.ufl_cell() element = finat.ufl.VectorElement("DG", cell, 0, dim=gdim) - vmesh = MeshGeometry.__new__(MeshGeometry, element) + vmesh = MeshGeometry.__new__(MeshGeometry, element, topology.comm) vmesh._init_topology(topology) # Save vertex reference coordinate (within reference cell) in function parent_tdim = topology._parent_mesh.ufl_cell().topological_dimension() diff --git a/firedrake/pointquery_utils.py b/firedrake/pointquery_utils.py index 238fcf4fd1..32df5b4454 100644 --- a/firedrake/pointquery_utils.py +++ b/firedrake/pointquery_utils.py @@ -34,7 +34,7 @@ def make_wrapper(function, **kwargs): def src_locate_cell(mesh, tolerance=None): src = ['#include '] - src.append(compile_coordinate_element(mesh.ufl_coordinate_element(), tolerance)) + src.append(compile_coordinate_element(mesh, tolerance)) src.append(make_wrapper(mesh.coordinates, forward_args=["void*", "double*", RealType_c+"*"], kernel_name="to_reference_coords_kernel", @@ -129,9 +129,10 @@ def to_reference_coords_newton_step(ufl_coordinate_element, parameters, x0_dtype # Set up UFL form cell = ufl_coordinate_element.cell domain = ufl.Mesh(ufl_coordinate_element) + gdim = domain.geometric_dimension() K = ufl.JacobianInverse(domain) x = ufl.SpatialCoordinate(domain) - x0_element = finat.ufl.VectorElement("Real", cell, 0) + x0_element = finat.ufl.VectorElement("Real", cell, 0, dim=gdim) x0 = ufl.Coefficient(ufl.FunctionSpace(domain, x0_element)) expr = ufl.dot(K, x - x0) @@ -192,11 +193,21 @@ def predicate(index): @PETSc.Log.EventDecorator() -def compile_coordinate_element(ufl_coordinate_element, contains_eps, parameters=None): +def compile_coordinate_element(mesh, contains_eps, parameters=None): """Generates C code for changing to reference coordinates. - :arg ufl_coordinate_element: UFL element of the coordinates - :returns: C code as string + Parameters + ---------- + mesh : MeshTopology + The mesh. + + contains_eps : float + The tolerance used to verify that a point is contained by a cell. + + Returns + ------- + str + A string of C code. """ if parameters is None: parameters = tsfc.default_parameters() @@ -204,25 +215,24 @@ def compile_coordinate_element(ufl_coordinate_element, contains_eps, parameters= _ = tsfc.default_parameters() _.update(parameters) parameters = _ + + ufl_coordinate_element = mesh.ufl_coordinate_element() # Create FInAT element element = tsfc.finatinterface.create_element(ufl_coordinate_element) - gdim, = ufl_coordinate_element.reference_value_shape - cell = ufl_coordinate_element.cell - extruded = isinstance(cell, ufl.TensorProductCell) code = { - "geometric_dimension": gdim, - "topological_dimension": cell.topological_dimension(), + "geometric_dimension": mesh.geometric_dimension(), + "topological_dimension": mesh.topological_dimension(), "celldist_l1_c_expr": celldist_l1_c_expr(element.cell, "X"), "to_reference_coords_newton_step": to_reference_coords_newton_step(ufl_coordinate_element, parameters), "init_X": init_X(element.cell, parameters), "max_iteration_count": 1 if is_affine(ufl_coordinate_element) else 16, "convergence_epsilon": 1e-12, - "dX_norm_square": dX_norm_square(cell.topological_dimension()), - "X_isub_dX": X_isub_dX(cell.topological_dimension()), - "extruded_arg": ", int const *__restrict__ layers" if extruded else "", - "extr_comment_out": "//" if extruded else "", - "non_extr_comment_out": "//" if not extruded else "", + "dX_norm_square": dX_norm_square(mesh.topological_dimension()), + "X_isub_dX": X_isub_dX(mesh.topological_dimension()), + "extruded_arg": ", int const *__restrict__ layers" if mesh.extruded else "", + "extr_comment_out": "//" if mesh.extruded else "", + "non_extr_comment_out": "//" if not mesh.extruded else "", "IntType": as_cstr(IntType), "ScalarType": ScalarType_c, "RealType": RealType_c, diff --git a/tests/output/test_io_backward_compat.py b/tests/output/test_io_backward_compat.py index 67cf0d6dff..0e43326a2a 100644 --- a/tests/output/test_io_backward_compat.py +++ b/tests/output/test_io_backward_compat.py @@ -150,7 +150,7 @@ def _get_mesh_and_V(params): def _get_expr(V): mesh = V.mesh() dim = mesh.geometric_dimension() - shape = V.ufl_element().value_shape + shape = V.value_shape if dim == 2: x, y = SpatialCoordinate(mesh) z = x + y @@ -240,9 +240,8 @@ def test_io_backward_compat_base_load(version, params): def _get_expr_timestepping(V, i): mesh = V.mesh() - element = V.ufl_element() x, y = SpatialCoordinate(mesh) - shape = element.value_shape + shape = V.value_shape if shape == (4, ): return as_vector([x + i, y + i, x * y + i, i * i]) else: diff --git a/tests/regression/test_2dcohomology.py b/tests/regression/test_2dcohomology.py index 73bf409538..85ad069924 100644 --- a/tests/regression/test_2dcohomology.py +++ b/tests/regression/test_2dcohomology.py @@ -51,7 +51,7 @@ def test_betti0(space, mesh): L = assemble(inner(nabla_grad(u), nabla_grad(v))*dx) - bc0 = DirichletBC(V0, Constant(0.0), 9) + bc0 = DirichletBC(V0, 0, 9) L0 = assemble(inner(nabla_grad(u), nabla_grad(v))*dx, bcs=[bc0]) u, s, v = linalg.svd(L.M.values) @@ -94,8 +94,8 @@ def test_betti1(space, mesh): L = assemble((inner(sigma, tau) - inner(u, rot(tau)) + inner(rot(sigma), v) + inner(div(u), div(v))) * dx) - bc0 = DirichletBC(W.sub(0), Constant(0.0), 9) - bc1 = DirichletBC(W.sub(1), Constant((0.0, 0.0)), 9) + bc0 = DirichletBC(W.sub(0), 0, 9) + bc1 = DirichletBC(W.sub(1), 0, 9) L0 = assemble((inner(sigma, tau) - inner(u, rot(tau)) + inner(rot(sigma), v) + inner(div(u), div(v))) * dx, bcs=[bc0, bc1]) @@ -155,7 +155,7 @@ def test_betti2(space, mesh): L = assemble((inner(sigma, tau) - inner(u, div(tau)) + inner(div(sigma), v))*dx) - bc1 = DirichletBC(W.sub(0), Constant((0.0, 0.0)), 9) + bc1 = DirichletBC(W.sub(0), 0, 9) L0 = assemble((inner(sigma, tau) - inner(u, div(tau)) + inner(div(sigma), v))*dx, bcs=[bc1]) dV1 = V1.dof_count diff --git a/tests/submesh/test_submesh_interpolate.py b/tests/submesh/test_submesh_interpolate.py index f2d94459e3..27c7414d4a 100644 --- a/tests/submesh/test_submesh_interpolate.py +++ b/tests/submesh/test_submesh_interpolate.py @@ -21,9 +21,9 @@ def _get_expr(V): x, y, z = SpatialCoordinate(m) else: raise NotImplementedError("Not implemented") - if V.ufl_element().value_shape == (): + if V.value_shape == (): return cos(x) + x * exp(y) + sin(z) - elif V.ufl_element().value_shape == (2, ): + elif V.value_shape == (2, ): return as_vector([cos(x), sin(y)]) @@ -50,7 +50,7 @@ def _test_submesh_interpolate_cell_cell(mesh, subdomain_cond, fe_fesub): # if there is no ambiguity on the subdomain boundary. # For testing, the following suffices. g.interpolate(f) - temp = Constant(999.) if V.ufl_element().value_shape == () else as_vector([Constant(999.)] * np.prod(V.ufl_element().value_shape, dtype=int)) + temp = Constant(999.*np.ones(V.value_shape)) g.interpolate(temp, subset=mesh.topology.cell_subset(label_value)) # pollute the data g.interpolate(gsub, subset=mesh.topology.cell_subset(label_value)) assert assemble(inner(g - f, g - f) * dx(label_value)).real < 1e-14