Merge branch 'main' into jpdean/mixed_mesh_hackathon_2

cpp/dolfinx/fem/CoordinateElement.h

@@ -78,7 +78,9 @@ class CoordinateElement
   /// compute. Use 0 for the basis functions only
   /// @param[in] num_points Number of points at which to evaluate the
   /// basis functions.
-  /// @return Shape of the array to be filled by `tabulate`.
+  /// @return Shape of the array to be filled by `tabulate`, where (0)
+  /// is derivative index, (1) is the point index, (2) is the basis
+  /// function index and (3) is the basis function component.
   std::array<std::size_t, 4> tabulate_shape(std::size_t nd,
                                             std::size_t num_points) const;
 

docker/Dockerfile.test-env

@@ -13,21 +13,21 @@
 #
 
 ARG ADIOS2_VERSION=2.10.2
-ARG DOXYGEN_VERSION=1_12_0
+ARG DOXYGEN_VERSION=1_13_1
 ARG GMSH_VERSION=4_13_1
 ARG HDF5_VERSION=1.14.5
-ARG KAHIP_VERSION=3.16
+ARG KAHIP_VERSION=3.17
 # NOTE: The NumPy version (https://pypi.org/project/numpy/#history)
 # should be pinned to the most recent NumPy release that is supported by
 # the most recent Numba release, see
 # https://numba.readthedocs.io/en/stable/user/installing.html#version-support-information
 ARG NUMPY_VERSION=2.0.2
-ARG PETSC_VERSION=3.22.1
-ARG SLEPC_VERSION=3.22.1
+ARG PETSC_VERSION=3.22.2
+ARG SLEPC_VERSION=3.22.2
 
 ARG MPICH_VERSION=4.2.3
 ARG OPENMPI_SERIES=5.0
-ARG OPENMPI_PATCH=5
+ARG OPENMPI_PATCH=6
 
 ########################################
 

python/dolfinx/fem/__init__.py

@@ -13,6 +13,7 @@
 from dolfinx.cpp.fem import compute_integration_domains as _compute_integration_domains
 from dolfinx.cpp.fem import create_interpolation_data as _create_interpolation_data
 from dolfinx.cpp.fem import create_sparsity_pattern as _create_sparsity_pattern
+from dolfinx.cpp.fem import discrete_curl as _discrete_curl
 from dolfinx.cpp.fem import discrete_gradient as _discrete_gradient
 from dolfinx.cpp.fem import interpolation_matrix as _interpolation_matrix
 from dolfinx.cpp.la import SparsityPattern
@@ -120,21 +121,38 @@ def create_interpolation_data(
     )
 
 
+def discrete_curl(V0: FunctionSpace, V1: FunctionSpace) -> _MatrixCSR:
+    """Assemble a discrete curl operator.
+
+    The discrete curl operator interpolates the curl of H(curl) finite
+    element function into a H(div) space.
+
+    Args:
+        V0: H1(curl) space to interpolate the curl from.
+        V1: H(div) space to interpolate into.
+
+    Returns:
+        Discrete curl operator.
+    """
+    return _MatrixCSR(_discrete_curl(V0._cpp_object, V1._cpp_object))
+
+
 def discrete_gradient(space0: FunctionSpace, space1: FunctionSpace) -> _MatrixCSR:
     """Assemble a discrete gradient operator.
 
-    The discrete gradient operator interpolates the gradient of
-    a H1 finite element function into a H(curl) space. It is assumed that
-    the H1 space uses an identity map and the H(curl) space uses a covariant Piola map.
+    The discrete gradient operator interpolates the gradient of a H1
+    finite element function into a H(curl) space. It is assumed that the
+    H1 space uses an identity map and the H(curl) space uses a covariant
+    Piola map.
 
     Args:
-        space0: H1 space to interpolate the gradient from
-        space1: H(curl) space to interpolate into
+        space0: H1 space to interpolate the gradient from.
+        space1: H(curl) space to interpolate into.
 
     Returns:
-        Discrete gradient operator
+        Discrete gradient operator.
     """
-    return _discrete_gradient(space0._cpp_object, space1._cpp_object)
+    return _MatrixCSR(_discrete_gradient(space0._cpp_object, space1._cpp_object))
 
 
 def interpolation_matrix(space0: FunctionSpace, space1: FunctionSpace) -> _MatrixCSR:
@@ -155,13 +173,14 @@ def compute_integration_domains(
 ):
     """Given an integral type and a set of entities compute integration entities.
 
-    This function returns a list `[(id, entities)]`. For cell integrals
-    `entities` are the cell indices. For exterior facet integrals,
-    `entities` is a list of `(cell_index, local_facet_index)` pairs. For
-    interior facet integrals, `entities` is a list of `(cell_index0,
-    local_facet_index0, cell_index1, local_facet_index1)`. `id` refers
-    to the subdomain id used in the definition of the integration
-    measures of the variational form.
+    This function returns a list ``[(id, entities)]``. For cell
+    integrals ``entities`` are the cell indices. For exterior facet
+    integrals, ``entities`` is a list of ``(cell_index,
+    local_facet_index)`` pairs. For interior facet integrals,
+    ``entities`` is a list of ``(cell_index0, local_facet_index0,
+    cell_index1, local_facet_index1)``. ``id`` refers to the subdomain
+    id used in the definition of the integration measures of the
+    variational form.
 
     Note:
         Owned mesh entities only are returned. Ghost entities are not

python/dolfinx/fem/petsc.py

@@ -31,6 +31,7 @@
 from dolfinx import la
 from dolfinx.cpp.fem import pack_coefficients as _pack_coefficients
 from dolfinx.cpp.fem import pack_constants as _pack_constants
+from dolfinx.cpp.fem.petsc import discrete_curl as _discrete_curl
 from dolfinx.cpp.fem.petsc import discrete_gradient as _discrete_gradient
 from dolfinx.cpp.fem.petsc import interpolation_matrix as _interpolation_matrix
 from dolfinx.fem import assemble as _assemble
@@ -60,6 +61,7 @@
     "create_vector",
     "create_vector_block",
     "create_vector_nest",
+    "discrete_curl",
     "discrete_gradient",
     "interpolation_matrix",
     "set_bc",
@@ -1043,6 +1045,19 @@ def J(self, x: PETSc.Vec, A: PETSc.Mat) -> None:
         A.assemble()
 
 
+def discrete_curl(space0: _FunctionSpace, space1: _FunctionSpace) -> PETSc.Mat:
+    """Assemble a discrete curl operator.
+
+    Args:
+        space0: H1 space to interpolate the gradient from.
+        space1: H(curl) space to interpolate into.
+
+    Returns:
+        Discrete curl operator.
+    """
+    return _discrete_curl(space0._cpp_object, space1._cpp_object)
+
+
 def discrete_gradient(space0: _FunctionSpace, space1: _FunctionSpace) -> PETSc.Mat:
     """Assemble a discrete gradient operator.
 

python/dolfinx/wrappers/assemble.cpp

@@ -139,6 +139,20 @@ void declare_discrete_operators(nb::module_& m)
           return A;
         });
 
+  m.def(
+      "discrete_curl",
+      [](const dolfinx::fem::FunctionSpace<U>& V0,
+         const dolfinx::fem::FunctionSpace<U>& V1)
+      {
+        dolfinx::la::SparsityPattern sp = create_sparsity(V0, V1);
+
+        // Build operator
+        dolfinx::la::MatrixCSR<T> A(sp);
+        dolfinx::fem::discrete_curl<U, T>(V0, V1, A.mat_set_values());
+        return A;
+      },
+      nb::arg("V0"), nb::arg("V1"));
+
   m.def(
       "discrete_gradient",
       [](const dolfinx::fem::FunctionSpace<U>& V0,

python/dolfinx/wrappers/petsc.cpp

@@ -1,4 +1,4 @@
-// Copyright (C) 2017-2023 Chris Richardson and Garth N. Wells
+// Copyright (C) 2017-2025 Chris Richardson and Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -10,6 +10,7 @@
 #include "caster_mpi.h"
 #include "caster_petsc.h"
 #include "pycoeff.h"
+#include <concepts>
 #include <dolfinx/common/IndexMap.h>
 #include <dolfinx/fem/DirichletBC.h>
 #include <dolfinx/fem/DofMap.h>
@@ -41,19 +42,58 @@
 namespace
 {
 // Declare assembler function that have multiple scalar types
-template <typename T, typename U>
+template <typename T, std::floating_point U>
 void declare_petsc_discrete_operators(nb::module_& m)
 {
   m.def(
-      "discrete_gradient",
+      "discrete_curl",
       [](const dolfinx::fem::FunctionSpace<U>& V0,
          const dolfinx::fem::FunctionSpace<U>& V1)
       {
         assert(V0.mesh());
         auto mesh = V0.mesh();
         assert(V1.mesh());
         assert(mesh == V1.mesh());
+
+        auto dofmap0 = V0.dofmap();
+        assert(dofmap0);
+        auto dofmap1 = V1.dofmap();
+        assert(dofmap1);
+
+        // Create and build  sparsity pattern
+        assert(dofmap0->index_map);
+        assert(dofmap1->index_map);
         MPI_Comm comm = mesh->comm();
+        dolfinx::la::SparsityPattern sp(
+            comm, {dofmap1->index_map, dofmap0->index_map},
+            {dofmap1->index_map_bs(), dofmap0->index_map_bs()});
+
+        int tdim = mesh->topology()->dim();
+        auto map = mesh->topology()->index_map(tdim);
+        assert(map);
+        std::vector<std::int32_t> c(map->size_local(), 0);
+        std::iota(c.begin(), c.end(), 0);
+        dolfinx::fem::sparsitybuild::cells(sp, {c, c}, {*dofmap1, *dofmap0});
+        sp.finalize();
+
+        // Build operator
+        Mat A = dolfinx::la::petsc::create_matrix(comm, sp);
+        MatSetOption(A, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
+        dolfinx::fem::discrete_curl<U, T>(
+            V0, V1, dolfinx::la::petsc::Matrix::set_fn(A, INSERT_VALUES));
+        return A;
+      },
+      nb::rv_policy::take_ownership, nb::arg("V0"), nb::arg("V1"));
+
+  m.def(
+      "discrete_gradient",
+      [](const dolfinx::fem::FunctionSpace<U>& V0,
+         const dolfinx::fem::FunctionSpace<U>& V1)
+      {
+        assert(V0.mesh());
+        auto mesh = V0.mesh();
+        assert(V1.mesh());
+        assert(mesh == V1.mesh());
 
         auto dofmap0 = V0.dofmap();
         assert(dofmap0);
@@ -63,6 +103,7 @@ void declare_petsc_discrete_operators(nb::module_& m)
         // Create and build  sparsity pattern
         assert(dofmap0->index_map);
         assert(dofmap1->index_map);
+        MPI_Comm comm = mesh->comm();
         dolfinx::la::SparsityPattern sp(
             comm, {dofmap1->index_map, dofmap0->index_map},
             {dofmap1->index_map_bs(), dofmap0->index_map_bs()});
@@ -94,7 +135,6 @@ void declare_petsc_discrete_operators(nb::module_& m)
         auto mesh = V0.mesh();
         assert(V1.mesh());
         assert(mesh == V1.mesh());
-        MPI_Comm comm = mesh->comm();
 
         auto dofmap0 = V0.dofmap();
         assert(dofmap0);
@@ -104,6 +144,7 @@ void declare_petsc_discrete_operators(nb::module_& m)
         // Create and build  sparsity pattern
         assert(dofmap0->index_map);
         assert(dofmap1->index_map);
+        MPI_Comm comm = mesh->comm();
         dolfinx::la::SparsityPattern sp(
             comm, {dofmap1->index_map, dofmap0->index_map},
             {dofmap1->index_map_bs(), dofmap0->index_map_bs()});