Geometry.h

// Copyright (C) 2006-2020 Anders Logg and Garth N. Wells
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <dolfinx/common/MPI.h>
#include <dolfinx/fem/CoordinateElement.h>
#include <dolfinx/graph/AdjacencyList.h>
#include <dolfinx/graph/scotch.h>
#include <functional>
#include <memory>
#include <vector>
#include <xtensor/xbuilder.hpp>
#include <xtensor/xtensor.hpp>
#include <xtensor/xview.hpp>
 
namespace dolfinx::common
{
class IndexMap;
}
 
namespace dolfinx::fem
{
class CoordinateElement;
}
 
namespace dolfinx::mesh
{
class Topology;
 
 
class Geometry
{
public:
  template <typename AdjacencyList32, typename Array, typename Vector64>
  Geometry(const std::shared_ptr<const common::IndexMap>& index_map,
           AdjacencyList32&& dofmap, const fem::CoordinateElement& element,
           Array&& x, Vector64&& input_global_indices)
      : _dim(x.shape(1)), _dofmap(std::forward<AdjacencyList32>(dofmap)),
        _index_map(index_map), _cmap(element), _x(std::forward<Array>(x)),
        _input_global_indices(std::forward<Vector64>(input_global_indices))
  {
    assert(_x.shape(1) > 0 and _x.shape(1) <= 3);
    if (_x.shape(0) != _input_global_indices.size())
      throw std::runtime_error("Size mis-match");
 
    // Make all geometry 3D
    if (_dim != 3)
    {
      xt::xtensor<double, 2> c
          = xt::zeros<double>({_x.shape(0), static_cast<std::size_t>(3)});
 
      // The below should work, but misbehaves with the Intel icpx compiler
      // xt::view(c, xt::all(), xt::range(0, _dim)) = _x;
      auto x_view = xt::view(c, xt::all(), xt::range(0, _dim));
      x_view.assign(_x);
 
      std::swap(c, _x);
    }
  }
 
  Geometry(const Geometry&) = default;
 
  Geometry(Geometry&&) = default;
 
  ~Geometry() = default;
 
  Geometry& operator=(const Geometry&) = delete;
 
  Geometry& operator=(Geometry&&) = default;
 
  int dim() const;
 
  const graph::AdjacencyList<std::int32_t>& dofmap() const;
 
  std::shared_ptr<const common::IndexMap> index_map() const;
 
  xt::xtensor<double, 2>& x();
 
  const xt::xtensor<double, 2>& x() const;
 
  const fem::CoordinateElement& cmap() const;
 
  const std::vector<std::int64_t>& input_global_indices() const;
 
private:
  // Geometric dimension
  int _dim;
 
  // Map per cell for extracting coordinate data
  graph::AdjacencyList<std::int32_t> _dofmap;
 
  // IndexMap for geometry 'dofmap'
  std::shared_ptr<const common::IndexMap> _index_map;
 
  // The coordinate element
  fem::CoordinateElement _cmap;
 
  // Coordinates for all points stored as a contiguous array
  xt::xtensor<double, 2> _x;
 
  // Global indices as provided on Geometry creation
  std::vector<std::int64_t> _input_global_indices;
};
 
mesh::Geometry
create_geometry(MPI_Comm comm, const Topology& topology,
                const fem::CoordinateElement& coordinate_element,
                const graph::AdjacencyList<std::int64_t>& cells,
                const xt::xtensor<double, 2>& x,
                const std::function<std::vector<int>(
                    const graph::AdjacencyList<std::int32_t>&)>& reorder_fn
                = nullptr);
 
} // namespace dolfinx::mesh