Further fixes to FiniteElement Python wrapper (#3554)

* Fix finite element signature

* Same issue with interpolation points

* Update API in tests

* Update demos

Author: jorgensd

python/demo/demo_axis.py

@@ -786,7 +786,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
     v_dg_el = element("DG", msh.basix_cell(), degree, shape=(3,), dtype=real_type)
     W = fem.functionspace(msh, v_dg_el)
     Es_dg = fem.Function(W)
-    Es_expr = fem.Expression(Esh, W.element.interpolation_points())
+    Es_expr = fem.Expression(Esh, W.element.interpolation_points)
     Es_dg.interpolate(Es_expr)
     with VTXWriter(msh.comm, "sols/Es.bp", Es_dg) as f:
         f.write(0.0)

python/demo/demo_elasticity.py

@@ -250,7 +250,7 @@ def σ(v):
 
 # +
 W = functionspace(msh, ("Discontinuous Lagrange", 0))
-sigma_vm_expr = Expression(sigma_vm, W.element.interpolation_points())
+sigma_vm_expr = Expression(sigma_vm, W.element.interpolation_points)
 sigma_vm_h = Function(W)
 sigma_vm_h.interpolate(sigma_vm_expr)
 # -

python/dolfinx/fem/element.py

@@ -225,7 +225,7 @@ def interpolation_points(self) -> npt.NDArray[np.floating]:
             the points will typically be the quadrature points used to evaluate moment degrees of
             freedom.
         """
-        return self._cpp_object.interpolation_points
+        return self._cpp_object.interpolation_points()
 
     @property
     def interpolation_ident(self) -> bool:

python/dolfinx/wrappers/fem.cpp

@@ -100,8 +100,7 @@ void declare_function_space(nb::module_& m, std::string type)
             [](dolfinx::fem::FiniteElement<T>* self,
                basix::FiniteElement<T>& element,
                std::optional<std::vector<std::size_t>> block_shape,
-               bool symmetric)
-            {
+               bool symmetric) {
               new (self) dolfinx::fem::FiniteElement<T>(element, block_shape,
                                                         symmetric);
             },
@@ -284,7 +283,7 @@ void declare_function_space(nb::module_& m, std::string type)
             nb::arg("x"), nb::arg("cell_permutations"), nb::arg("dim"))
         .def_prop_ro("needs_dof_transformations",
                      &dolfinx::fem::FiniteElement<T>::needs_dof_transformations)
-        .def("signature", &dolfinx::fem::FiniteElement<T>::signature);
+        .def_prop_ro("signature", &dolfinx::fem::FiniteElement<T>::signature);
   }
 }
 

python/test/unit/fem/test_discrete_operators.py

@@ -169,7 +169,7 @@ def test_gradient_interpolation(cell_type, p, q):
     u = Function(V, uvec, dtype=dtype)
     u.interpolate(lambda x: 2 * x[0] ** p + 3 * x[1] ** p)
 
-    grad_u = Expression(ufl.grad(u), W.element.interpolation_points(), dtype=dtype)
+    grad_u = Expression(ufl.grad(u), W.element.interpolation_points, dtype=dtype)
     w_expr = Function(W, dtype=dtype)
     w_expr.interpolate(grad_u)
 

python/test/unit/fem/test_dof_permuting.py

@@ -166,7 +166,7 @@ def test_dof_positions(cell_type, space_type):
     entities = {i: {} for i in range(1, tdim)}
     for cell in range(coord_dofs.shape[0]):
         # Push coordinates forward
-        X = V.element.interpolation_points()
+        X = V.element.interpolation_points
         xg = x_g[coord_dofs[cell], :tdim]
         x = cmap.push_forward(X, xg)
 

python/test/unit/fem/test_dofmap.py

@@ -327,7 +327,7 @@ def test_higher_order_coordinate_map(points, celltype, order):
     mesh = create_mesh(MPI.COMM_WORLD, cells, points, domain)
 
     V = functionspace(mesh, ("Lagrange", 2))
-    X = V.element.interpolation_points()
+    X = V.element.interpolation_points
     coord_dofs = mesh.geometry.dofmap
     x_g = mesh.geometry.x
     cmap = mesh.geometry.cmap
@@ -402,7 +402,7 @@ def test_higher_order_tetra_coordinate_map(order):
     )
     mesh = create_mesh(MPI.COMM_WORLD, cells, points, domain)
     V = functionspace(mesh, ("Lagrange", order))
-    X = V.element.interpolation_points()
+    X = V.element.interpolation_points
     x_dofs = mesh.geometry.dofmap
     x_g = mesh.geometry.x
 

python/test/unit/fem/test_expression.py

@@ -48,7 +48,7 @@ def test_rank0(dtype):
     f.interpolate(lambda x: x[0] ** 2 + 2.0 * x[1] ** 2)
 
     ufl_expr = ufl.grad(f)
-    points = vdP1.element.interpolation_points()
+    points = vdP1.element.interpolation_points
 
     compiled_expr = Expression(ufl_expr, points, dtype=dtype)
     num_cells = mesh.topology.index_map(2).size_local
@@ -98,7 +98,7 @@ def test_rank1_hdiv(dtype):
     RT1 = functionspace(mesh, ("RT", 2))
     f = ufl.TrialFunction(RT1)
 
-    points = vdP1.element.interpolation_points()
+    points = vdP1.element.interpolation_points
     compiled_expr = Expression(f, points, dtype=dtype)
     num_cells = mesh.topology.index_map(2).size_local
     array_evaluated = compiled_expr.eval(mesh, np.arange(num_cells, dtype=np.int32))
@@ -355,7 +355,7 @@ def test_expression_eval_cells_subset(dtype):
     u = dolfinx.fem.Function(V, dtype=dtype)
     u.x.array[:] = dofs_to_cells
     u.x.scatter_forward()
-    e = dolfinx.fem.Expression(u, V.element.interpolation_points())
+    e = dolfinx.fem.Expression(u, V.element.interpolation_points)
 
     # Test eval on single cell
     for c in range(cells_imap.size_local):
@@ -387,8 +387,8 @@ def test_expression_comm(dtype):
     mesh = create_unit_square(MPI.COMM_WORLD, 4, 4, dtype=xtype)
     v = Constant(mesh, dtype(1))
     u = Function(functionspace(mesh, ("Lagrange", 1)), dtype=dtype)
-    Expression(v, u.function_space.element.interpolation_points(), comm=MPI.COMM_WORLD)
-    Expression(v, u.function_space.element.interpolation_points(), comm=MPI.COMM_SELF)
+    Expression(v, u.function_space.element.interpolation_points, comm=MPI.COMM_WORLD)
+    Expression(v, u.function_space.element.interpolation_points, comm=MPI.COMM_SELF)
 
 
 def compute_exterior_facet_entities(mesh, facets):

python/test/unit/fem/test_function_space.py

@@ -123,7 +123,7 @@ def test_sub(Q, W):
     assert W.element.num_sub_elements == X.element.num_sub_elements
     assert W.element.space_dimension == X.element.space_dimension
     assert W.value_shape == X.value_shape
-    assert W.element.interpolation_points().shape == X.element.interpolation_points().shape
+    assert W.element.interpolation_points.shape == X.element.interpolation_points.shape
     assert W.element == X.element
 
 

python/test/unit/fem/test_interpolation.py

@@ -618,7 +618,7 @@ def test_nedelec_spatial(order, dim):
     # The expression (x,y,z) is contained in the N1curl function space
     # order>1
     f_ex = x
-    f = Expression(f_ex, V.element.interpolation_points())
+    f = Expression(f_ex, V.element.interpolation_points)
     u.interpolate(f)
     assert np.abs(assemble_scalar(form(ufl.inner(u - f_ex, u - f_ex) * ufl.dx))) == pytest.approx(
         0, abs=1e-10
@@ -628,7 +628,7 @@ def test_nedelec_spatial(order, dim):
     # order
     V2 = functionspace(mesh, ("N2curl", 1))
     w = Function(V2)
-    f2 = Expression(f_ex, V2.element.interpolation_points())
+    f2 = Expression(f_ex, V2.element.interpolation_points)
     w.interpolate(f2)
     assert np.abs(assemble_scalar(form(ufl.inner(w - f_ex, w - f_ex) * ufl.dx))) == pytest.approx(0)
 
@@ -650,7 +650,7 @@ def test_vector_interpolation_spatial(order, dim, affine):
 
     # The expression (x,y,z)^n is contained in space
     f = ufl.as_vector([x[i] ** order for i in range(dim)])
-    u.interpolate(Expression(f, V.element.interpolation_points()))
+    u.interpolate(Expression(f, V.element.interpolation_points))
     assert np.abs(assemble_scalar(form(ufl.inner(u - f, u - f) * ufl.dx))) == pytest.approx(0)
 
 
@@ -663,7 +663,7 @@ def test_2D_lagrange_to_curl(order):
     u1 = Function(W)
     u1.interpolate(lambda x: x[0])
     f = ufl.as_vector((u0, u1))
-    f_expr = Expression(f, V.element.interpolation_points())
+    f_expr = Expression(f, V.element.interpolation_points)
     u.interpolate(f_expr)
     x = ufl.SpatialCoordinate(mesh)
     f_ex = ufl.as_vector((-x[1], x[0]))
@@ -679,7 +679,7 @@ def test_de_rahm_2D(order):
     g = ufl.grad(w)
     Q = functionspace(mesh, ("N2curl", order - 1))
     q = Function(Q)
-    q.interpolate(Expression(g, Q.element.interpolation_points()))
+    q.interpolate(Expression(g, Q.element.interpolation_points))
     x = ufl.SpatialCoordinate(mesh)
     g_ex = ufl.as_vector((1 + x[1], 4 * x[1] + x[0]))
     assert np.abs(assemble_scalar(form(ufl.inner(q - g_ex, q - g_ex) * ufl.dx))) == pytest.approx(
@@ -692,7 +692,7 @@ def test_de_rahm_2D(order):
     def curl2D(u):
         return ufl.as_vector((ufl.Dx(u[1], 0), -ufl.Dx(u[0], 1)))
 
-    v.interpolate(Expression(curl2D(ufl.grad(w)), V.element.interpolation_points()))
+    v.interpolate(Expression(curl2D(ufl.grad(w)), V.element.interpolation_points))
     h_ex = ufl.as_vector((1, -1))
     assert np.abs(assemble_scalar(form(ufl.inner(v - h_ex, v - h_ex) * ufl.dx))) == pytest.approx(
         0, abs=np.sqrt(np.finfo(mesh.geometry.x.dtype).eps)
@@ -720,7 +720,7 @@ def test_interpolate_subset(order, dim, affine, callable_):
     x = ufl.SpatialCoordinate(mesh)
     f = x[1] ** order
     if not callable_:
-        expr = Expression(f, V.element.interpolation_points())
+        expr = Expression(f, V.element.interpolation_points)
         u.interpolate(expr, cells_local)
     else:
         u.interpolate(lambda x: x[1] ** order, cells_local)
@@ -760,7 +760,7 @@ def test_interpolate_callable_subset(bound):
     u0, u1 = Function(V), Function(V)
     x = ufl.SpatialCoordinate(mesh)
     f = x[0]
-    expr = Expression(f, V.element.interpolation_points())
+    expr = Expression(f, V.element.interpolation_points)
     u0.interpolate(lambda x: x[0], cells_local)
     u1.interpolate(expr, cells_local)
     assert np.allclose(u0.x.array, u1.x.array, rtol=1.0e-6, atol=1.0e-6)
@@ -1114,7 +1114,7 @@ def modified_grad(x):
     V_sub = functionspace(submesh, ("N2curl", 1))
     u_sub = Function(V_sub)
 
-    parent_expr = Expression(ufl.grad(u), V_sub.element.interpolation_points())
+    parent_expr = Expression(ufl.grad(u), V_sub.element.interpolation_points)
 
     # Map from parent to sub mesh
 
@@ -1136,7 +1136,7 @@ def modified_grad(x):
 
     # Map exact solution (based on quadrature points) back to parent mesh
     sub_vec = ufl.as_vector((-0.2 * u_sub_exact[1], 0.1 * u_sub_exact[0]))
-    sub_expr = Expression(sub_vec, W.element.interpolation_points())
+    sub_expr = Expression(sub_vec, W.element.interpolation_points)
 
     # Mapping back needs to be restricted to the subset of cells in the submesh
     w.interpolate(sub_expr, cells=sub_to_parent, expr_mesh=submesh, cell_map=parent_to_sub)

python/test/unit/fem/test_petsc_discrete_operators.py

@@ -103,7 +103,7 @@ def test_gradient_interpolation_petsc(self, cell_type, p, q):
         u = Function(V)
         u.interpolate(lambda x: 2 * x[0] ** p + 3 * x[1] ** p)
 
-        grad_u = Expression(ufl.grad(u), W.element.interpolation_points())
+        grad_u = Expression(ufl.grad(u), W.element.interpolation_points)
         w_expr = Function(W)
         w_expr.interpolate(grad_u)