Remove gdim from element input and update FunctionSpace C++ API (#105)

* Fixes related to FEniCS/basix#772 #702

* More gdims to remove

Author: jorgensd

cpp/PeriodicConstraint.h

@@ -64,8 +64,7 @@ dolfinx_mpc::mpc_data<T> _create_periodic_condition(
       parents_glob[i] = parents_glob[i] * bs + parent_rems[i];
     return parents_glob;
   };
-  if (const std::size_t value_size
-      = V.element()->value_size() / V.element()->block_size();
+  if (const std::size_t value_size = V.value_size() / V.element()->block_size();
       value_size > 1)
     throw std::runtime_error(
         "Periodic conditions for vector valued spaces are not "

cpp/mpc_helpers.h

@@ -233,7 +233,10 @@ create_extended_functionspace(const dolfinx::fem::FunctionSpace<U>& V,
       old_dofmap.element_dof_layout(), new_index_map, old_dofmap.bs(),
       std::move(flattened_dofmap), old_dofmap.bs());
 
-  return dolfinx::fem::FunctionSpace(V.mesh(), element, new_dofmap);
+  return dolfinx::fem::FunctionSpace(
+      V.mesh(), element, new_dofmap,
+      dolfinx::fem::compute_value_shape(element, V.mesh()->topology()->dim(),
+                                        V.mesh()->geometry().dim()));
 }
 
 } // namespace dolfinx_mpc

cpp/utils.h

@@ -1018,7 +1018,7 @@ evaluate_basis_functions(const dolfinx::fem::FunctionSpace<U>& V,
 
   assert(basis_shape[2]
          == std::size_t(element->space_dimension() / bs_element));
-  assert(basis_shape[3] == std::size_t(element->value_size() / bs_element));
+  assert(basis_shape[3] == std::size_t(V.value_size() / bs_element));
   std::array<std::size_t, 3> reference_shape
       = {basis_shape[1], basis_shape[2], basis_shape[3]};
   std::vector<U> output_basis(std::reduce(

python/benchmarks/bench_contact_3D.py

@@ -92,7 +92,7 @@ def over_plane(p0, p1, p2):
         # Concatenate points and cells
         points = np.vstack([mesh0.geometry.x, mesh1.geometry.x])
         cells = np.vstack([cells0, cells1])
-        domain = Mesh(element("Lagrange", ct.name, 1, shape=(points.shape[1],), gdim=points.shape[1]))
+        domain = Mesh(element("Lagrange", ct.name, 1, shape=(points.shape[1],)))
         # Rotate mesh
         points = np.dot(r_matrix, points.T).T.astype(default_real_type)
 

python/demos/create_and_export_mesh.py

@@ -394,7 +394,7 @@ def over_line(p0, p1):
         cells1 += mesh0.geometry.x.shape[0]
 
         cells = np.vstack([cells0, cells1])
-        domain = ufl.Mesh(element("Lagrange", ct.name, 1, shape=(points.shape[1],), gdim=points.shape[1]))
+        domain = ufl.Mesh(element("Lagrange", ct.name, 1, shape=(points.shape[1],)))
         mesh = _mesh.create_mesh(MPI.COMM_SELF, cells, points, domain)
         tdim = mesh.topology.dim
         fdim = tdim - 1
@@ -514,7 +514,7 @@ def over_plane(p0, p1, p2):
         cells1 += mesh0.geometry.x.shape[0]
 
         cells = np.vstack([cells0, cells1])
-        domain = ufl.Mesh(element("Lagrange", ct.name, 1, shape=(points.shape[1],), gdim=points.shape[1]))
+        domain = ufl.Mesh(element("Lagrange", ct.name, 1, shape=(points.shape[1],)))
         mesh = _mesh.create_mesh(MPI.COMM_SELF, cells, points, domain)
         tdim = mesh.topology.dim
         fdim = tdim - 1

python/demos/demo_elasticity_disconnect.py

@@ -201,8 +201,7 @@ def sigma(v):
 # Write solution to file
 V_out = functionspace(mesh, basix.ufl.element("Lagrange", mesh.topology.cell_name(),
                       mesh.geometry.cmap.degree,
-                      lagrange_variant=basix.LagrangeVariant(mesh.geometry.cmap.variant),
-                      gdim=mesh.geometry.dim, shape=(V.dofmap.bs,)))
+                      lagrange_variant=basix.LagrangeVariant(mesh.geometry.cmap.variant), shape=(V.dofmap.bs,)))
 u_out = Function(V_out)
 u_out.interpolate(u_h)
 u_out.name = "uh"

python/demos/demo_stokes.py

@@ -123,8 +123,8 @@ def create_mesh_gmsh(L: int = 2, H: int = 1, res: float = 0.1, theta: float = np
 # The next step is the create the function spaces for the fluid velocit and pressure.
 # We will use a mixed-formulation, and we use `basix.ufl` to create the Taylor-Hood finite element pair
 
-P2 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 2, gdim=mesh.geometry.dim, shape=(mesh.geometry.dim, ))
-P1 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1, gdim=mesh.geometry.dim)
+P2 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim, ))
+P1 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1)
 
 TH = basix.ufl.mixed_element([P2, P1])
 W = fem.functionspace(mesh, TH)

python/dolfinx_mpc/multipointconstraint.py

@@ -249,7 +249,7 @@ def create_slip_constraint(self, space: _fem.FunctionSpace, facet_marker: Tuple[
                 cellname = mesh.ufl_cell().cellname()
                 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
                 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
-                me = basix.ufl.mixed_element([Ve, Qe], gdim=mesh.geometry.dim)
+                me = basix.ufl.mixed_element([Ve, Qe])
                 W = dolfinx.fem.functionspace(mesh, me)
                 mpc = MultiPointConstraint(W)
                 n_space, _ = W.sub(0).collapse()
@@ -265,7 +265,7 @@ def create_slip_constraint(self, space: _fem.FunctionSpace, facet_marker: Tuple[
                 cellname = mesh.ufl_cell().cellname()
                 Ve = basix.ufl.element(basix.ElementFamily.P, cellname , 2, shape=(mesh.geometry.dim,))
                 Qe = basix.ufl.element(basix.ElementFamily.P, cellname , 1)
-                me = basix.ufl.mixed_element([Ve, Qe], gdim=mesh.geometry.dim)
+                me = basix.ufl.mixed_element([Ve, Qe])
                 W = dolfinx.fem.functionspace(mesh, me)
                 mpc = MultiPointConstraint(W)
                 n_space, _ = W.sub(0).collapse()