Favor exterior_facet_indices over locate_entities_boundary for re…

…trieving complete boundary (#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 (#3284) * Data independent form compilation for Python interface (#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 <[email protected]> * Apply suggestions from code review Co-authored-by: Chris Richardson <[email protected]> * Apply suggestions from code review * Apply suggestions from code review Co-authored-by: Igor Baratta <[email protected]> * Apply suggestions from code review * Add documentation and default scalar type --------- Co-authored-by: Chris Richardson <[email protected]> Co-authored-by: Igor Baratta <[email protected]> * Bump docker/build-push-action from 5 to 6 (#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] <[email protected]> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> * Improve check for facet->cell connectivity (#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 (#3289) * Remove temporary workarounds * Use dict comprehension * Export missing compile time constants to python (#3285) Co-authored-by: papel <[email protected]> * Update python/demo/demo_biharmonic.py Co-authored-by: Jørgen Schartum Dokken <[email protected]> * Apply suggestions --------- Signed-off-by: dependabot[bot] <[email protected]> Co-authored-by: Chris Richardson <[email protected]> Co-authored-by: Jørgen Schartum Dokken <[email protected]> Co-authored-by: Igor Baratta <[email protected]> Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> Co-authored-by: papel <[email protected]>
FEniCS · Jul 2, 2024 · 02e9a4d · 02e9a4d
1 parent b6864c0
commit 02e9a4d
Show file tree

Hide file tree

Showing 7 changed files with 26 additions and 86 deletions.
diff --git a/cpp/demo/biharmonic/main.cpp b/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);

diff --git a/python/demo/demo_biharmonic.py b/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

diff --git a/python/demo/demo_hdg.py b/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

diff --git a/python/demo/demo_navier-stokes.py b/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]

diff --git a/python/test/unit/fem/test_assemble_submesh.py b/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]
 

diff --git a/python/test/unit/fem/test_bcs.py b/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(
         [

diff --git a/python/test/unit/fem/test_fem_pipeline.py b/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))]