Mixed topology

cpp/demo/hyperelasticity/main.cpp

@@ -222,7 +222,7 @@ int main(int argc, char* argv[])
     // Construct appropriate Basix element for stress
     constexpr auto family = basix::element::family::P;
     const auto cell_type
-        = mesh::cell_type_to_basix_type(mesh->topology().cell_type());
+        = mesh::cell_type_to_basix_type(mesh->topology().cell_types()[0]);
     constexpr int k = 0;
     constexpr bool discontinuous = true;
 

cpp/demo/interpolation_different_meshes/main.cpp

@@ -33,13 +33,13 @@ int main(int argc, char* argv[])
                          mesh::CellType::hexahedron));
 
     basix::FiniteElement element_tet = basix::element::create_lagrange(
-        mesh::cell_type_to_basix_type(mesh_tet->topology().cell_type()), 1,
+        mesh::cell_type_to_basix_type(mesh_tet->topology().cell_types()[0]), 1,
         basix::element::lagrange_variant::equispaced, false);
     auto V_tet = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(mesh_tet, element_tet, 3));
 
     basix::FiniteElement element_hex = basix::element::create_lagrange(
-        mesh::cell_type_to_basix_type(mesh_hex->topology().cell_type()), 2,
+        mesh::cell_type_to_basix_type(mesh_hex->topology().cell_types()[0]), 2,
         basix::element::lagrange_variant::equispaced, false);
     auto V_hex = std::make_shared<fem::FunctionSpace<double>>(
         fem::create_functionspace(mesh_hex, element_hex, 3));

cpp/demo/mixed_topology/CMakeLists.txt

@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 3.16)
+project(mixed_topology_testing)
+
+# Find DOLFINx config file
+find_package(DOLFINX REQUIRED)
+
+# Set C++20 standard
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+include(CheckSymbolExists)
+set(CMAKE_REQUIRED_INCLUDES ${PETSC_INCLUDE_DIRS})
+check_symbol_exists(PETSC_USE_COMPLEX petscsystypes.h PETSC_SCALAR_COMPLEX)
+
+# Add target to compile UFL files
+if(PETSC_SCALAR_COMPLEX EQUAL 1)
+  set(SCALAR_TYPE "--scalar_type=double _Complex")
+endif()
+
+set(CMAKE_INCLUDE_CURRENT_DIR ON)
+
+add_executable(${PROJECT_NAME} main.cpp)
+target_link_libraries(${PROJECT_NAME} dolfinx)

cpp/demo/mixed_topology/main.cpp

@@ -0,0 +1,149 @@
+#include <dolfinx/common/IndexMap.h>
+#include <dolfinx/common/MPI.h>
+#include <dolfinx/graph/AdjacencyList.h>
+#include <dolfinx/mesh/Geometry.h>
+#include <dolfinx/mesh/Mesh.h>
+#include <dolfinx/mesh/Topology.h>
+#include <dolfinx/mesh/cell_types.h>
+#include <iostream>
+#include <vector>
+
+// Note: this demo is not currently intended to provide a fully functional
+// example of using a mixed-topology mesh, but shows only the
+// basic constrution. Experimental.
+
+int main(int argc, char* argv[])
+{
+  dolfinx::init_logging(argc, argv);
+  MPI_Init(&argc, &argv);
+
+  const int nx_s = 2;
+  const int nx_t = 2;
+  const int ny = 4;
+
+  const int num_s = nx_s * ny;
+  const int num_t = 2 * nx_t * ny;
+
+  std::vector<double> x;
+  for (int i = 0; i < nx_s + nx_t + 1; ++i)
+  {
+    for (int j = 0; j < ny + 1; ++j)
+    {
+      x.push_back(i);
+      x.push_back(j);
+      x.push_back(0);
+    }
+  }
+
+  std::vector<std::int64_t> cells;
+  std::vector<std::int32_t> offsets{0};
+  for (int i = 0; i < nx_s + nx_t; ++i)
+  {
+    for (int j = 0; j < ny; ++j)
+    {
+      const int v_0 = j + i * (ny + 1);
+      const int v_1 = v_0 + 1;
+      const int v_2 = v_0 + ny + 1;
+      const int v_3 = v_0 + ny + 2;
+
+      if (i < nx_s)
+      {
+        cells.push_back(v_0);
+        cells.push_back(v_1);
+        cells.push_back(v_3);
+        cells.push_back(v_2);
+
+        offsets.push_back(cells.size());
+      }
+      else
+      {
+        cells.push_back(v_0);
+        cells.push_back(v_1);
+        cells.push_back(v_2);
+
+        offsets.push_back(cells.size());
+
+        cells.push_back(v_1);
+        cells.push_back(v_2);
+        cells.push_back(v_3);
+
+        offsets.push_back(cells.size());
+      }
+    }
+  }
+
+  dolfinx::graph::AdjacencyList<std::int64_t> cells_list(cells, offsets);
+
+  std::vector<std::int64_t> original_global_index(num_s + num_t);
+  std::iota(original_global_index.begin(), original_global_index.end(), 0);
+
+  std::vector<int> ghost_owners;
+
+  std::vector<std::int32_t> cell_group_offsets{0, num_s, num_s + num_t,
+                                               num_s + num_t, num_s + num_t};
+
+  std::vector<std::int64_t> boundary_vertices;
+  for (int j = 0; j < ny + 1; ++j)
+  {
+    boundary_vertices.push_back(j);
+    boundary_vertices.push_back(j + (nx_s + nx_t + 1) * ny);
+  }
+  for (int i = 0; i < nx_s + nx_t + 1; ++i)
+  {
+    boundary_vertices.push_back((ny + 1) * i);
+    boundary_vertices.push_back(ny + (ny + 1) * i);
+  }
+
+  std::sort(boundary_vertices.begin(), boundary_vertices.end());
+  boundary_vertices.erase(
+      std::unique(boundary_vertices.begin(), boundary_vertices.end()),
+      boundary_vertices.end());
+
+  std::vector<dolfinx::mesh::CellType> cell_types{
+      dolfinx::mesh::CellType::quadrilateral,
+      dolfinx::mesh::CellType::triangle};
+  std::vector<dolfinx::fem::CoordinateElement> elements;
+  for (auto ct : cell_types)
+    elements.push_back(dolfinx::fem::CoordinateElement(ct, 1));
+
+  {
+    auto topo = dolfinx::mesh::create_topology(
+        MPI_COMM_WORLD, cells_list, original_global_index, ghost_owners,
+        cell_types, cell_group_offsets, boundary_vertices);
+
+    auto topo_cells = topo.connectivity(2, 0);
+
+    for (int i = 0; i < topo_cells->num_nodes(); ++i)
+    {
+      std::cout << i << " [";
+      for (auto q : topo_cells->links(i))
+        std::cout << q << " ";
+      std::cout << "]\n";
+    }
+
+    topo.create_connectivity(1, 0);
+
+    auto topo_facets = topo.connectivity(1, 0);
+    for (int i = 0; i < topo_facets->num_nodes(); ++i)
+    {
+      std::cout << i << " [";
+      for (auto q : topo_facets->links(i))
+        std::cout << q << " ";
+      std::cout << "]\n";
+    }
+
+    auto geom = dolfinx::mesh::create_geometry(MPI_COMM_WORLD, topo, elements,
+                                               cells_list, x, 2);
+
+    dolfinx::mesh::Mesh mesh(MPI_COMM_WORLD, topo, geom);
+    std::cout << "num cells = " << mesh.topology().index_map(2)->size_local()
+              << "\n";
+    for (auto q : mesh.topology().entity_group_offsets(2))
+      std::cout << q << " ";
+    std::cout << "\n";
+  }
+
+  MPI_Finalize();
+
+  return 0;
+}

cpp/dolfinx/fem/DirichletBC.cpp

@@ -182,10 +182,15 @@ fem::locate_dofs_topological(mesh::Topology& topology, const DofMap& dofmap,
                              int dim, std::span<const std::int32_t> entities,
                              bool remote)
 {
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+  {
+    throw std::runtime_error("Multiple cell types in DirichletBC");
+  }
+
   // Prepare an element - local dof layout for dofs on entities of the
   // entity_dim
-  const int num_cell_entities
-      = mesh::cell_num_entities(topology.cell_type(), dim);
+  const int num_cell_entities = mesh::cell_num_entities(cell_types.back(), dim);
   std::vector<std::vector<int>> entity_dofs;
   for (int i = 0; i < num_cell_entities; ++i)
   {
@@ -298,9 +303,13 @@ std::array<std::vector<std::int32_t>, 2> fem::locate_dofs_topological(
   // Check that dof layouts are the same
   assert(dofmap0.element_dof_layout() == dofmap1.element_dof_layout());
 
+  auto cell_types = topology.cell_types();
+  if (cell_types.size() > 1)
+  {
+    throw std::runtime_error("Multiple cell types in DirichletBC");
+  }
   // Build vector of local dofs for each cell entity
-  const int num_cell_entities
-      = mesh::cell_num_entities(topology.cell_type(), dim);
+  const int num_cell_entities = mesh::cell_num_entities(cell_types.back(), dim);
   std::vector<std::vector<int>> entity_dofs;
   for (int i = 0; i < num_cell_entities; ++i)
   {

cpp/dolfinx/fem/DofMap.cpp

@@ -233,8 +233,8 @@ std::pair<DofMap, std::vector<std::int32_t>> DofMap::collapse(
       // Create new element dof layout and reset parent
       ElementDofLayout collapsed_dof_layout = layout.copy();
 
-      auto [_index_map, bs, dofmap]
-          = build_dofmap_data(comm, topology, collapsed_dof_layout, reorder_fn);
+      auto [_index_map, bs, dofmap] = build_dofmap_data(
+          comm, topology, {collapsed_dof_layout}, reorder_fn);
       auto index_map
           = std::make_shared<common::IndexMap>(std::move(_index_map));
       return DofMap(layout, index_map, bs, std::move(dofmap), bs);

cpp/dolfinx/fem/Expression.h

@@ -141,7 +141,12 @@ class Expression
     // Prepare cell geometry
     const graph::AdjacencyList<std::int32_t>& x_dofmap
         = _mesh->geometry().dofmap();
-    const std::size_t num_dofs_g = _mesh->geometry().cmap().dim();
+
+    // Get geometry data
+    auto cmaps = _mesh->geometry().cmaps();
+    assert(cmaps.size() == 1);
+
+    const std::size_t num_dofs_g = cmaps.back().dim();
     std::span<const U> x_g = _mesh->geometry().x();
 
     // Create data structures used in evaluation

cpp/dolfinx/fem/FiniteElement.h

@@ -168,7 +168,7 @@ class FiniteElement
   /// @return True if the map is the identity
   bool map_ident() const noexcept;
 
-  /// @brief Points on the reference cell at which an expression need to
+  /// @brief Points on the reference cell at which an expression needs to
   /// be evaluated in order to interpolate the expression in the finite
   /// element space.
   ///

cpp/dolfinx/fem/Function.h

@@ -393,15 +393,20 @@ class Function
     const std::size_t tdim = mesh->topology().dim();
     auto map = mesh->topology().index_map(tdim);
 
+    // Get coordinate map
+    if (mesh->geometry().cmaps().size() > 1)
+    {
+      throw std::runtime_error(
+          "Function with multiple geometry maps not implemented.");
+    }
+    const CoordinateElement& cmap = mesh->geometry().cmaps()[0];
+
     // Get geometry data
     const graph::AdjacencyList<std::int32_t>& x_dofmap
         = mesh->geometry().dofmap();
-    const std::size_t num_dofs_g = mesh->geometry().cmap().dim();
+    const std::size_t num_dofs_g = cmap.dim();
     std::span<const U> x_g = mesh->geometry().x();
 
-    // Get coordinate map
-    const CoordinateElement& cmap = mesh->geometry().cmap();
-
     // Get element
     std::shared_ptr<const FiniteElement> element = _function_space->element();
     assert(element);

cpp/dolfinx/fem/FunctionSpace.h

@@ -216,12 +216,16 @@ class FunctionSpace
     const auto [X, Xshape] = _element->interpolation_points();
 
     // Get coordinate map
-    const CoordinateElement& cmap = _mesh->geometry().cmap();
+    if (_mesh->geometry().cmaps().size() > 1)
+    {
+      throw std::runtime_error("Mixed topology not supported");
+    }
+    const CoordinateElement& cmap = _mesh->geometry().cmaps()[0];
 
     // Prepare cell geometry
     const graph::AdjacencyList<std::int32_t>& x_dofmap
         = _mesh->geometry().dofmap();
-    const std::size_t num_dofs_g = _mesh->geometry().cmap().dim();
+    const std::size_t num_dofs_g = cmap.dim();
     std::span<const T> x_g = _mesh->geometry().x();
 
     // Array to hold coordinates to return

cpp/dolfinx/fem/assemble_matrix_impl.h

@@ -48,7 +48,7 @@ void assemble_cells(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Iterate over active cells
@@ -144,7 +144,7 @@ void assemble_exterior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in assembly
@@ -240,7 +240,7 @@ void assemble_interior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Data structures used in assembly
@@ -445,7 +445,11 @@ void assemble_matrix(
     else
       get_perm = [](std::size_t) { return 0; };
 
-    int num_cell_facets = mesh::cell_num_entities(mesh->topology().cell_type(),
+    auto cell_types = mesh->topology().cell_types();
+    if (cell_types.size() > 1)
+      throw std::runtime_error("Multiple cell types in the assembler.");
+
+    int num_cell_facets = mesh::cell_num_entities(cell_types.back(),
                                                   mesh->topology().dim() - 1);
     const std::vector<int> c_offsets = a.coefficient_offsets();
     for (int i : a.integral_ids(IntegralType::interior_facet))

cpp/dolfinx/fem/assemble_scalar_impl.h

@@ -34,7 +34,7 @@ T assemble_cells(const mesh::Geometry<scalar_value_type_t<T>>& geometry,
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -74,7 +74,7 @@ T assemble_exterior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -117,7 +117,7 @@ T assemble_interior_facets(
 
   // Prepare cell geometry
   const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
-  const std::size_t num_dofs_g = geometry.cmap().dim();
+  const std::size_t num_dofs_g = geometry.cmaps()[0].dim();
   auto x = geometry.x();
 
   // Create data structures used in assembly
@@ -199,7 +199,12 @@ T assemble_scalar(
     const std::vector<std::uint8_t>& perms
         = mesh->topology().get_facet_permutations();
 
-    int num_cell_facets = mesh::cell_num_entities(mesh->topology().cell_type(),
+    auto cell_types = mesh->topology().cell_types();
+    if (cell_types.size() > 1)
+    {
+      throw std::runtime_error("MUltiple cell types in the assembler");
+    }
+    int num_cell_facets = mesh::cell_num_entities(cell_types.back(),
                                                   mesh->topology().dim() - 1);
     const std::vector<int> c_offsets = M.coefficient_offsets();
     for (int i : M.integral_ids(IntegralType::interior_facet))