Expose function MeshTags.find to Python

cpp/dolfinx/mesh/MeshTags.h

@@ -79,7 +79,7 @@ class MeshTags
 
   /// @brief Find all entities with a given tag value
   /// @param[in] value The value
-  /// @return Indices of tagged entities
+  /// @return Indices of tagged entities. The indices are sorted.
   std::vector<std::int32_t> find(const T value) const
   {
     int n = std::count(_values.begin(), _values.end(), value);

python/demo/demo_pyvista.py

@@ -94,7 +94,7 @@ def plot_scalar():
 
 def plot_meshtags():
 
-    msh = create_unit_square(MPI.COMM_WORLD, 12, 12, cell_type=CellType.quadrilateral)
+    msh = create_unit_square(MPI.COMM_WORLD, 25, 25, cell_type=CellType.quadrilateral)
 
     # We continue using the mesh from the previous section, and find all
     # cells satisfying the condition below
@@ -128,8 +128,7 @@ def in_circle(x):
     # We can also visualize subsets of data, by creating a smaller topology,
     # only consisting of those entities that has value one in the
     # mesh tags
-    cells, types, x = plot.create_vtk_mesh(
-        msh, entities=cell_tags.indices[cell_tags.values == 1])
+    cells, types, x = plot.create_vtk_mesh(msh, entities=cell_tags.find(1))
 
     # We add this grid to the second plotter
     sub_grid = pyvista.UnstructuredGrid(cells, types, x)
@@ -169,12 +168,11 @@ def in_circle(x):
     # We start by interpolating a discontinuous function into a second order
     # discontinuous Lagrange space. Note that we use the `cell_tags` from
     # the previous section to get the cells for each of the regions
-    cells0 = cell_tags.indices[cell_tags.values == 0]
-    cells1 = cell_tags.indices[cell_tags.values == 1]
     V = FunctionSpace(msh, ("Discontinuous Lagrange", 2))
     u = Function(V, dtype=np.float64)
-    u.interpolate(lambda x: x[0], cells0)
-    u.interpolate(lambda x: x[1] + 1, cells1)
+    u.interpolate(lambda x: x[0], cell_tags.find(0))
+    u.interpolate(lambda x: x[1] + 1, cell_tags.find(1))
+    u.x.scatter_forward()
 
     # To get a topology that has a 1-1 correspondence with the
     # degrees-of-freedom in the function space, we call

python/dolfinx/wrappers/mesh.cpp

@@ -123,7 +123,9 @@ void declare_meshtags(py::module& m, std::string type)
                                return py::array_t<std::int32_t>(
                                    self.indices().size(), self.indices().data(),
                                    py::cast(self));
-                             });
+                             })
+      .def("find", [](dolfinx::mesh::MeshTags<T>& self, T value)
+           { return as_pyarray(self.find(value)); });
 
   m.def("create_meshtags",
         [](const std::shared_ptr<const dolfinx::mesh::Mesh>& mesh, int dim,