Extend geometry scalar typing (#2611)

* Work on Python geometry type * Type updates * Type updates * Type update * Simplifications * Set x data layout * Add concepts constrain on function * More mesh float typing * Tidy up * Remove duplications * Interface fix * Doc improvements * Template fides over type * Fix typo * Doc fix * Doc update * Doc fixes * Test fix and extension * Add geometry type test * Remove function * Fix test * Range fix * Verbose tests * Fix test bug * Undo CI change * Simplifications * Re-enable test * Expand Fides test * Fix typo
FEniCS · Apr 4, 2023 · 4fb398f · 4fb398f
1 parent ea67cef
commit 4fb398f
Show file tree

Hide file tree

Showing 26 changed files with 1,109 additions and 1,343 deletions.
diff --git a/cpp/doc/source/fem.rst b/cpp/doc/source/fem.rst
@@ -93,9 +93,6 @@ Sparsity pattern construction
 .. doxygenfunction:: dolfinx::fem::create_sparsity_pattern(const Form<T, double>&)
    :project: DOLFINx
 
-.. doxygenfunction:: dolfinx::fem::create_sparsity_pattern(const mesh::Topology &topology, const std::array<std::reference_wrapper<const DofMap>, 2> &dofmaps, const std::set<IntegralType> &integrals)
-   :project: DOLFINx
-
 
 PETSc helpers
 -------------

diff --git a/cpp/dolfinx/fem/Form.h b/cpp/dolfinx/fem/Form.h
@@ -65,22 +65,25 @@ template <typename T, std::floating_point U = dolfinx::scalar_value_type_t<T>>
 class Form
 {
 public:
+  /// Scalar type
+  using scalar_type = T;
+
   /// @brief Create a finite element form.
   ///
-  /// @note User applications will normally call a fem::Form builder
-  /// function rather using this interface directly.
+  /// @note User applications will normally call a builder function
+  /// rather using this interface directly.
   ///
   /// @param[in] V Function spaces for the form arguments
-  /// @param[in] integrals The integrals in the form. The first key is
-  /// the domain type. For each key there is a list of tuples (domain id,
+  /// @param[in] integrals Integrals in the form. The first key is the
+  /// domain type. For each key there is a list of tuples (domain id,
   /// integration kernel, entities).
   /// @param[in] coefficients
   /// @param[in] constants Constants in the Form
   /// @param[in] needs_facet_permutations Set to true is any of the
   /// integration kernels require cell permutation data
-  /// @param[in] mesh The mesh of the domain. This is required when
-  /// there are not argument functions from which the mesh can be
-  /// extracted, e.g. for functionals
+  /// @param[in] mesh Mesh of the domain. This is required when there
+  /// are no argument functions from which the mesh can be extracted,
+  /// e.g. for functionals.
   Form(const std::vector<std::shared_ptr<const FunctionSpace<U>>>& V,
        const std::map<IntegralType,
                       std::vector<std::tuple<
@@ -325,9 +328,6 @@ class Form
     return _constants;
   }
 
-  /// Scalar type (T)
-  using scalar_type = T;
-
 private:
   using kern = std::function<void(T*, const T*, const T*,
                                   const scalar_value_type_t<T>*, const int*,

diff --git a/cpp/dolfinx/fem/assembler.h b/cpp/dolfinx/fem/assembler.h
@@ -243,7 +243,8 @@ void apply_lifting(
 /// local index.
 template <typename T, std::floating_point U>
 void assemble_matrix(
-    auto mat_add, const Form<T, U>& a, std::span<const T> constants,
+    la::MatSet<T> auto mat_add, const Form<T, U>& a,
+    std::span<const T> constants,
     const std::map<std::pair<IntegralType, int>,
                    std::pair<std::span<const T>, int>>& coefficients,
     std::span<const std::int8_t> dof_marker0,

diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h
@@ -72,7 +72,9 @@ Mat create_matrix_block(
 
   std::shared_ptr mesh = V[0][0]->mesh();
   assert(mesh);
-  const int tdim = mesh->topology()->dim();
+  auto topology = mesh->topology();
+  assert(topology);
+  const int tdim = topology->dim();
 
   // Build sparsity pattern for each block
   std::vector<std::vector<std::unique_ptr<la::SparsityPattern>>> patterns(
@@ -94,22 +96,58 @@ Mat create_matrix_block(
         // Build sparsity pattern for block
         assert(V[0][row]->dofmap());
         assert(V[1][col]->dofmap());
-        std::array<const std::reference_wrapper<const DofMap>, 2> dofmaps{
+        std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
             *V[0][row]->dofmap(), *V[1][col]->dofmap()};
         assert(patterns[row].back());
         auto& sp = patterns[row].back();
         assert(sp);
+
         if (form->num_integrals(IntegralType::cell) > 0)
-          sparsitybuild::cells(*sp, *mesh->topology(), dofmaps);
+        {
+          auto map = topology->index_map(tdim);
+          assert(map);
+          std::vector<std::int32_t> c(map->size_local(), 0);
+          std::iota(c.begin(), c.end(), 0);
+          sparsitybuild::cells(*sp, c, dofmaps);
+        }
+
         if (form->num_integrals(IntegralType::interior_facet) > 0)
         {
+          // Loop over owned facets
           mesh->topology_mutable()->create_entities(tdim - 1);
-          sparsitybuild::interior_facets(*sp, *mesh->topology(), dofmaps);
+          auto f_to_c = topology->connectivity(tdim - 1, tdim);
+          if (!f_to_c)
+          {
+            throw std::runtime_error(
+                "Facet-cell connectivity has not been computed.");
+          }
+          auto map = topology->index_map(tdim - 1);
+          assert(map);
+          std::vector<std::int32_t> facets;
+          facets.reserve(2 * map->size_local());
+          for (int f = 0; f < map->size_local(); ++f)
+            if (auto cells = f_to_c->links(f); cells.size() == 2)
+              facets.insert(facets.end(), {cells[0], cells[1]});
+          sparsitybuild::interior_facets(*sp, facets, dofmaps);
         }
+
         if (form->num_integrals(IntegralType::exterior_facet) > 0)
         {
+          // Loop over owned facets
           mesh->topology_mutable()->create_entities(tdim - 1);
-          sparsitybuild::exterior_facets(*sp, *mesh->topology(), dofmaps);
+          auto connectivity = topology->connectivity(tdim - 1, tdim);
+          if (!connectivity)
+          {
+            throw std::runtime_error(
+                "Facet-cell connectivity has not been computed.");
+          }
+          auto map = topology->index_map(tdim - 1);
+          assert(map);
+          std::vector<std::int32_t> cells;
+          for (int f = 0; f < map->size_local(); ++f)
+            if (auto c = connectivity->links(f); c.size() == 1)
+              cells.push_back(c[0]);
+          sparsitybuild::cells(*sp, cells, dofmaps);
         }
       }
       else

diff --git a/cpp/dolfinx/fem/sparsitybuild.cpp b/cpp/dolfinx/fem/sparsitybuild.cpp
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019 Garth N. Wells
+// Copyright (C) 2007-2023 Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -16,85 +16,24 @@ using namespace dolfinx::fem;
 
 //-----------------------------------------------------------------------------
 void sparsitybuild::cells(
-    la::SparsityPattern& pattern, const mesh::Topology& topology,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps)
+    la::SparsityPattern& pattern, std::span<const std::int32_t> cells,
+    std::array<std::reference_wrapper<const DofMap>, 2> dofmaps)
 {
-  const int D = topology.dim();
-  auto cells = topology.connectivity(D, 0);
-  assert(cells);
-  for (int c = 0; c < cells->num_nodes(); ++c)
-  {
-    pattern.insert(dofmaps[0].get().cell_dofs(c),
-                   dofmaps[1].get().cell_dofs(c));
-  }
-}
-//-----------------------------------------------------------------------------
-void sparsitybuild::cells(
-    la::SparsityPattern& pattern, const std::span<const std::int32_t>& cells,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps)
-{
-  for (std::int32_t c : cells)
-  {
-    pattern.insert(dofmaps[0].get().cell_dofs(c),
-                   dofmaps[1].get().cell_dofs(c));
-  }
-}
-//-----------------------------------------------------------------------------
-void sparsitybuild::interior_facets(
-    la::SparsityPattern& pattern, const mesh::Topology& topology,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps)
-{
-  const int D = topology.dim();
-  if (!topology.connectivity(D - 1, 0))
-    throw std::runtime_error("Topology facets have not been created.");
-
-  auto connectivity = topology.connectivity(D - 1, D);
-  if (!connectivity)
-    throw std::runtime_error("Facet-cell connectivity has not been computed.");
-
-  // Array to store macro-dofs, if required (for interior facets)
-  std::array<std::vector<std::int32_t>, 2> macro_dofs;
-
-  // Loop over owned facets
-  auto map = topology.index_map(D - 1);
-  assert(map);
-  const std::int32_t num_facets = map->size_local();
-  for (int f = 0; f < num_facets; ++f)
-  {
-    // Get cells incident with facet
-    auto cells = connectivity->links(f);
-
-    // Proceed to next facet if only connection
-    if (cells.size() == 1)
-      continue;
-
-    // Tabulate dofs for each dimension on macro element
-    assert(cells.size() == 2);
-    const int cell0 = cells[0];
-    const int cell1 = cells[1];
-    for (std::size_t i = 0; i < 2; i++)
-    {
-      auto cell_dofs0 = dofmaps[i].get().cell_dofs(cell0);
-      auto cell_dofs1 = dofmaps[i].get().cell_dofs(cell1);
-      macro_dofs[i].resize(cell_dofs0.size() + cell_dofs1.size());
-      std::copy(cell_dofs0.begin(), cell_dofs0.end(), macro_dofs[i].begin());
-      std::copy(cell_dofs1.begin(), cell_dofs1.end(),
-                std::next(macro_dofs[i].begin(), cell_dofs0.size()));
-    }
-
-    pattern.insert(macro_dofs[0], macro_dofs[1]);
-  }
+  const DofMap& map0 = dofmaps[0].get();
+  const DofMap& map1 = dofmaps[1].get();
+  for (auto c : cells)
+    pattern.insert(map0.cell_dofs(c), map1.cell_dofs(c));
 }
 //-----------------------------------------------------------------------------
 void sparsitybuild::interior_facets(
-    la::SparsityPattern& pattern, const std::span<const std::int32_t>& facets,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps)
+    la::SparsityPattern& pattern, std::span<const std::int32_t> facets,
+    std::array<std::reference_wrapper<const DofMap>, 2> dofmaps)
 {
   std::array<std::vector<std::int32_t>, 2> macro_dofs;
-  for (std::size_t index = 0; index < facets.size(); index += 4)
+  for (std::size_t index = 0; index < facets.size(); index += 2)
   {
-    const int cell_0 = facets[index];
-    const int cell_1 = facets[index + 2];
+    int cell_0 = facets[index];
+    int cell_1 = facets[index + 1];
     for (std::size_t i = 0; i < 2; ++i)
     {
       auto cell_dofs_0 = dofmaps[i].get().cell_dofs(cell_0);
@@ -104,49 +43,7 @@ void sparsitybuild::interior_facets(
       std::copy(cell_dofs_1.begin(), cell_dofs_1.end(),
                 std::next(macro_dofs[i].begin(), cell_dofs_0.size()));
     }
-
     pattern.insert(macro_dofs[0], macro_dofs[1]);
   }
 }
 //-----------------------------------------------------------------------------
-void sparsitybuild::exterior_facets(
-    la::SparsityPattern& pattern, const mesh::Topology& topology,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps)
-{
-  const int D = topology.dim();
-  if (!topology.connectivity(D - 1, 0))
-    throw std::runtime_error("Topology facets have not been created.");
-
-  auto connectivity = topology.connectivity(D - 1, D);
-  if (!connectivity)
-    throw std::runtime_error("Facet-cell connectivity has not been computed.");
-
-  // Loop over owned facets
-  auto map = topology.index_map(D - 1);
-  assert(map);
-  const std::int32_t num_facets = map->size_local();
-  for (int f = 0; f < num_facets; ++f)
-  {
-    // Proceed to next facet if we have an interior facet
-    if (connectivity->num_links(f) == 2)
-      continue;
-
-    auto cells = connectivity->links(f);
-    assert(cells.size() == 1);
-    pattern.insert(dofmaps[0].get().cell_dofs(cells[0]),
-                   dofmaps[1].get().cell_dofs(cells[0]));
-  }
-}
-//-----------------------------------------------------------------------------
-void sparsitybuild::exterior_facets(
-    la::SparsityPattern& pattern, const std::span<const std::int32_t>& facets,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps)
-{
-  for (std::size_t index = 0; index < facets.size(); index += 2)
-  {
-    std::int32_t cell = facets[index];
-    pattern.insert(dofmaps[0].get().cell_dofs(cell),
-                   dofmaps[1].get().cell_dofs(cell));
-  }
-}
-//-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/fem/sparsitybuild.h b/cpp/dolfinx/fem/sparsitybuild.h
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019 Garth N. Wells
+// Copyright (C) 2007-2023 Garth N. Wells
 //
 // This file is part of DOLFINx (https://www.fenicsproject.org)
 //
@@ -28,80 +28,26 @@ class DofMap;
 /// Support for building sparsity patterns from degree-of-freedom maps.
 namespace sparsitybuild
 {
-
-/// @brief Iterate over cells and insert entries into sparsity pattern
-///
-/// @param[in,out] pattern The sparsity pattern to insert into
-/// @param[in] topology The mesh topology to build the sparsity pattern
-/// over
-/// @param[in] dofmaps The dofmap to use in building the sparsity
-/// pattern
-/// @note The sparsity pattern is not finalised
-void cells(
-    la::SparsityPattern& pattern, const mesh::Topology& topology,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps);
-
-/// @brief Iterate over cells and insert entries into sparsity pattern
+/// @brief Iterate over cells and insert entries into sparsity pattern.
 ///
 /// @param[in,out] pattern The sparsity pattern to insert into
 /// @param[in] cells The cell indices
-/// @param[in] dofmaps The dofmap to use in building the sparsity
-/// pattern
+/// @param[in] dofmaps Dofmaps to used in building the sparsity pattern
 /// @note The sparsity pattern is not finalised
-void cells(
-    la::SparsityPattern& pattern, const std::span<const std::int32_t>& cells,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps);
+void cells(la::SparsityPattern& pattern, std::span<const std::int32_t> cells,
+           std::array<std::reference_wrapper<const DofMap>, 2> dofmaps);
 
 /// @brief Iterate over interior facets and insert entries into sparsity
 /// pattern.
 ///
-/// @param[in,out] pattern The sparsity pattern to insert into
-/// @param[in] topology The mesh topology to build the sparsity pattern
-/// over
+/// @param[in,out] pattern Sparsity pattern to insert into
+/// @param[in] facets Facets as `(cell0, cell1)` pairs for each facet.
 /// @param[in] dofmaps The dofmap to use in building the sparsity
-/// pattern
-/// @note The sparsity pattern is not finalised
-void interior_facets(
-    la::SparsityPattern& pattern, const mesh::Topology& topology,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps);
-
-/// @brief Iterate over interior facets and insert entries into sparsity
 /// pattern.
-///
-/// @param[in,out] pattern The sparsity pattern to insert into
-/// @param[in] facets The facets as `(cell, local_index)` pairs, where
-/// `local_index` is the index of the facet relative to the `cell`
-/// @param[in] dofmaps The dofmap to use in building the sparsity
-/// pattern
-/// @note The sparsity pattern is not finalised
+/// @note The sparsity pattern is not finalised.
 void interior_facets(
-    la::SparsityPattern& pattern, const std::span<const std::int32_t>& facets,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps);
-
-/// @brief Iterate over exterior facets and insert entries into sparsity
-/// pattern.
-///
-/// @param[in,out] pattern The sparsity pattern to insert into
-/// @param[in] topology The mesh topology to build the sparsity pattern
-/// over
-/// @param[in] dofmaps The dofmap to use in building the sparsity
-/// pattern
-/// @note The sparsity pattern is not finalised
-void exterior_facets(
-    la::SparsityPattern& pattern, const mesh::Topology& topology,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps);
-
-/// @brief Iterate over exterior facets and insert entries into sparsity
-/// pattern.
-///
-/// @param[in,out] pattern The sparsity pattern to insert into
-/// @param[in] facets The facets as (cell, local_index) pairs
-/// @param[in] dofmaps The dofmap to use in building the sparsity
-/// pattern
-/// @note The sparsity pattern is not finalised
-void exterior_facets(
-    la::SparsityPattern& pattern, const std::span<const std::int32_t>& facets,
-    const std::array<const std::reference_wrapper<const DofMap>, 2>& dofmaps);
+    la::SparsityPattern& pattern, std::span<const std::int32_t> facets,
+    std::array<std::reference_wrapper<const DofMap>, 2> dofmaps);
 
 } // namespace sparsitybuild
 } // namespace dolfinx::fem