Favor exterior_facet_indices over locate_entities_boundary for re…

…trieving complete boundary (FEniCS#3283)

* Introduce optionally empty marker for boundary entities

* Adapt to optional boundary marker

* Make use of simplifications in python code

* Make use of simplifications in cpp code

* Add xfail for now (FEniCS#3284)

* Data independent form compilation for Python interface (FEniCS#3263)

* Create mesh independent form compilator for Python interface

* Ruff formatting

* Mypy

* Parameterize test over all dtypes

* Ruff format

* Fix compile options handling

* Fix coeff and constant name mapping

* Ruff

* Remove extra definition of basix element and check that mesh c_el is consistent

* Add support for linear and bilinear forms

* Fix typo

* Apply suggestions from code review

Co-authored-by: Chris Richardson <chris@bpi.cam.ac.uk>

* Apply suggestions from code review

Co-authored-by: Chris Richardson <chris@bpi.cam.ac.uk>

* Apply suggestions from code review

* Apply suggestions from code review

Co-authored-by: Igor Baratta <igorbaratta@gmail.com>

* Apply suggestions from code review

* Add documentation and default scalar type

---------

Co-authored-by: Chris Richardson <chris@bpi.cam.ac.uk>
Co-authored-by: Igor Baratta <igorbaratta@gmail.com>

* Bump docker/build-push-action from 5 to 6 (FEniCS#3279)

Bumps [docker/build-push-action](https://github.com/docker/build-push-action) from 5 to 6.
- [Release notes](https://github.com/docker/build-push-action/releases)
- [Commits](docker/build-push-action@v5...v6)

---
updated-dependencies:
- dependency-name: docker/build-push-action
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>

* Improve check for facet->cell connectivity (FEniCS#3281)

* Improve check

* Update utils.cpp

* Revert optional and use

* Ruff

* Typo

* another

* Fix biharmonic

* Fix python biharmonic

* Fix python demo hdg

* last?

* fix test fem pipeline

* fem pipeline the last?

* Remove workarounds following upstream updates for CI (FEniCS#3289)

* Remove temporary workarounds

* Use dict comprehension

* Export missing compile time constants to python (FEniCS#3285)

Co-authored-by: papel <xDNdnQoyTcrrMdWCVZft4YLqkaRWnRb@protonmail.com>

* Update python/demo/demo_biharmonic.py

Co-authored-by: Jørgen Schartum Dokken <dokken92@gmail.com>

* Apply suggestions

---------

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: Chris Richardson <chris@bpi.cam.ac.uk>
Co-authored-by: Jørgen Schartum Dokken <dokken92@gmail.com>
Co-authored-by: Igor Baratta <igorbaratta@gmail.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
Co-authored-by: papel <xDNdnQoyTcrrMdWCVZft4YLqkaRWnRb@protonmail.com>

cpp/demo/biharmonic/main.cpp

@@ -207,23 +207,7 @@ int main(int argc, char* argv[])
     //  follows:
 
     //  Define boundary condition
-    auto facets = mesh::locate_entities_boundary(
-        *mesh, 1,
-        [](auto x)
-        {
-          constexpr U eps = 1.0e-6;
-          std::vector<std::int8_t> marker(x.extent(1), false);
-          for (std::size_t p = 0; p < x.extent(1); ++p)
-          {
-            U x0 = x(0, p);
-            U x1 = x(1, p);
-            if (std::abs(x0) < eps or std::abs(x0 - 1) < eps)
-              marker[p] = true;
-            if (std::abs(x1) < eps or std::abs(x1 - 1) < eps)
-              marker[p] = true;
-          }
-          return marker;
-        });
+    auto facets = mesh::exterior_facet_indices(*mesh->topology());
     const auto bdofs = fem::locate_dofs_topological(
         *V->mesh()->topology_mutable(), *V->dofmap(), 1, facets);
     auto bc = std::make_shared<const fem::DirichletBC<T>>(0.0, bdofs, V);

python/demo/demo_biharmonic.py

@@ -126,8 +126,6 @@
 from mpi4py import MPI
 
 # +
-import numpy as np
-
 import dolfinx
 import ufl
 from dolfinx import fem, io, mesh, plot
@@ -165,14 +163,9 @@
 # function that returns `True` for points `x` on the boundary and
 # `False` otherwise.
 
-facets = mesh.locate_entities_boundary(
-    msh,
-    dim=1,
-    marker=lambda x: np.isclose(x[0], 0.0)
-    | np.isclose(x[0], 1.0)
-    | np.isclose(x[1], 0.0)
-    | np.isclose(x[1], 1.0),
-)
+tdim = msh.topology.dim
+msh.topology.create_connectivity(tdim - 1, tdim)
+facets = mesh.exterior_facet_indices(msh.topology)
 
 # We now find the degrees-of-freedom that are associated with the
 # boundary facets using {py:func}`locate_dofs_topological

python/demo/demo_hdg.py

@@ -82,19 +82,6 @@ def u_e(x):
     return u_e
 
 
-def boundary(x):
-    """Boundary marker."""
-    lr = np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0)
-    tb = np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)
-    lrtb = lr | tb
-    if tdim == 2:
-        return lrtb
-    else:
-        assert tdim == 3
-        fb = np.isclose(x[2], 0.0) | np.isclose(x[2], 1.0)
-        return lrtb | fb
-
-
 comm = MPI.COMM_WORLD
 rank = comm.rank
 dtype = PETSc.ScalarType
@@ -187,7 +174,7 @@ def boundary(x):
 
 # Apply Dirichlet boundary conditions
 # We begin by locating the boundary facets of msh
-msh_boundary_facets = mesh.locate_entities_boundary(msh, fdim, boundary)
+msh_boundary_facets = mesh.exterior_facet_indices(msh.topology)
 # Since the boundary condition is enforced in the facet space, we must
 # use the mesh_to_facet_mesh map to get the corresponding facets in
 # facet_mesh

python/demo/demo_navier-stokes.py

@@ -261,15 +261,6 @@ def f_expr(x):
     return np.vstack((np.zeros_like(x[0]), np.zeros_like(x[0])))
 
 
-def boundary_marker(x):
-    return (
-        np.isclose(x[0], 0.0)
-        | np.isclose(x[0], 1.0)
-        | np.isclose(x[1], 0.0)
-        | np.isclose(x[1], 1.0)
-    )
-
-
 # -
 
 # We define some simulation parameters
@@ -344,7 +335,7 @@ def jump(phi, n):
 L = fem.form([L_0, L_1])
 
 # Boundary conditions
-boundary_facets = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, boundary_marker)
+boundary_facets = mesh.exterior_facet_indices(msh.topology)
 boundary_vel_dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, boundary_facets)
 bc_u = fem.dirichletbc(u_D, boundary_vel_dofs)
 bcs = [bc_u]

python/test/unit/fem/test_assemble_submesh.py

@@ -20,6 +20,7 @@
     create_submesh,
     create_unit_cube,
     create_unit_square,
+    exterior_facet_indices,
     locate_entities,
     locate_entities_boundary,
     meshtags,
@@ -301,14 +302,8 @@ def test_mixed_dom_codim_1(n, k):
 
     # Create a submesh of the boundary
     tdim = msh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        msh,
-        tdim - 1,
-        lambda x: np.isclose(x[0], 0.0)
-        | np.isclose(x[0], 1.0)
-        | np.isclose(x[1], 0.0)
-        | np.isclose(x[1], 1.0),
-    )
+    msh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(msh.topology)
 
     smsh, smsh_to_msh = create_submesh(msh, tdim - 1, boundary_facets)[:2]
 

python/test/unit/fem/test_bcs.py

@@ -27,7 +27,7 @@
     locate_dofs_topological,
     set_bc,
 )
-from dolfinx.mesh import CellType, create_unit_cube, create_unit_square, locate_entities_boundary
+from dolfinx.mesh import CellType, create_unit_cube, create_unit_square, exterior_facet_indices
 from ufl import dx, inner
 
 
@@ -119,9 +119,8 @@ def test_constant_bc_constructions():
     V2 = functionspace(msh, ("Lagrange", 1, (gdim, gdim)))
 
     tdim = msh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        msh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
-    )
+    msh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(msh.topology)
     boundary_dofs0 = locate_dofs_topological(V0, tdim - 1, boundary_facets)
     boundary_dofs1 = locate_dofs_topological(V1, tdim - 1, boundary_facets)
     boundary_dofs2 = locate_dofs_topological(V2, tdim - 1, boundary_facets)
@@ -166,10 +165,8 @@ def test_constant_bc(mesh_factory):
     V = functionspace(mesh, ("Lagrange", 1))
     c = default_scalar_type(2)
     tdim = mesh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        mesh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
-    )
-
+    mesh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(mesh.topology)
     boundary_dofs = locate_dofs_topological(V, tdim - 1, boundary_facets)
 
     u_bc = Function(V)
@@ -205,9 +202,8 @@ def test_vector_constant_bc(mesh_factory):
     V = functionspace(mesh, ("Lagrange", 1, (gdim,)))
     assert V.num_sub_spaces == gdim
     c = np.arange(1, mesh.geometry.dim + 1, dtype=default_scalar_type)
-    boundary_facets = locate_entities_boundary(
-        mesh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
-    )
+    mesh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(mesh.topology)
 
     # Set using sub-functions
     Vs = [V.sub(i).collapse()[0] for i in range(V.num_sub_spaces)]
@@ -252,9 +248,8 @@ def test_sub_constant_bc(mesh_factory):
     V = functionspace(mesh, ("Lagrange", 1, (gdim,)))
     c = Constant(mesh, default_scalar_type(3.14))
     tdim = mesh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        mesh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
-    )
+    mesh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(mesh.topology)
 
     for i in range(V.num_sub_spaces):
         Vi = V.sub(i).collapse()[0]
@@ -288,9 +283,8 @@ def test_mixed_constant_bc(mesh_factory):
     func, args = mesh_factory
     mesh = func(*args)
     tdim = mesh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        mesh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
-    )
+    mesh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(mesh.topology)
     TH = mixed_element(
         [
             element("Lagrange", mesh.basix_cell(), 1, dtype=default_real_type),
@@ -330,9 +324,8 @@ def test_mixed_blocked_constant():
     Dirichlet BC based on a vector valued Constant."""
     mesh = create_unit_square(MPI.COMM_WORLD, 4, 4)
     tdim = mesh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        mesh, tdim - 1, lambda x: np.ones(x.shape[1], dtype=bool)
-    )
+    mesh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(mesh.topology)
 
     TH = mixed_element(
         [

python/test/unit/fem/test_fem_pipeline.py

@@ -35,7 +35,6 @@
     create_unit_cube,
     create_unit_square,
     exterior_facet_indices,
-    locate_entities_boundary,
 )
 from ufl import (
     CellDiameter,
@@ -250,9 +249,9 @@ def test_petsc_curl_curl_eigenvalue(family, order):
     a = inner(ufl.curl(u), ufl.curl(v)) * dx
     b = inner(u, v) * dx
 
-    boundary_facets = locate_entities_boundary(
-        mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True, dtype=bool)
-    )
+    tdim = mesh.topology.dim
+    mesh.topology.create_connectivity(tdim - 1, tdim)
+    boundary_facets = exterior_facet_indices(mesh.topology)
     boundary_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets)
 
     zero_u = Function(V)
@@ -377,9 +376,7 @@ def b(tau_S, v):
 
     # Strong (Dirichlet) boundary condition
     tdim = mesh.topology.dim
-    boundary_facets = locate_entities_boundary(
-        mesh, tdim - 1, lambda x: np.full(x.shape[1], True, dtype=bool)
-    )
+    boundary_facets = exterior_facet_indices(mesh.topology)
     boundary_dofs = locate_dofs_topological((V.sub(1), V_1), tdim - 1, boundary_facets)
 
     bcs = [dirichletbc(zero_u, boundary_dofs, V.sub(1))]