Remove xtensor from VTK helper interface (#2272)

* Remove some xtensor

* Updates

* Wrapper update

* Tidy up

* Remove comments

* Simplification

* Simplify

* Improve function name

* Add check

* Simplify

* Doc fix

* Smalle simplification

cpp/dolfinx/io/ADIOS2Writers.cpp

@@ -21,7 +21,6 @@
 #include <dolfinx/mesh/utils.h>
 #include <pugixml.hpp>
 #include <xtensor/xtensor.hpp>
-#include <xtensor/xview.hpp>
 
 using namespace dolfinx;
 using namespace dolfinx::io;
@@ -183,10 +182,13 @@ void vtx_write_mesh(adios2::IO& io, adios2::Engine& engine,
   // Pack mesh 'nodes'. Output is written as [N0, v0_0,...., v0_N0, N1,
   // v1_0,...., v1_N1,....], where N is the number of cell nodes and v0,
   // etc, is the node index
-  xt::xtensor<std::int64_t, 2> cells({num_cells, num_nodes + 1});
-  xt::view(cells, xt::all(), xt::xrange(std::size_t(1), cells.shape(1)))
-      = io::extract_vtk_connectivity(mesh);
-  xt::view(cells, xt::all(), 0) = num_nodes;
+  const auto [vtkcells, shape] = io::extract_vtk_connectivity(mesh);
+  std::vector<std::int64_t> cells(shape[0] * (shape[1] + 1), shape[1]);
+  for (std::size_t c = 0; c < shape[0]; ++c)
+  {
+    std::copy_n(std::next(vtkcells.begin(), c * shape[1]), shape[1],
+                std::next(cells.begin(), c * (shape[1] + 1)));
+  }
 
   // Put topology (nodes)
   adios2::Variable<std::int64_t> local_topology = define_variable<std::int64_t>(
@@ -223,31 +225,34 @@ void vtx_write_mesh_from_space(adios2::IO& io, adios2::Engine& engine,
   const int tdim = mesh->topology().dim();
 
   // Get a VTK mesh with points at the 'nodes'
-  const auto [x, x_id, x_ghost, vtk] = io::vtk_mesh_from_space(V);
+  auto [x, xshape, x_id, x_ghost, vtk, vtkshape] = io::vtk_mesh_from_space(V);
 
-  std::uint32_t num_dofs = x.shape(0);
+  std::uint32_t num_dofs = xshape[0];
   std::uint32_t num_cells = mesh->topology().index_map(tdim)->size_local();
 
   // -- Pack mesh 'nodes'. Output is written as [N0, v0_0,...., v0_N0, N1,
   // v1_0,...., v1_N1,....], where N is the number of cell nodes and v0,
   // etc, is the node index.
-  std::vector<std::size_t> shape = {num_cells, vtk.shape(1) + 1};
+  std::vector<std::size_t> shape = {num_cells, vtkshape[1] + 1};
 
   // Create vector, setting all entries to nodes per cell (vtk.shape(1))
-  std::vector<std::uint64_t> vtk_topology(shape[0] * shape[1], vtk.shape(1));
+  std::vector<std::uint64_t> vtk_topology(shape[0] * shape[1], vtkshape[1]);
 
   // Set the [v0_0,...., v0_N0, v1_0,...., v1_N1,....] data
-  auto _vtk = xt::adapt(vtk_topology, shape);
-  xt::view(_vtk, xt::all(), xt::range(1, _vtk.shape(1)))
-      = xt::view(vtk, xt::range(0, num_cells, xt::all()));
+  std::vector<std::int64_t> cells(shape[0] * (shape[1] + 1), shape[1]);
+  for (std::size_t c = 0; c < shape[0]; ++c)
+  {
+    std::copy_n(std::next(vtk.begin(), c * vtkshape[1]), vtkshape[1],
+                std::next(vtk_topology.begin(), c * shape[1]));
+  }
 
   // Define ADIOS2 variables for geometry, topology, celltypes and
   // corresponding VTK data
   adios2::Variable<double> local_geometry
       = define_variable<double>(io, "geometry", {}, {}, {num_dofs, 3});
   adios2::Variable<std::uint64_t> local_topology
       = define_variable<std::uint64_t>(io, "connectivity", {}, {},
-                                       {num_cells, vtk.shape(1) + 1});
+                                       {num_cells, vtkshape[1] + 1});
   adios2::Variable<std::uint32_t> cell_type
       = define_variable<std::uint32_t>(io, "types");
   adios2::Variable<std::uint32_t> vertices = define_variable<std::uint32_t>(
@@ -548,7 +553,7 @@ void fides_write_mesh(adios2::IO& io, adios2::Engine& engine,
   const int tdim = topology.dim();
   const std::int32_t num_cells = topology.index_map(tdim)->size_local();
   const int num_nodes = geometry.cmap().dim();
-  const xt::xtensor<std::int64_t, 2> cells = extract_vtk_connectivity(mesh);
+  const auto [cells, shape] = io::extract_vtk_connectivity(mesh);
 
   // "Put" topology data in the result in the ADIOS2 file
   adios2::Variable<std::int64_t> local_topology = define_variable<std::int64_t>(

cpp/dolfinx/io/VTKFile.cpp

@@ -23,9 +23,6 @@
 #include <span>
 #include <sstream>
 #include <string>
-#include <xtensor/xbuilder.hpp>
-#include <xtensor/xcomplex.hpp>
-#include <xtensor/xview.hpp>
 
 using namespace dolfinx;
 
@@ -185,7 +182,8 @@ void add_data(const std::string& name, int rank,
 
 /// Add mesh geometry and topology data to a pugixml node. This function
 /// adds the Points and Cells nodes to the input node.
-/// @param[in] x Coordinates of the points
+/// @param[in] x Coordinates of the points, row-major storage
+/// @param[in] xshape The shape of `x`
 /// @param[in] x_id Unique global index for each point
 /// @param[in] x_ghost Flag indicating if a point is a owned (0) or is a
 /// ghost (1)
@@ -194,10 +192,12 @@ void add_data(const std::string& name, int rank,
 /// @param[in] celltype The cell type
 /// @param[in] tdim Topological dimension of the cells
 /// @param[in,out] piece_node The XML node to add data to
-void add_mesh(const xt::xtensor<double, 2>& x,
+void add_mesh(const std::span<const double>& x,
+              std::array<std::size_t, 2> /*xshape*/,
               const std::span<const std::int64_t> x_id,
               const std::span<const std::uint8_t> x_ghost,
-              const xt::xtensor<std::int64_t, 2>& cells,
+              const std::span<const std::int64_t>& cells,
+              std::array<std::size_t, 2> cshape,
               const common::IndexMap& cellmap, mesh::CellType celltype,
               int tdim, pugi::xml_node& piece_node)
 {
@@ -232,8 +232,8 @@ void add_mesh(const xt::xtensor<double, 2>& x,
   offsets_node.append_attribute("format") = "ascii";
   {
     std::stringstream ss;
-    int num_nodes = cells.shape(1);
-    for (std::size_t i = 0; i < cells.shape(0); ++i)
+    int num_nodes = cshape[1];
+    for (std::size_t i = 0; i < cshape[0]; ++i)
       ss << (i + 1) * num_nodes << " ";
     offsets_node.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
   }
@@ -245,7 +245,7 @@ void add_mesh(const xt::xtensor<double, 2>& x,
   int vtk_celltype = io::cells::get_vtk_cell_type(celltype, tdim);
   {
     std::stringstream ss;
-    for (std::size_t c = 0; c < cells.shape(0); ++c)
+    for (std::size_t c = 0; c < cshape[0]; ++c)
       ss << vtk_celltype << " ";
     type_node.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
   }
@@ -263,7 +263,7 @@ void add_mesh(const xt::xtensor<double, 2>& x,
     std::stringstream ss;
     for (std::int32_t c = 0; c < cellmap.size_local(); ++c)
       ss << 0 << " ";
-    for (std::size_t c = cellmap.size_local(); c < cells.shape(0); ++c)
+    for (std::size_t c = cellmap.size_local(); c < cshape[0]; ++c)
       ss << 1 << " ";
     ghost_cell_node.append_child(pugi::node_pcdata).set_value(ss.str().c_str());
   }
@@ -411,34 +411,40 @@ void write_function(
   pugi::xml_node grid_node_vtu = vtk_node_vtu.append_child("UnstructuredGrid");
 
   // Build mesh data using first FunctionSpace
-  xt::xtensor<double, 2> x;
+  std::vector<double> x;
+  std::array<std::size_t, 2> xshape;
   std::vector<std::int64_t> x_id;
   std::vector<std::uint8_t> x_ghost;
-  xt::xtensor<std::int32_t, 2> cells;
+  std::vector<std::int64_t> cells;
+  std::array<std::size_t, 2> cshape;
   if (is_cellwise(*V0))
   {
-    cells = io::extract_vtk_connectivity(*mesh0);
+    std::vector<std::int64_t> tmp;
+    std::tie(tmp, cshape) = io::extract_vtk_connectivity(*mesh0);
+    cells.assign(tmp.begin(), tmp.end());
     const mesh::Geometry& geometry = mesh0->geometry();
-    x = xt::adapt(geometry.x().data(), geometry.x().size(), xt::no_ownership(),
-                  std::vector({geometry.x().size() / 3, std::size_t(3)}));
+    x.assign(geometry.x().begin(), geometry.x().end());
+    xshape = {geometry.x().size() / 3, 3};
     x_id = geometry.input_global_indices();
     auto xmap = geometry.index_map();
     assert(xmap);
-    x_ghost.resize(x.shape(0), 0);
+    x_ghost.resize(xshape[0], 0);
     std::fill(std::next(x_ghost.begin(), xmap->size_local()), x_ghost.end(), 1);
   }
   else
-    std::tie(x, x_id, x_ghost, cells) = io::vtk_mesh_from_space(*V0);
+    std::tie(x, xshape, x_id, x_ghost, cells, cshape)
+        = io::vtk_mesh_from_space(*V0);
 
   // Add "Piece" node and required metadata
   pugi::xml_node piece_node = grid_node_vtu.append_child("Piece");
-  piece_node.append_attribute("NumberOfPoints") = x.shape(0);
-  piece_node.append_attribute("NumberOfCells") = cells.shape(0);
+  piece_node.append_attribute("NumberOfPoints") = xshape[0];
+  piece_node.append_attribute("NumberOfCells") = cshape[0];
 
   // Add mesh data to "Piece" node
   int tdim = mesh0->topology().dim();
-  add_mesh(x, x_id, x_ghost, cells, *mesh0->topology().index_map(tdim),
-           mesh0->topology().cell_type(), mesh0->topology().dim(), piece_node);
+  add_mesh(x, xshape, x_id, x_ghost, cells, cshape,
+           *mesh0->topology().index_map(tdim), mesh0->topology().cell_type(),
+           mesh0->topology().dim(), piece_node);
 
   // FIXME: is this actually setting the first?
   // Set last scalar/vector/tensor Functions in u to be the 'active'
@@ -475,9 +481,9 @@ void write_function(
       assert(!data_node.empty());
       auto dofmap = V->dofmap();
       int bs = dofmap->bs();
-      std::vector<T> data(cells.shape(0) * num_comp, 0);
+      std::vector<T> data(cshape[0] * num_comp, 0);
       auto u_vector = _u.get().x()->array();
-      for (std::size_t c = 0; c < cells.shape(0); ++c)
+      for (std::size_t c = 0; c < cshape[0]; ++c)
       {
         std::span<const std::int32_t> dofs = dofmap->cell_dofs(c);
         for (std::size_t i = 0; i < dofs.size(); ++i)
@@ -541,7 +547,7 @@ void write_function(
         // Interpolate on each cell
         auto u_vector = _u.get().x()->array();
         std::vector<T> u(u_vector.size());
-        for (std::size_t c = 0; c < cells.shape(0); ++c)
+        for (std::size_t c = 0; c < cshape[0]; ++c)
         {
           std::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(c);
           std::span<const std::int32_t> dofs = dofmap->cell_dofs(c);
@@ -765,15 +771,13 @@ void io::VTKFile::write(const mesh::Mesh& mesh, double time)
   piece_node.append_attribute("NumberOfCells") = num_cells;
 
   // Add mesh data to "Piece" node
-  xt::xtensor<std::int64_t, 2> cells = extract_vtk_connectivity(mesh);
-  xt::xtensor<double, 2> x
-      = xt::adapt(geometry.x().data(), geometry.x().size(), xt::no_ownership(),
-                  std::vector({geometry.x().size() / 3, std::size_t(3)}));
-  std::vector<std::uint8_t> x_ghost(x.shape(0), 0);
+  const auto [cells, cshape] = extract_vtk_connectivity(mesh);
+  std::array<std::size_t, 2> xshape = {geometry.x().size() / 3, 3};
+  std::vector<std::uint8_t> x_ghost(xshape[0], 0);
   std::fill(std::next(x_ghost.begin(), xmap->size_local()), x_ghost.end(), 1);
-  add_mesh(x, geometry.input_global_indices(), x_ghost, cells,
-           *topology.index_map(tdim), topology.cell_type(), topology.dim(),
-           piece_node);
+  add_mesh(geometry.x(), xshape, geometry.input_global_indices(), x_ghost,
+           cells, cshape, *topology.index_map(tdim), topology.cell_type(),
+           topology.dim(), piece_node);
 
   // Create filepath for a .vtu file
   auto create_vtu_path = [file_root = _filename.parent_path(),

cpp/dolfinx/io/XDMFFile.cpp

@@ -346,8 +346,10 @@ XDMFFile::read_meshtags(const std::shared_ptr<const mesh::Mesh>& mesh,
   mesh::CellType cell_type = mesh::to_type(cell_type_str.first);
 
   // Permute entities from VTK to DOLFINx ordering
-  xt::xtensor<std::int64_t, 2> entities1 = io::cells::compute_permutation(
-      entities, io::cells::perm_vtk(cell_type, entities.shape(1)));
+  std::vector<std::int64_t> entities1 = io::cells::apply_permutation(
+      std::span(entities.data(), entities.size()),
+      {entities.shape(0), entities.shape(1)},
+      io::cells::perm_vtk(cell_type, entities.shape(1)));
 
   std::pair<std::vector<std::int32_t>, std::vector<std::int32_t>>
       entities_values = xdmf_utils::distribute_entity_data(

cpp/dolfinx/io/cells.cpp

@@ -10,7 +10,6 @@
 #include <dolfinx/mesh/cell_types.h>
 #include <numeric>
 #include <stdexcept>
-#include <xtensor/xview.hpp>
 
 using namespace dolfinx;
 namespace
@@ -364,17 +363,21 @@ io::cells::transpose(const std::vector<std::uint8_t>& map)
   return transpose;
 }
 //-----------------------------------------------------------------------------
-xt::xtensor<std::int64_t, 2>
-io::cells::compute_permutation(const xt::xtensor<std::int64_t, 2>& cells,
-                               const std::vector<std::uint8_t>& p)
+std::vector<std::int64_t>
+io::cells::apply_permutation(const std::span<const std::int64_t>& cells,
+                             std::array<std::size_t, 2> shape,
+                             const std::span<const std::uint8_t>& p)
 {
+  assert(cells.size() == shape[0] * shape[1]);
+  assert(shape[1] == p.size());
+
   LOG(INFO) << "IO permuting cells";
-  xt::xtensor<std::int64_t, 2> cells_new(cells.shape());
-  for (std::size_t c = 0; c < cells_new.shape(0); ++c)
+  std::vector<std::int64_t> cells_new(cells.size());
+  for (std::size_t c = 0; c < shape[0]; ++c)
   {
-    auto cell = xt::row(cells, c);
-    auto cell_new = xt::row(cells_new, c);
-    for (std::size_t i = 0; i < cell_new.shape(0); ++i)
+    auto cell = cells.subspan(c * shape[1], shape[1]);
+    std::span cell_new(cells_new.data() + c * shape[1], shape[1]);
+    for (std::size_t i = 0; i < shape[1]; ++i)
       cell_new[i] = cell[p[i]];
   }
   return cells_new;

cpp/dolfinx/io/cells.h

@@ -6,10 +6,11 @@
 
 #pragma once
 
+#include <array>
 #include <cstdint>
 #include <dolfinx/mesh/cell_types.h>
+#include <span>
 #include <vector>
-#include <xtensor/xtensor.hpp>
 
 /// Functions for the re-ordering of input mesh topology to the DOLFINx
 /// ordering, and transpose orderings for file output.
@@ -106,14 +107,17 @@ std::vector<std::uint8_t> transpose(const std::vector<std::uint8_t>& map);
 
 /// Permute cell topology by applying a permutation array for each cell
 /// @param[in] cells Array of cell topologies, with each row
-/// representing a cell
+/// representing a cell (row-major storage)
+/// @param[in] shape The shape of the `cells` array
 /// @param[in] p The permutation array that maps `a_p[i] = a[p[i]]`,
 /// where `a_p` is the permuted array
 /// @return Permuted cell topology, where for a cell `v_new[i] =
-/// v_old[map[i]]`
-xt::xtensor<std::int64_t, 2>
-compute_permutation(const xt::xtensor<std::int64_t, 2>& cells,
-                    const std::vector<std::uint8_t>& p);
+/// v_old[map[i]]`. The storage is row-major and the shape is the same
+/// as `cells`.
+std::vector<std::int64_t>
+apply_permutation(const std::span<const std::int64_t>& cells,
+                  std::array<std::size_t, 2> shape,
+                  const std::span<const std::uint8_t>& p);
 
 /// Get VTK cell identifier
 /// @param[in] cell The cell type