assembler.h

Go to the documentation of this file.

// Copyright (C) 2018-2026 Garth N. Wells and Jørgen S. Dokken
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include "Function.h"
#include "FunctionSpace.h"
#include "assemble_expression_impl.h"
#include "assemble_matrix_impl.h"
#include "assemble_scalar_impl.h"
#include "assemble_vector_impl.h"
#include "pack.h"
#include "traits.h"
#include "utils.h"
#include <algorithm>
#include <basix/mdspan.hpp>
#include <cstdint>
#include <dolfinx/common/types.h>
#include <dolfinx/mesh/EntityMap.h>
#include <memory>
#include <optional>
#include <span>
#include <stdexcept>
#include <vector>

 
namespace dolfinx::fem
{
template <dolfinx::scalar T, std::floating_point U>
class DirichletBC;
template <dolfinx::scalar T, std::floating_point U>
class Expression;
template <dolfinx::scalar T, std::floating_point U>
class Form;
template <std::floating_point T>
class FunctionSpace;

template <dolfinx::scalar T, std::floating_point U>
void tabulate_expression(
    std::span<T> values, const fem::Expression<T, U>& e,
    md::mdspan<const T, md::dextents<std::size_t, 2>> coeffs,
    std::span<const T> constants, const mesh::Mesh<U>& mesh,
    fem::MDSpan2 auto entities,
    std::optional<
        std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
        element)
{
  // Check that domain is the same as mesh of the expression
  if (e.coordinate_element_hash() != mesh.geometry().cmaps().front().hash())
  {
    throw std::invalid_argument(
        "Expression was created on a different mesh. Cannot tabulate.");
  }
  auto [X, Xshape] = e.X();
  impl::tabulate_expression(values, e.kernel(), Xshape, e.value_size(), coeffs,
                            constants, mesh, entities, element);
}

template <dolfinx::scalar T, std::floating_point U>
void tabulate_expression(std::span<T> values, const fem::Expression<T, U>& e,
                         const mesh::Mesh<U>& mesh, fem::MDSpan2 auto entities)
{
  // Check that domain is the same as mesh of the expression
  if (e.coordinate_element_hash() != mesh.geometry().cmaps().front().hash())
  {
    throw std::invalid_argument(
        "Expression was created on a different mesh. Cannot tabulate.");
  }
 
  std::optional<
      std::pair<std::reference_wrapper<const FiniteElement<U>>, std::size_t>>
      element = std::nullopt;
  if (auto V = e.argument_space(); V)
  {
    std::size_t num_argument_dofs
        = V->dofmap()->element_dof_layout().num_dofs() * V->dofmap()->bs();
    assert(V->element());
    element = {std::cref(*V->element()), num_argument_dofs};
  }
 
  std::vector<int> coffsets = e.coefficient_offsets();
  const std::vector<std::shared_ptr<const Function<T, U>>>& coefficients
      = e.coefficients();
  std::vector<T> coeffs(entities.extent(0) * coffsets.back());
  int cstride = coffsets.back();
  {
    std::vector<std::reference_wrapper<const Function<T, U>>> c;
    std::ranges::transform(coefficients, std::back_inserter(c),
                           [](auto c) -> const Function<T, U>& { return *c; });
    fem::pack_coefficients(c, mesh, entities, e.entity_maps(), coffsets,
                           std::span(coeffs));
  }
  std::vector<T> constants = fem::pack_constants(e);
 
  tabulate_expression<T, U>(
      values, e, md::mdspan(coeffs.data(), entities.extent(0), cstride),
      std::span<const T>(constants), mesh, entities, element);
}
 
// -- Helper functions -----------------------------------------------------

template <dolfinx::scalar T>
std::map<std::pair<IntegralType, int>, std::pair<std::span<const T>, int>>
make_coefficients_span(const std::map<std::pair<IntegralType, int>,
                                      std::pair<std::vector<T>, int>>& coeffs)
{
  using Key = typename std::remove_reference_t<decltype(coeffs)>::key_type;
  std::map<Key, std::pair<std::span<const T>, int>> c;
  std::ranges::transform(
      coeffs, std::inserter(c, c.end()),
      [](auto& e) -> typename decltype(c)::value_type
      { return {e.first, {e.second.first, e.second.second}}; });
  return c;
}
 
// -- Scalar ----------------------------------------------------------------

template <dolfinx::scalar T, std::floating_point U>
T assemble_scalar(
    const Form<T, U>& M, std::span<const T> constants,
    const std::map<std::pair<IntegralType, int>,
                   std::pair<std::span<const T>, int>>& coefficients)
{
  using mdspanx3_t
      = md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>>;
 
  std::shared_ptr<const mesh::Mesh<U>> mesh = M.mesh();
  assert(mesh);
  std::span x = mesh->geometry().x();
 
  // Accumulate contributions from each cell type
  const int num_cell_types = mesh->topology()->cell_types().size();
  T val = 0;
  for (int cell_type_idx = 0; cell_type_idx < num_cell_types; ++cell_type_idx)
  {
    // Geometry dofmap and data
    md::mdspan<const std::int32_t, md::dextents<std::size_t, 2>> x_dofmap
        = mesh->geometry().dofmaps().at(cell_type_idx);
    val += impl::assemble_scalar(M, x_dofmap,
                                 mdspanx3_t(x.data(), x.size() / 3, 3),
                                 constants, coefficients, cell_type_idx);
  }
  return val;
}

template <dolfinx::scalar T, std::floating_point U>
T assemble_scalar(const Form<T, U>& M)
{
  const std::vector<T> constants = pack_constants(M);
  auto coefficients = allocate_coefficient_storage(M);
  pack_coefficients(M, coefficients);
  return assemble_scalar(M, std::span(constants),
                         make_coefficients_span(coefficients));
}
 
// -- Vectors ----------------------------------------------------------------

// template <dolfinx::scalar T, std::floating_point U>
template <typename V, std::floating_point U,
          dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
  requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
void assemble_vector(
    V&& b, const Form<T, U>& L, std::span<const T> constants,
    const std::map<std::pair<IntegralType, int>,
                   std::pair<std::span<const T>, int>>& coefficients)
{
  impl::assemble_vector(b, L, constants, coefficients);
}

// template <dolfinx::scalar T, std::floating_point U>
// void assemble_vector(std::span<T> b, const Form<T, U>& L)
template <typename V, std::floating_point U,
          dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
  requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
void assemble_vector(V&& b, const Form<T, U>& L)
{
  auto coefficients = allocate_coefficient_storage(L);
  pack_coefficients(L, coefficients);
  const std::vector<T> constants = pack_constants(L);
  assemble_vector(b, L, std::span(constants),
                  make_coefficients_span(coefficients));
}

template <typename V,
          std::floating_point U
          = scalar_value_t<typename std::remove_cvref_t<V>::value_type>,
          dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
  requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
void apply_lifting(
    V&& b,
    const std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>>&
        a,
    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::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
    const std::vector<std::span<const T>>& x0, T alpha)
{
  // If all forms are null, there is nothing to do
  if (std::ranges::all_of(a, [](auto ai) { return !ai; }))
    return;
 
  common::Timer t("[Apply lifting]");
 
  if (!x0.empty() and x0.size() != a.size())
  {
    throw std::invalid_argument(
        "Mismatch in size between x0 and bilinear form in assembler.");
  }
 
  if (a.size() != bcs1.size())
  {
    throw std::invalid_argument(
        "Mismatch in size between a and bcs in assembler.");
  }
 
  // Reused across iterations so `assign` below can recycle the
  // existing buffer instead of reallocating for every block.
  std::vector<std::int8_t> bc_markers1;
  std::vector<T> bc_values1;
  for (std::size_t j = 0; j < a.size(); ++j)
  {
    if (a[j] and !bcs1[j].empty())
    {
      assert(a[j]->get().function_spaces().at(0));
      auto V1 = a[j]->get().function_spaces()[1];
      assert(V1);
 
      const int bs0 = a[j]->get().function_spaces()[0]->dofmaps().front()->bs();
      const int bs1 = V1->dofmaps().front()->bs();
 
      std::span<const T> _x0;
      if (!x0.empty())
        _x0 = x0[j];
 
      std::shared_ptr<const DofMap> dofmap = V1->dofmaps().front();
      auto map1 = dofmap->index_map;
      const int map_bs1 = dofmap->index_map_bs();
      assert(map1);
      const int crange = map_bs1 * (map1->size_local() + map1->num_ghosts());
      bc_markers1.assign(crange, false);
      bc_values1.assign(crange, 0);
      for (auto& bc : bcs1[j])
      {
        bc.get().mark_dofs(bc_markers1);
        bc.get().set(bc_values1, std::nullopt, 1);
      }
 
      if (bs0 == 1 and bs1 == 1)
      {
        impl::lift_bc(b, a[j]->get(), std::integral_constant<int, 1>{},
                      std::integral_constant<int, 1>{}, constants[j], coeffs[j],
                      std::span<const T>(bc_values1), bc_markers1, _x0, alpha);
      }
      else if (bs0 == 3 and bs1 == 3)
      {
        impl::lift_bc(b, a[j]->get(), std::integral_constant<int, 3>{},
                      std::integral_constant<int, 3>{}, constants[j], coeffs[j],
                      std::span<const T>(bc_values1), bc_markers1, _x0, alpha);
      }
      else
      {
        impl::lift_bc(b, a[j]->get(), bs0, bs1, constants[j], coeffs[j],
                      std::span<const T>(bc_values1), bc_markers1, _x0, alpha);
      }
    }
  }
}

template <typename V,
          std::floating_point U
          = scalar_value_t<typename std::remove_cvref_t<V>::value_type>,
          dolfinx::scalar T = typename std::remove_cvref_t<V>::value_type>
  requires std::is_same_v<typename std::remove_cvref_t<V>::value_type, T>
void apply_lifting(
    V&& b,
    const std::vector<std::optional<std::reference_wrapper<const Form<T, U>>>>&
        a,
    const std::vector<
        std::vector<std::reference_wrapper<const DirichletBC<T, U>>>>& bcs1,
    const std::vector<std::span<const T>>& x0, T alpha)
{
  std::vector<
      std::map<std::pair<IntegralType, int>, std::pair<std::vector<T>, int>>>
      coeffs;
  std::vector<std::vector<T>> constants;
  for (const auto& _a : a)
  {
    if (_a)
    {
      auto coefficients = allocate_coefficient_storage(_a->get());
      pack_coefficients(_a->get(), coefficients);
      coeffs.push_back(coefficients);
      constants.push_back(pack_constants(_a->get()));
    }
    else
    {
      coeffs.emplace_back();
      constants.emplace_back();
    }
  }
 
  std::vector<std::span<const T>> _constants(constants.begin(),
                                             constants.end());
  std::vector<std::map<std::pair<IntegralType, int>,
                       std::pair<std::span<const T>, int>>>
      _coeffs;
  std::ranges::transform(coeffs, std::back_inserter(_coeffs),
                         [](auto& c) { return make_coefficients_span(c); });
 
  apply_lifting(b, a, _constants, _coeffs, bcs1, x0, alpha);
}
 
// -- Matrices ---------------------------------------------------------------

template <dolfinx::scalar T, std::floating_point U>
void assemble_matrix(
    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,
    std::span<const std::int8_t> dof_marker1)
 
{
  common::Timer t_assm("[Assemble Matrix]");
  using mdspanx3_t
      = md::mdspan<const U, md::extents<std::size_t, md::dynamic_extent, 3>>;
 
  std::shared_ptr<const mesh::Mesh<U>> mesh = a.mesh();
  assert(mesh);
  std::span x = mesh->geometry().x();
  impl::assemble_matrix<false>(mat_add, a,
                               mdspanx3_t(x.data(), x.size() / 3, 3), constants,
                               coefficients, dof_marker0, dof_marker1);
}

template <dolfinx::scalar T, std::floating_point U>
void assemble_matrix(
    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,
    const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
{
  // Index maps for dof ranges
  // NOTE: For mixed-topology meshes, there will be multiple DOF maps,
  // but the index maps are the same.
  auto map0 = a.function_spaces().at(0)->dofmaps().front()->index_map;
  auto map1 = a.function_spaces().at(1)->dofmaps().front()->index_map;
  auto bs0 = a.function_spaces().at(0)->dofmaps().front()->index_map_bs();
  auto bs1 = a.function_spaces().at(1)->dofmaps().front()->index_map_bs();
 
  // Build dof markers
  std::vector<std::int8_t> dof_marker0, dof_marker1;
  assert(map0);
  std::int32_t dim0 = bs0 * (map0->size_local() + map0->num_ghosts());
  assert(map1);
  std::int32_t dim1 = bs1 * (map1->size_local() + map1->num_ghosts());
  for (std::size_t k = 0; k < bcs.size(); ++k)
  {
    assert(bcs[k].get().function_space());
    if (a.function_spaces().at(0)->contains(*bcs[k].get().function_space()))
    {
      dof_marker0.resize(dim0, false);
      bcs[k].get().mark_dofs(dof_marker0);
    }
 
    if (a.function_spaces().at(1)->contains(*bcs[k].get().function_space()))
    {
      dof_marker1.resize(dim1, false);
      bcs[k].get().mark_dofs(dof_marker1);
    }
  }
 
  // Assemble
  fem::assemble_matrix(mat_add, a, constants, coefficients, dof_marker0,
                       dof_marker1);
}

template <dolfinx::scalar T, std::floating_point U>
void assemble_matrix(
    auto mat_add, const Form<T, U>& a,
    const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs)
{
  // Prepare constants and coefficients
  const std::vector<T> constants = pack_constants(a);
  auto coefficients = allocate_coefficient_storage(a);
  pack_coefficients(a, coefficients);
 
  // Assemble
  assemble_matrix(mat_add, a, std::span(constants),
                  make_coefficients_span(coefficients), bcs);
}

template <dolfinx::scalar T, std::floating_point U>
void assemble_matrix(auto mat_add, const Form<T, U>& a,
                     std::span<const std::int8_t> dof_marker0,
                     std::span<const std::int8_t> dof_marker1)
 
{
  // Prepare constants and coefficients
  const std::vector<T> constants = pack_constants(a);
  auto coefficients = allocate_coefficient_storage(a);
  pack_coefficients(a, coefficients);
 
  // Assemble
  impl::assemble_matrix<false>(mat_add, a, std::span(constants),
                               make_coefficients_span(coefficients),
                               dof_marker0, dof_marker1);
}

template <dolfinx::scalar T>
void set_diagonal(auto set_fn, std::span<const std::int32_t> rows,
                  T diagonal = 1.0)
{
  for (std::size_t i = 0; i < rows.size(); ++i)
  {
    std::span diag_span(&diagonal, 1);
    set_fn(rows.subspan(i, 1), rows.subspan(i, 1), diag_span);
  }
}

template <dolfinx::scalar T, std::floating_point U>
void set_diagonal(
    auto set_fn, const FunctionSpace<U>& V,
    const std::vector<std::reference_wrapper<const DirichletBC<T, U>>>& bcs,
    T diagonal = 1.0)
{
  spdlog::debug("Set diagonal");
  for (auto& bc : bcs)
  {
    if (V.contains(*bc.get().function_space()))
    {
      const auto [dofs, range] = bc.get().dof_indices();
      set_diagonal(set_fn, dofs.first(range), diagonal);
    }
  }
}
 
} // namespace dolfinx::fem