Select Git revision
Connectivity.hpp
Connectivity.hpp 17.56 KiB
#ifndef CONNECTIVITY_HPP
#define CONNECTIVITY_HPP
#include <PugsAssert.hpp>
#include <PugsMacros.hpp>
#include <PugsUtils.hpp>
#include <PugsTraits.hpp>
#include <TinyVector.hpp>
#include <ItemValue.hpp>
#include <IConnectivity.hpp>
#include <ConnectivityMatrix.hpp>
#include <ItemToItemMatrix.hpp>
#include <ConnectivityComputer.hpp>
#include <SubItemValuePerItem.hpp>
#include <algorithm>
#include <vector>
#include <CellType.hpp>
#include <CSRGraph.hpp>
#include <ItemType.hpp>
#include <RefId.hpp>
#include <RefItemList.hpp>
#include <SynchronizerManager.hpp>
#include <iostream>
#include <set>
class ConnectivityDescriptor;
template <size_t Dim>
class Connectivity final : public IConnectivity
{
public:
PUGS_INLINE
static std::shared_ptr<Connectivity<Dim>> build(const ConnectivityDescriptor&);
private:
constexpr static auto& itemTId = ItemTypeId<Dim>::itemTId;
public:
static constexpr size_t Dimension = Dim;
PUGS_INLINE
size_t
dimension() const final
{
return Dimension;
}
private:
ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1];
WeakCellValue<const CellType> m_cell_type;
WeakCellValue<const int> m_cell_global_index;
WeakCellValue<const int> m_cell_number;
WeakFaceValue<const int> m_face_number;
WeakEdgeValue<const int> m_edge_number;
WeakNodeValue<const int> m_node_number;
WeakCellValue<const int> m_cell_owner;
WeakCellValue<const bool> m_cell_is_owned;
WeakFaceValue<const int> m_face_owner;
WeakFaceValue<const bool> m_face_is_owned;
WeakEdgeValue<const int> m_edge_owner;
WeakEdgeValue<const bool> m_edge_is_owned;
WeakNodeValue<const int> m_node_owner;
WeakNodeValue<const bool> m_node_is_owned;
WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;
WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
WeakFaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces;
WeakCellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells;
WeakEdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges;
WeakNodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes;
WeakCellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells;
WeakFaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces;
WeakNodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes;
WeakCellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells;
WeakEdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges;
WeakFaceValuePerNode<const unsigned short> m_node_local_numbers_in_their_faces;
ConnectivityComputer m_connectivity_computer;
std::vector<RefCellList> m_ref_cell_list_vector;
std::vector<RefFaceList> m_ref_face_list_vector;
std::vector<RefEdgeList> m_ref_edge_list_vector;
std::vector<RefNodeList> m_ref_node_list_vector;
WeakCellValue<const double> m_inv_cell_nb_nodes;
void _computeCellFaceAndFaceNodeConnectivities();
template <typename SubItemValuePerItemType>
PUGS_INLINE const SubItemValuePerItemType&
_lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const
{
using ReversedItemOfItem = typename SubItemValuePerItemType::ItemOfItemType::Reversed;
if (not sub_item_value_per_item.isBuilt()) {
const_cast<SubItemValuePerItemType&>(sub_item_value_per_item) =
m_connectivity_computer.computeLocalItemNumberInChildItem<ReversedItemOfItem>(*this);
}
return sub_item_value_per_item;
}
friend class ConnectivityComputer;
PUGS_INLINE
const ConnectivityMatrix&
_getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const
{
const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
if (not connectivity_matrix.isBuilt()) {
const_cast<ConnectivityMatrix&>(connectivity_matrix) =
m_connectivity_computer.computeConnectivityMatrix(*this, item_type_0, item_type_1);
}
return connectivity_matrix;
}
public: PUGS_INLINE
const auto&
cellType() const
{
return m_cell_type;
}
PUGS_INLINE
const auto&
cellNumber() const
{
return m_cell_number;
}
PUGS_INLINE
const auto&
faceNumber() const
{
return m_face_number;
}
PUGS_INLINE
const auto&
edgeNumber() const
{
return m_edge_number;
}
PUGS_INLINE
const auto&
nodeNumber() const
{
return m_node_number;
}
template <ItemType item_type>
PUGS_INLINE const auto&
number() const
{
if constexpr (item_type == ItemType::cell) {
return m_cell_number;
} else if constexpr (item_type == ItemType::face) {
return m_face_number;
} else if constexpr (item_type == ItemType::edge) {
return m_edge_number;
} else if constexpr (item_type == ItemType::node) {
return m_node_number;
} else {
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
PUGS_INLINE
const auto&
cellOwner() const
{
return m_cell_owner;
}
PUGS_INLINE
const auto&
faceOwner() const
{
return m_face_owner;
}
PUGS_INLINE
const auto&
edgeOwner() const
{ std::cerr << __FILE__ << ':' << __LINE__ << ": edge owner not built\n";
std::terminate();
return m_edge_owner;
}
PUGS_INLINE
const auto&
nodeOwner() const
{
return m_node_owner;
}
template <ItemType item_type>
PUGS_INLINE const auto&
owner() const
{
if constexpr (item_type == ItemType::cell) {
return m_cell_owner;
} else if constexpr (item_type == ItemType::face) {
return m_face_owner;
} else if constexpr (item_type == ItemType::edge) {
return m_edge_owner;
} else if constexpr (item_type == ItemType::node) {
return m_node_owner;
} else {
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
PUGS_INLINE
const auto&
cellIsOwned() const
{
return m_cell_is_owned;
}
PUGS_INLINE
const auto&
faceIsOwned() const
{
return m_face_is_owned;
}
PUGS_INLINE
const auto&
edgeIsOwned() const
{
std::cerr << __FILE__ << ':' << __LINE__ << ": edge is owned not built\n";
std::terminate();
return m_edge_is_owned;
}
PUGS_INLINE
const auto&
nodeIsOwned() const
{
return m_node_is_owned;
}
template <ItemType item_type>
PUGS_INLINE const auto&
isOwned() const
{
if constexpr (item_type == ItemType::cell) {
return m_cell_is_owned;
} else if constexpr (item_type == ItemType::face) {
return m_face_is_owned;
} else if constexpr (item_type == ItemType::edge) {
return m_edge_is_owned;
} else if constexpr (item_type == ItemType::node) { return m_node_is_owned;
} else {
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
PUGS_INLINE
const bool&
isConnectivityMatrixBuilt(const ItemType& item_type_0, const ItemType& item_type_1) const
{
const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
return connectivity_matrix.isBuilt();
}
template <ItemType source_item_type, ItemType target_item_type>
PUGS_INLINE auto
getItemToItemMatrix() const
{
return ItemToItemMatrix<source_item_type, target_item_type>(_getMatrix(source_item_type, target_item_type));
}
PUGS_INLINE
auto
cellToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
}
PUGS_INLINE
auto
cellToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
}
PUGS_INLINE
auto
cellToNodeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
}
PUGS_INLINE
auto
faceToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
}
PUGS_INLINE
auto
faceToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
}
PUGS_INLINE
auto
faceToNodeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
}
PUGS_INLINE
auto
edgeToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
}
PUGS_INLINE auto
edgeToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
}
PUGS_INLINE
auto
edgeToNodeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
}
PUGS_INLINE
auto
nodeToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
}
PUGS_INLINE
auto
nodeToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
}
PUGS_INLINE
auto
nodeToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>();
}
PUGS_INLINE
const auto&
cellFaceIsReversed() const
{
static_assert(Dimension > 1, "reversed faces makes no sense in dimension 1");
return m_cell_face_is_reversed;
}
PUGS_INLINE
const auto&
faceEdgeIsReversed() const
{
static_assert(Dimension > 2, "reversed edges makes no sense in dimension 1 or 2");
return m_face_edge_is_reversed;
}
PUGS_INLINE
const auto&
cellLocalNumbersInTheirNodes() const
{
return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes);
}
PUGS_INLINE
const auto&
cellLocalNumbersInTheirEdges() const
{
if constexpr (Dimension > 2) {
return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges);
} else {
return cellLocalNumbersInTheirFaces();
}
}
PUGS_INLINE
const auto& cellLocalNumbersInTheirFaces() const
{
if constexpr (Dimension > 1) {
return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces);
} else {
return cellLocalNumbersInTheirNodes();
}
}
PUGS_INLINE
const auto&
faceLocalNumbersInTheirCells() const
{
if constexpr (Dimension > 1) {
return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells);
} else {
return nodeLocalNumbersInTheirCells();
}
}
PUGS_INLINE
const auto&
faceLocalNumbersInTheirEdges() const
{
static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_edges);
}
PUGS_INLINE
const auto&
faceLocalNumbersInTheirNodes() const
{
static_assert(Dimension > 1, "this function has no meaning in 1d");
return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes);
}
PUGS_INLINE
const auto&
edgeLocalNumbersInTheirCells() const
{
if constexpr (Dimension > 2) {
return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells);
} else {
return faceLocalNumbersInTheirCells();
}
}
PUGS_INLINE
const auto&
edgeLocalNumbersInTheirFaces() const
{
static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_faces);
}
PUGS_INLINE
const auto&
edgeLocalNumbersInTheirNodes() const
{
static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_nodes);
}
PUGS_INLINE
const auto&
nodeLocalNumbersInTheirCells() const
{
return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells);
}
PUGS_INLINE
const auto&
nodeLocalNumbersInTheirEdges() const
{
static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_edges);
}
PUGS_INLINE
const auto&
nodeLocalNumbersInTheirFaces() const
{
static_assert(Dimension > 1, "this function has no meaning in 1d");
return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
}
template <ItemType item_type>
size_t
numberOfRefItemList() const
{
if constexpr (item_type == ItemType::cell) {
return m_ref_cell_list_vector.size();
} else if constexpr (item_type == ItemType::face) {
return m_ref_face_list_vector.size();
} else if constexpr (item_type == ItemType::edge) {
return m_ref_edge_list_vector.size();
} else if constexpr (item_type == ItemType::node) {
return m_ref_node_list_vector.size();
} else {
static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
}
}
template <ItemType item_type>
const RefItemList<item_type>&
refItemList(const size_t& i) const
{
if constexpr (item_type == ItemType::cell) {
return m_ref_cell_list_vector[i];
} else if constexpr (item_type == ItemType::face) {
return m_ref_face_list_vector[i];
} else if constexpr (item_type == ItemType::edge) {
return m_ref_edge_list_vector[i];
} else if constexpr (item_type == ItemType::node) {
return m_ref_node_list_vector[i];
} else {
static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
}
}
template <ItemType item_type>
void
addRefItemList(const RefItemList<item_type>& ref_item_list)
{
if constexpr (item_type == ItemType::cell) {
m_ref_cell_list_vector.push_back(ref_item_list);
} else if constexpr (item_type == ItemType::face) {
m_ref_face_list_vector.push_back(ref_item_list);
} else if constexpr (item_type == ItemType::edge) {
m_ref_edge_list_vector.push_back(ref_item_list);
} else if constexpr (item_type == ItemType::node) {
m_ref_node_list_vector.push_back(ref_item_list);
} else {
static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
}
}
PUGS_INLINE
CSRGraph
cellToCellGraph() const {
std::vector<std::set<int>> cell_cells(this->numberOfCells());
if constexpr (true) {
const auto& face_to_cell_matrix = this->faceToCellMatrix();
for (FaceId l = 0; l < this->numberOfFaces(); ++l) {
const auto& face_to_cell = face_to_cell_matrix[l];
if (face_to_cell.size() > 1) {
const CellId cell_0 = face_to_cell[0];
const CellId cell_1 = face_to_cell[1];
cell_cells[cell_0].insert(cell_1);
cell_cells[cell_1].insert(cell_0);
}
}
} else {
const auto& node_to_cell_matrix = this->nodeToCellMatrix();
for (NodeId l = 0; l < this->numberOfNodes(); ++l) {
const auto& node_to_cell = node_to_cell_matrix[l];
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_0 = node_to_cell[i_cell];
for (size_t j_cell = 0; j_cell < i_cell; ++j_cell) {
const CellId cell_1 = node_to_cell[j_cell];
cell_cells[cell_0].insert(cell_1);
cell_cells[cell_1].insert(cell_0);
}
}
}
}
Array<int> entries(this->numberOfCells() + 1);
entries[0] = 0;
for (size_t j = 0; j < this->numberOfCells(); ++j) {
entries[j + 1] = entries[j] + cell_cells[j].size();
}
Array<int> neighbors(entries[this->numberOfCells()]);
{
size_t k = 0;
for (size_t j = 0; j < this->numberOfCells(); ++j) {
for (CellId cell_id : cell_cells[j]) {
neighbors[k] = m_cell_global_index[cell_id];
++k;
}
}
}
return CSRGraph(entries, neighbors);
}
PUGS_INLINE
size_t
numberOfNodes() const final
{
const auto& node_to_cell_matrix = this->_getMatrix(ItemType::node, ItemType::cell);
return node_to_cell_matrix.numRows();
}
PUGS_INLINE
size_t
numberOfEdges() const final
{
const auto& edge_to_node_matrix = this->_getMatrix(ItemType::edge, ItemType::node);
return edge_to_node_matrix.numRows();
}
PUGS_INLINE
size_t
numberOfFaces() const final
{
const auto& face_to_node_matrix = this->_getMatrix(ItemType::face, ItemType::cell); return face_to_node_matrix.numRows();
}
PUGS_INLINE
size_t
numberOfCells() const final
{
const auto& cell_to_node_matrix = this->_getMatrix(ItemType::cell, ItemType::node);
return cell_to_node_matrix.numRows();
}
CellValue<const double>
invCellNbNodes() const
{
if (not m_inv_cell_nb_nodes.isBuilt()) {
const auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
WeakCellValue<double> inv_cell_nb_nodes(*this);
parallel_for(
this->numberOfCells(), PUGS_LAMBDA(const CellId& j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
inv_cell_nb_nodes[j] = 1. / cell_nodes.length;
});
const_cast<WeakCellValue<const double>&>(m_inv_cell_nb_nodes) = inv_cell_nb_nodes;
}
return m_inv_cell_nb_nodes;
}
Connectivity(const Connectivity&) = delete;
private:
Connectivity();
void _buildFrom(const ConnectivityDescriptor& descriptor);
public:
~Connectivity()
{
auto& manager = SynchronizerManager::instance();
manager.deleteConnectivitySynchronizer(this);
}
};
template <size_t Dim>
PUGS_INLINE std::shared_ptr<Connectivity<Dim>>
Connectivity<Dim>::build(const ConnectivityDescriptor& descriptor)
{
std::shared_ptr<Connectivity<Dim>> connectivity_ptr(new Connectivity<Dim>);
connectivity_ptr->_buildFrom(descriptor);
return connectivity_ptr;
}
using Connectivity3D = Connectivity<3>;
using Connectivity2D = Connectivity<2>;
using Connectivity1D = Connectivity<1>;
#endif // CONNECTIVITY_HPP