Function test_matrix_apply

cpp/test/matrix.cpp:66–120

}

void test_matrix_apply()
{
  MPI_Comm comm = MPI_COMM_WORLD;
  auto part = mesh::create_cell_partitioner(mesh::GhostMode::none, 2);
  auto mesh = std::make_shared<mesh::Mesh<double>>(
      mesh::create_box(comm, {{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}}, {12, 12, 12},
                       mesh::CellType::tetrahedron, part));

  auto element = basix::create_element<double>(
      basix::element::family::P, basix::cell::type::tetrahedron, 2,
      basix::element::lagrange_variant::unset,
      basix::element::dpc_variant::unset, false);

  auto V = std::make_shared<fem::FunctionSpace<double>>(
      fem::create_functionspace<double>(
          mesh, std::make_shared<fem::FiniteElement<double>>(element)));

  // Prepare and set Constants for the bilinear form
  auto kappa = std::make_shared<fem::Constant<double>>(2.0);
  auto ui = std::make_shared<fem::Function<double, double>>(V);

  // Define variational forms
  auto a = std::make_shared<fem::Form<double, double>>(
      fem::create_form<double, double>(*form_poisson_a, {V, V}, {},
                                       {{"kappa", kappa}}, {}, {}));

  // Create sparsity pattern
  la::SparsityPattern sp = fem::create_sparsity_pattern(*a);
  sp.finalize();

  // Assemble matrix
  la::MatrixCSR<double> A(sp);
  fem::assemble_matrix(A.mat_add_values(), *a, {});
  A.scatter_rev();
  CHECK((V->dofmap()->index_map->size_local() == A.num_owned_rows()));

  // Get compatible vectors
  auto col_map = A.index_map(1);

  la::Vector<double> x(col_map, 1);
  la::Vector<double> y(col_map, 1);

  std::size_t col_size = col_map->size_local() + col_map->num_ghosts();
  CHECK(x.array().size() == col_size);

  // Fill x vector with 1 (Constant)
  std::ranges::fill(x.array(), 1);

  // Matrix A represents the action of the Laplace operator, so when
  // applied to a constant vector the result should be zero
  A.mult(x, y);

  std::ranges::for_each(y.array(),
                        [](auto a) { REQUIRE(std::abs(a) < 1e-13); });
}

void test_matrix_cast()
{