Add mixed-topology support to `locate_dofs_topological` by jpdean · Pull Request #4167 · FEniCS/dolfinx

jpdean · 2026-04-21T11:13:33Z

This PR adds mixed-topology support tolocate_dofs_topological.

jorgensd · 2026-04-21T17:10:21Z

+        throw std::runtime_error("Cell-to-entity connectivity of dimension "
+                                 + std::to_string(tdim) + "->"
+                                 + std::to_string(dim) + " doesn't exist");


Suggested change

throw std::runtime_error("Cell-to-entity connectivity of dimension "

+ std::to_string(tdim) + "->"

+ std::to_string(dim) + " doesn't exist");

throw std::runtime_error(std::format("Cell-to-entity connectivity of dimension {}->{} doesn't exist", tdim, dim));

jorgensd · 2026-04-21T17:12:36Z

+/// `topology.entity_types(dim)`
+/// @returns A list of `(cell_index, entity_index, cell_type_index)` tuples
+/// for each input entity. `cell_type_index` is the index of the cell type in
+/// `topology.cell_types()`.


Heads up, @chrisrichardson is rewriting all of these to have i inputs, as in: #4169

but not in Topology just yet, so get this merged first!

jorgensd · 2026-04-21T17:12:57Z

  const int tdim = topology.dim();
-  auto e_to_c = topology.connectivity(dim, tdim);
-  if (!e_to_c)
+  std::vector<mesh::CellType> cell_types = topology.cell_types();


Suggested change

std::vector<mesh::CellType> cell_types = topology.cell_types();

const std::vector<mesh::CellType>& cell_types = topology.cell_types();

jorgensd · 2026-04-21T17:14:43Z

+      c_to_e_connectivities[i]
+          = topology.connectivity({tdim, i}, {dim, entity_type_index});
+
+      // NOTE: If e_to_c exists, c_to_e should also exist.


Is this assumption always true? I didn't think that was the case, i.e. (tdim, 0) exists at construction of a topology, but (0, tdim) does not.

true - maybe just reword this comment to state that both are required.

I think, with the current connectivity implementation in FEniCSx, this statement is always true because of the direction I stated it in. Your example, Jørgen, is assuming that the reverse statement also holds, which is not the case. However, I am happy to clarify the comment! The code doesn't actually rely on it anyway, I've included a check just in case!

jorgensd · 2026-04-21T17:15:56Z

+  // Get connectivities for each cell type.
+  std::int32_t num_entities_local = 0;
+  bool any_connectivity = false;
+  std::vector<std::shared_ptr<const dolfinx::graph::AdjacencyList<int32_t>>>


Should we rather use a reserve and pushback?

jorgensd · 2026-04-21T17:20:00Z

+    {
+      auto e_to_c = e_to_c_connectivities[i];
+
+      if (e_to_c and e_to_c->num_links(e) > 0)


Can it happen that an entity is not connected to a cell?, i.e. e_to_c->num_links(e) == 0?

jorgensd · 2026-04-21T17:22:41Z

+    {
+      throw std::runtime_error(
+          "Entity " + std::to_string(e)
+          + " is not connected to any cells or is out of bounds.");


Out of bounds is covered in Line 93 to 99?

jorgensd · 2026-04-21T17:23:36Z

+      throw std::runtime_error(
+          "Entity " + std::to_string(e)
+          + " is not connected to any cells or is out of bounds.");


Suggested change

throw std::runtime_error(

"Entity " + std::to_string(e)

+ " is not connected to any cells or is out of bounds.");

throw std::runtime_error(std::format("Entity {} is not connected to any cells or is out of bounds.", e));

jorgensd · 2026-04-21T17:24:47Z

+  const std::size_t num_cell_types = topology.cell_types().size();
+  if (num_cell_types != 1)


Suggested change

const std::size_t num_cell_types = topology.cell_types().size();

if (num_cell_types != 1)

if (topology.cell_types().size() != 1)

jorgensd · 2026-04-21T17:26:34Z

-    throw std::runtime_error(
-        "Cell-to-entity connectivity has not been computed. Missing dims "
-        + std::to_string(tdim) + "->" + std::to_string(dim));
+    throw std::runtime_error("No connectivity from entities of dimension "


Use std::format here as well.

jorgensd · 2026-04-21T17:27:59Z

  // the index map into account as they can differ
-  const int bs = dofmap.bs();
-  const int element_bs = dofmap.element_dof_layout().block_size();
+  const int bs = dofmaps.front().get().bs();


We should probably add a comment here that the block size of either dofmap should be the same, as pointed out by @chrisrichardson in #4165

jorgensd · 2026-04-21T17:28:26Z

+    std::span<const std::int32_t> entities, int entity_type_index, bool remote)
+{
+  std::vector<mesh::CellType> cell_types = topology.cell_types();
+  const int num_cell_types = cell_types.size();


See warning from sonarcloud.

jorgensd · 2026-04-21T17:36:51Z

  {
    // Work with blocks
-    for (auto [cell, entity_local_index] : entity_indices)
+    for (auto [cell, entity_local_index, cell_type_idx] : entity_indices)


Should we use c rather than cell to avoid the shadowing warning?

jorgensd · 2026-04-21T17:38:02Z

      // dofs
-      std::span<const std::int32_t> cell_dofs = dofmap.cell_dofs(cell);
-      for (int index : entity_dofs[entity_local_index])
+      std::span<const std::int32_t> cell_dofs


Seems like we should call dofmaps[cell_type+idx].get() once rather than twice.

jorgensd · 2026-04-21T17:39:11Z

    // found by other processes, e.g. a vertex dof on this process that
    // has no connected facets on the boundary.
-    auto map = dofmap.index_map;
+    auto map = dofmaps.front().get().index_map;


Same here and below, call dofmaps.front().get() once as a separate variable?

jorgensd

Some minor things. In general it looks good!:)

jpdean added 25 commits July 10, 2025 14:39

Get something working

af3d836

Pass entity type idx

226ed8e

Expose ent type idx

5b15fc5

Hacks to check BC is being set

4d70b07

Fix connectivity bug

5459f6a

Merge branch 'main' into jpdean/mixed_locate_dofs_topo

39696b3

Improve logic and fix bug

62576ab

Improve locate_dofs_topological

975aada

Get non-blocked case working

fb42204

Change name back

334cdd9

Put name back

6223c72

Tidy

96810c1

Uncomment code

60595f7

Change arg order

d26925a

Use reference wrapper

e631dfa

Use ref wrapper and add back old interface

555c1e4

Add FIXME

90143f1

Merge branch 'main' into jpdean/mixed_locate_dofs_topo

ec8cc1f

Ruff

1baa389

Fix for test

1bda932

Fix FIXME

24cff86

Merge branch 'main' into jpdean/mixed_locate_dofs_topo

dc82029

Docs

80ee1c1

Use front()

6b2b268

Fix when errors are thrown

4baf424

jpdean requested review from chrisrichardson and jorgensd April 21, 2026 13:53

Index type

df9bf8e

jorgensd reviewed Apr 21, 2026

View reviewed changes

jorgensd requested changes Apr 21, 2026

View reviewed changes

	std::vector<mesh::CellType> cell_types = topology.cell_types();
	const std::vector<mesh::CellType>& cell_types = topology.cell_types();

		const std::size_t num_cell_types = topology.cell_types().size();
		if (num_cell_types != 1)

	const std::size_t num_cell_types = topology.cell_types().size();
	if (num_cell_types != 1)
	if (topology.cell_types().size() != 1)

Uh oh!

Conversation

jpdean commented Apr 21, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jorgensd left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants