Merge branch 'main' into checkpointing

cpp/demo/biharmonic/CMakeLists.txt

@@ -45,9 +45,7 @@ add_custom_command(
 
 set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
-add_executable(
-  ${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/biharmonic.c
-)
+add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/biharmonic.c)
 target_link_libraries(${PROJECT_NAME} dolfinx)
 
 # Do not throw error for 'multi-line comments' (these are typical in rst which

cpp/demo/biharmonic/main.cpp

@@ -191,9 +191,9 @@ int main(int argc, char* argv[])
 
     //  Define variational forms
     auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_biharmonic_a, {V, V}, {}, {{"alpha", alpha}}, {}));
+        *form_biharmonic_a, {V, V}, {}, {{"alpha", alpha}}, {}, {}));
     auto L = std::make_shared<fem::Form<T>>(
-        fem::create_form<T>(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}));
+        fem::create_form<T>(*form_biharmonic_L, {V}, {{"f", f}}, {}, {}, {}));
 
     //  Now, the Dirichlet boundary condition ($u = 0$) can be
     //  created using the class {cpp:class}`DirichletBC`. A

cpp/demo/codim_0_assembly/main.cpp

@@ -0,0 +1,173 @@
+// # Mixed assembly with a function mesh on a subset of cells
+//
+// This demo illustrates how to:
+//
+// * Create a submesh of co-dimension 0
+// * Assemble a mixed formulation with function spaces defined on the sub mesh
+// and parent mesh
+
+#include "mixed_codim0.h"
+#include <basix/finite-element.h>
+#include <cmath>
+#include <dolfinx.h>
+#include <dolfinx/fem/Constant.h>
+#include <dolfinx/fem/petsc.h>
+#include <dolfinx/la/MatrixCSR.h>
+#include <dolfinx/la/SparsityPattern.h>
+#include <utility>
+#include <vector>
+
+using namespace dolfinx;
+using T = PetscScalar;
+using U = typename dolfinx::scalar_value_type_t<T>;
+
+int main(int argc, char* argv[])
+{
+  dolfinx::init_logging(argc, argv);
+  PetscInitialize(&argc, &argv, nullptr, nullptr);
+
+  {
+    // Create mesh and function space
+    auto part = mesh::create_cell_partitioner(mesh::GhostMode::shared_facet);
+    auto mesh = std::make_shared<mesh::Mesh<U>>(
+        mesh::create_rectangle<U>(MPI_COMM_WORLD, {{{0.0, 0.0}, {2.0, 1.0}}},
+                                  {1, 4}, mesh::CellType::quadrilateral, part));
+
+    auto element = basix::create_element<U>(
+        basix::element::family::P, basix::cell::type::quadrilateral, 1,
+        basix::element::lagrange_variant::unset,
+        basix::element::dpc_variant::unset, false);
+
+    auto V = std::make_shared<fem::FunctionSpace<U>>(
+        fem::create_functionspace(mesh, element, {}));
+
+    // Next we find all cells of the mesh with y<0.5
+    const int tdim = mesh->topology()->dim();
+    auto marked_cells = mesh::locate_entities(
+        *mesh, tdim,
+        [](auto x)
+        {
+          using U = typename decltype(x)::value_type;
+          constexpr U eps = 1.0e-8;
+          std::vector<std::int8_t> marker(x.extent(1), false);
+          for (std::size_t p = 0; p < x.extent(1); ++p)
+          {
+            auto y = x(1, p);
+            if (std::abs(y) <= 0.5 + eps)
+              marker[p] = true;
+          }
+          return marker;
+        });
+    // We create a MeshTags object where we mark these cells with 2, and any
+    // other cell with 1
+    auto cell_map = mesh->topology()->index_map(tdim);
+    std::size_t num_cells_local
+        = mesh->topology()->index_map(tdim)->size_local()
+          + mesh->topology()->index_map(tdim)->num_ghosts();
+    std::vector<std::int32_t> cells(num_cells_local);
+    std::iota(cells.begin(), cells.end(), 0);
+    std::vector<std::int32_t> values(cell_map->size_local(), 1);
+    std::for_each(marked_cells.begin(), marked_cells.end(),
+                  [&values](auto& c) { values[c] = 2; });
+    dolfinx::mesh::MeshTags<std::int32_t> cell_marker(mesh->topology(), tdim,
+                                                      cells, values);
+
+    std::shared_ptr<mesh::Mesh<U>> submesh;
+    std::vector<std::int32_t> submesh_to_mesh;
+    {
+      auto [_submesh, _submesh_to_mesh, v_map, g_map]
+          = mesh::create_submesh(*mesh, tdim, cell_marker.find(2));
+      submesh = std::make_shared<mesh::Mesh<U>>(std::move(_submesh));
+      submesh_to_mesh = std::move(_submesh_to_mesh);
+    }
+
+    // We create the function space used for the trial space
+    auto W = std::make_shared<fem::FunctionSpace<U>>(
+        fem::create_functionspace(submesh, element, {}));
+
+    // A mixed-domain form has functions defined over different meshes. The mesh
+    // associated with the measure (dx, ds, etc.) is called the integration
+    // domain. To assemble mixed-domain forms, maps must be provided taking
+    // entities in the integration domain to entities on each mesh in the form.
+    // Since one of our forms has a measure defined over `mesh` and involves a
+    // function defined over `submesh`, we must provide a map from entities in
+    // `mesh` to entities in `submesh`. This is simply the "inverse" of
+    // `submesh_to_mesh`.
+    std::vector<std::int32_t> mesh_to_submesh(num_cells_local, -1);
+    for (std::size_t i = 0; i < submesh_to_mesh.size(); ++i)
+      mesh_to_submesh[submesh_to_mesh[i]] = i;
+
+    std::shared_ptr<const mesh::Mesh<U>> const_ptr = submesh;
+    std::map<std::shared_ptr<const mesh::Mesh<U>>,
+             std::span<const std::int32_t>>
+        entity_maps
+        = {{const_ptr, std::span<const std::int32_t>(mesh_to_submesh.data(),
+                                                     mesh_to_submesh.size())}};
+
+    // Next we compute the integration entities on the integration domain `mesh`
+    std::map<
+        fem::IntegralType,
+        std::vector<std::pair<std::int32_t, std::span<const std::int32_t>>>>
+        subdomain_map = {};
+    auto integration_entities = fem::compute_integration_domains(
+        fem::IntegralType::cell, *mesh->topology(), cell_marker.find(2), tdim);
+
+    subdomain_map[fem::IntegralType::cell].push_back(
+        {3, std::span<const std::int32_t>(integration_entities.data(),
+                                          integration_entities.size())});
+
+    // We can now create the bi-linear form
+    auto a_mixed = std::make_shared<fem::Form<T>>(
+        fem::create_form<T>(*form_mixed_codim0_a_mixed, {V, W}, {}, {},
+                            subdomain_map, entity_maps, V->mesh()));
+
+    la::SparsityPattern sp_mixed = fem::create_sparsity_pattern(*a_mixed);
+    sp_mixed.finalize();
+    la::MatrixCSR<double> A_mixed(sp_mixed);
+    fem::assemble_matrix(A_mixed.mat_add_values(), *a_mixed, {});
+    A_mixed.scatter_rev();
+
+    auto a = std::make_shared<fem::Form<T>>(
+        fem::create_form<T>(*form_mixed_codim0_a, {W, W}, {}, {}, {}, {}));
+    la::SparsityPattern sp = fem::create_sparsity_pattern(*a);
+    sp.finalize();
+
+    la::MatrixCSR<double> A(sp);
+    fem::assemble_matrix(A.mat_add_values(), *a, {});
+    A.scatter_rev();
+
+    std::vector<T> A_mixed_flattened = A_mixed.to_dense();
+    std::stringstream cc;
+    cc.precision(3);
+    cc << "A_mixed:" << std::endl;
+
+    std::size_t num_owned_rows = V->dofmap()->index_map->size_local();
+    std::size_t num_sub_cols = W->dofmap()->index_map->size_local()
+                               + W->dofmap()->index_map->num_ghosts();
+    for (std::size_t i = 0; i < num_owned_rows; i++)
+    {
+      for (std::size_t j = 0; j < num_sub_cols; j++)
+      {
+        cc << A_mixed_flattened[i * num_sub_cols + j] << " ";
+      }
+      cc << std::endl;
+    }
+
+    std::size_t num_owned_sub_rows = W->dofmap()->index_map->size_local();
+    std::vector<T> A_flattened = A.to_dense();
+    cc << "A" << std::endl;
+    for (std::size_t i = 0; i < num_owned_sub_rows; i++)
+    {
+      for (std::size_t j = 0; j < num_sub_cols; j++)
+      {
+        cc << A_flattened[i * num_sub_cols + j] << " ";
+      }
+      cc << std::endl;
+    }
+    std::cout << cc.str() << std::endl;
+  }
+
+  PetscFinalize();
+
+  return 0;
+}

cpp/demo/codim_0_assembly/mixed_codim0.py

@@ -0,0 +1,35 @@
+# This demo aims to illustrate how to assemble a matrix with a trial function
+# defined on a submesh of co-dimension 0, and a test function defined on the parent mesh
+from basix.ufl import element
+from ufl import (
+    FunctionSpace,
+    Mesh,
+    TestFunction,
+    TrialFunction,
+    dx,
+)
+
+cell = "quadrilateral"
+coord_element = element("Lagrange", cell, 1, shape=(2,))
+mesh = Mesh(coord_element)
+
+# We define the function space and test function on the full mesh
+e = element("Lagrange", cell, 1)
+V = FunctionSpace(mesh, e)
+v = TestFunction(V)
+
+# Next we define the sub-mesh
+submesh = Mesh(coord_element)
+W = FunctionSpace(submesh, e)
+p = TrialFunction(W)
+
+# And finally we define a "mass matrix" on the submesh, with the test function
+# of the parent mesh. The integration domain is the parent mesh, but we restrict integration
+# to all cells marked with subdomain_id=3, which will indicate what cells of our mesh is part
+# of the submesh
+a_mixed = p * v * dx(domain=mesh, subdomain_id=3)
+
+q = TestFunction(W)
+a = p * q * dx(domain=submesh)
+
+forms = [a_mixed, a]

cpp/demo/hyperelasticity/CMakeLists.txt

@@ -40,9 +40,7 @@ else()
 
   set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
-  add_executable(
-    ${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c
-  )
+  add_executable(${PROJECT_NAME} main.cpp ${CMAKE_CURRENT_BINARY_DIR}/hyperelasticity.c)
   target_link_libraries(${PROJECT_NAME} dolfinx)
 
   # Do not throw error for 'multi-line comments' (these are typical in rst which

cpp/demo/hyperelasticity/main.cpp

@@ -175,10 +175,10 @@ int main(int argc, char* argv[])
     auto u = std::make_shared<fem::Function<T>>(V);
     auto a = std::make_shared<fem::Form<T>>(
         fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{"u", u}},
-                            {{"B", B}, {"T", traction}}, {}));
+                            {{"B", B}, {"T", traction}}, {}, {}));
     auto L = std::make_shared<fem::Form<T>>(
         fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{"u", u}},
-                            {{"B", B}, {"T", traction}}, {}));
+                            {{"B", B}, {"T", traction}}, {}, {}));
 
     auto u_rotation = std::make_shared<fem::Function<T>>(V);
     u_rotation->interpolate(

cpp/demo/poisson/main.cpp

@@ -140,9 +140,9 @@ int main(int argc, char* argv[])
 
     // Define variational forms
     auto a = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_poisson_a, {V, V}, {}, {{"kappa", kappa}}, {}));
+        *form_poisson_a, {V, V}, {}, {{"kappa", kappa}}, {}, {}));
     auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
-        *form_poisson_L, {V}, {{"f", f}, {"g", g}}, {}, {}));
+        *form_poisson_L, {V}, {{"f", f}, {"g", g}}, {}, {}, {}));
 
     //  Now, the Dirichlet boundary condition ($u = 0$) can be created
     //  using the class {cpp:class}`DirichletBC`. A

cpp/demo/poisson_matrix_free/main.cpp

@@ -151,12 +151,12 @@ void solver(MPI_Comm comm)
 
   // Define variational forms
   auto L = std::make_shared<fem::Form<T, U>>(
-      fem::create_form<T>(*form_poisson_L, {V}, {}, {{"f", f}}, {}));
+      fem::create_form<T>(*form_poisson_L, {V}, {}, {{"f", f}}, {}, {}));
 
   // Action of the bilinear form "a" on a function ui
   auto ui = std::make_shared<fem::Function<T, U>>(V);
   auto M = std::make_shared<fem::Form<T, U>>(
-      fem::create_form<T>(*form_poisson_M, {V}, {{"ui", ui}}, {{}}, {}));
+      fem::create_form<T>(*form_poisson_M, {V}, {{"ui", ui}}, {{}}, {}, {}));
 
   // Define boundary condition
   auto u_D = std::make_shared<fem::Function<T, U>>(V);
@@ -231,7 +231,7 @@ void solver(MPI_Comm comm)
   // Compute L2 error (squared) of the solution vector e = (u - u_d, u
   // - u_d)*dx
   auto E = std::make_shared<fem::Form<T>>(fem::create_form<T, U>(
-      *form_poisson_E, {}, {{"uexact", u_D}, {"usol", u}}, {}, {}, mesh));
+      *form_poisson_E, {}, {{"uexact", u_D}, {"usol", u}}, {}, {}, {}, mesh));
   T error = fem::assemble_scalar(*E);
   if (dolfinx::MPI::rank(comm) == 0)
   {

cpp/dolfinx/fem/DofMap.cpp

@@ -78,7 +78,8 @@ fem::DofMap build_collapsed_dofmap(const DofMap& dofmap_view,
   }
 
   // Create a map from old dofs to new dofs
-  std::vector<std::int32_t> old_to_new(dofs_view.back() + bs_view, -1);
+  std::size_t array_size = dofs_view.empty() ? 0 : dofs_view.back() + bs_view;
+  std::vector<std::int32_t> old_to_new(array_size, -1);
   for (std::size_t new_idx = 0; new_idx < sub_imap_to_imap.size(); ++new_idx)
   {
     for (int k = 0; k < bs_view; ++k)
@@ -215,7 +216,6 @@ std::pair<DofMap, std::vector<std::int32_t>> DofMap::collapse(
     reorder_fn = [](const graph::AdjacencyList<std::int32_t>& g)
     { return graph::reorder_gps(g); };
   }
-
   // Create new dofmap
   auto create_subdofmap = [](MPI_Comm comm, auto index_map_bs, auto& layout,
                              auto& topology, auto& reorder_fn, auto& dmap)
@@ -227,7 +227,6 @@ 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, dofmaps] = build_dofmap_data(
           comm, topology, {collapsed_dof_layout}, reorder_fn);
       auto index_map

cpp/dolfinx/fem/assemble_vector_impl.h

@@ -1061,8 +1061,6 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
 /// @param[in,out] b The vector to be modified
 /// @param[in] a The bilinear forms, where a[j] is the form that
 /// generates A_j
-/// @param[in] x_dofmap Mesh geometry dofmap
-/// @param[in] x Mesh coordinates
 /// @param[in] constants Constants that appear in `a`
 /// @param[in] coeffs Coefficients that appear in `a`
 /// @param[in] bcs1 List of boundary conditions for each block, i.e.
@@ -1073,15 +1071,13 @@ void lift_bc(std::span<T> b, const Form<T, U>& a, mdspan2_t x_dofmap,
 template <dolfinx::scalar T, std::floating_point U>
 void apply_lifting(
     std::span<T> b, const std::vector<std::shared_ptr<const Form<T, U>>> a,
-    mdspan2_t x_dofmap, std::span<const scalar_value_type_t<T>> x,
     const std::vector<std::span<const T>>& constants,
     const std::vector<std::map<std::pair<IntegralType, int>,
                                std::pair<std::span<const T>, int>>>& coeffs,
     const std::vector<std::vector<std::shared_ptr<const DirichletBC<T, U>>>>&
         bcs1,
     const std::vector<std::span<const T>>& x0, T scale)
 {
-  // FIXME: make changes to reactivate this check
   if (!x0.empty() and x0.size() != a.size())
   {
     throw std::runtime_error(
@@ -1100,6 +1096,13 @@ void apply_lifting(
     std::vector<T> bc_values1;
     if (a[j] and !bcs1[j].empty())
     {
+      // Extract data from mesh
+      std::shared_ptr<const mesh::Mesh<U>> mesh = a[j]->mesh();
+      if (!mesh)
+        throw std::runtime_error("Unable to extract a mesh.");
+      mdspan2_t x_dofmap = mesh->geometry().dofmap();
+      auto x = mesh->geometry().x();
+
       assert(a[j]->function_spaces().at(0));
 
       auto V1 = a[j]->function_spaces()[1];

cpp/dolfinx/fem/assembler.h

@@ -161,31 +161,7 @@ void apply_lifting(
   if (std::ranges::all_of(a, [](auto ptr) { return ptr == nullptr; }))
     return;
 
-  std::shared_ptr<const mesh::Mesh<U>> mesh;
-  for (auto& a_i : a)
-  {
-    if (a_i and !mesh)
-      mesh = a_i->mesh();
-    if (a_i and mesh and a_i->mesh() != mesh)
-      throw std::runtime_error("Mismatch between meshes.");
-  }
-
-  if (!mesh)
-    throw std::runtime_error("Unable to extract a mesh.");
-
-  if constexpr (std::is_same_v<U, scalar_value_type_t<T>>)
-  {
-    impl::apply_lifting<T>(b, a, mesh->geometry().dofmap(),
-                           mesh->geometry().x(), constants, coeffs, bcs1, x0,
-                           scale);
-  }
-  else
-  {
-    auto x = mesh->geometry().x();
-    std::vector<scalar_value_type_t<T>> _x(x.begin(), x.end());
-    impl::apply_lifting<T>(b, a, mesh->geometry().dofmap(), _x, constants,
-                           coeffs, bcs1, x0, scale);
-  }
+  impl::apply_lifting<T>(b, a, constants, coeffs, bcs1, x0, scale);
 }
 
 /// Modify b such that: