Select Git revision
Connectivity.hpp
Stéphane Del Pino authored
This number is just an identifier and should not be used for anything else Also prepared global index item structure. This has to be a contiguous numbering starting from 0. This will have to take into account ghost cells.
Connectivity.hpp 16.78 KiB
#ifndef CONNECTIVITY_HPP
#define CONNECTIVITY_HPP
#include <PastisMacros.hpp>
#include <PastisAssert.hpp>
#include <PastisOStream.hpp>
#include <PastisUtils.hpp>
#include <TinyVector.hpp>
#include <ItemValue.hpp>
#include <IConnectivity.hpp>
#include <ConnectivityMatrix.hpp>
#include <ItemToItemMatrix.hpp>
#include <SubItemValuePerItem.hpp>
#include <ConnectivityComputer.hpp>
#include <vector>
#include <unordered_map>
#include <algorithm>
#include <CellType.hpp>
#include <CSRGraph.hpp>
#include <RefId.hpp>
#include <ItemType.hpp>
#include <RefNodeList.hpp>
#include <RefFaceList.hpp>
#include <tuple>
#include <algorithm>
#include <set>
template <size_t Dimension>
class Connectivity;
template <size_t Dimension>
class ConnectivityFace;
template<>
class ConnectivityFace<1>
{
public:
friend struct Hash;
struct Hash
{
size_t operator()(const ConnectivityFace& f) const;
};
};
template<>
class ConnectivityFace<2>
{
public:
friend struct Hash;
struct Hash
{
size_t operator()(const ConnectivityFace& f) const {
size_t hash = 0;
hash ^= std::hash<unsigned int>()(f.m_node0_id);
hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
return hash;
}
};
unsigned int m_node0_id;
unsigned int m_node1_id;
friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
{
os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
return os;
}
PASTIS_INLINE
bool operator==(const ConnectivityFace& f) const
{
return ((m_node0_id == f.m_node0_id) and
(m_node1_id == f.m_node1_id));
}
PASTIS_INLINE
bool operator<(const ConnectivityFace& f) const
{
return ((m_node0_id<f.m_node0_id) or
((m_node0_id == f.m_node0_id) and
(m_node1_id<f.m_node1_id)));
}
PASTIS_INLINE
ConnectivityFace& operator=(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace& operator=(ConnectivityFace&&) = default;
PASTIS_INLINE
ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
{
Assert(given_node_id_list.size()==2);
#warning rework this dirty constructor
const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]);
m_node0_id = min;
m_node1_id = max;
}
PASTIS_INLINE
ConnectivityFace(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace(ConnectivityFace&&) = default;
PASTIS_INLINE
~ConnectivityFace() = default;
};
template <>
class ConnectivityFace<3>
{
private:
friend class Connectivity<3>;
friend struct Hash;
struct Hash
{
size_t operator()(const ConnectivityFace& f) const {
size_t hash = 0;
for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
}
return hash;
}
};
bool m_reversed;
std::vector<NodeId::base_type> m_node_id_list;
friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
{
for (auto id : f.m_node_id_list) {
os << id << ' ';
}
return os;
}
PASTIS_INLINE
const bool& reversed() const
{
return m_reversed;
}
PASTIS_INLINE
const std::vector<unsigned int>& nodeIdList() const
{
return m_node_id_list;
}
PASTIS_INLINE
std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
{
const auto min_id = std::min_element(node_list.begin(), node_list.end());
const int shift = std::distance(node_list.begin(), min_id);
std::vector<unsigned int> rotated_node_list(node_list.size());
if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
for (size_t i=0; i<node_list.size(); ++i) {
rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
m_reversed = true;
}
} else {
for (size_t i=0; i<node_list.size(); ++i) {
rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
}
}
return rotated_node_list;
}
PASTIS_INLINE
ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
: m_reversed(false),
m_node_id_list(_sort(given_node_id_list))
{
;
}
public:
bool operator==(const ConnectivityFace& f) const
{
if (m_node_id_list.size() == f.nodeIdList().size()) {
for (size_t j=0; j<m_node_id_list.size(); ++j) {
if (m_node_id_list[j] != f.nodeIdList()[j]) {
return false;
}
}
return true;
}
return false;
}
PASTIS_INLINE
bool operator<(const ConnectivityFace& f) const
{
const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
for (size_t i=0; i<min_nb_nodes; ++i) {
if (m_node_id_list[i] < f.m_node_id_list[i]) return true;
if (m_node_id_list[i] != f.m_node_id_list[i]) return false; }
return m_node_id_list.size() < f.m_node_id_list.size();
}
PASTIS_INLINE
ConnectivityFace& operator=(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace& operator=(ConnectivityFace&&) = default;
PASTIS_INLINE
ConnectivityFace(const ConnectivityFace&) = default;
PASTIS_INLINE
ConnectivityFace(ConnectivityFace&&) = default;
PASTIS_INLINE
ConnectivityFace() = delete;
PASTIS_INLINE
~ConnectivityFace() = default;
};
template <size_t Dimension>
class Connectivity final
: public IConnectivity
{
private:
constexpr static auto& itemTId = ItemTypeId<Dimension>::itemTId;
public:
static constexpr size_t dimension = Dimension;
private:
ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
CellValue<const CellType> m_cell_type;
CellValue<const int> m_cell_global_index;
CellValue<const int> m_cell_number;
FaceValuePerCell<const bool> m_cell_face_is_reversed;
NodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
EdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
FaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces;
CellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells;
EdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges;
NodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes;
CellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells;
FaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces;
NodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes;
CellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells;
EdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges;
FaceValuePerNode<const unsigned short> m_node_local_numbers_in_their_faces;
ConnectivityComputer m_connectivity_computer;
std::vector<RefFaceList> m_ref_face_list;
std::vector<RefNodeList> m_ref_node_list;
CellValue<const double> m_inv_cell_nb_nodes;
using Face = ConnectivityFace<Dimension>;
std::unordered_map<Face, FaceId, typename Face::Hash> m_face_number_map;
void _computeCellFaceAndFaceNodeConnectivities();
template <typename SubItemValuePerItemType>
PASTIS_INLINE
const SubItemValuePerItemType&
_lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const
{
if (not sub_item_value_per_item.isBuilt()) {
const_cast<SubItemValuePerItemType&>(sub_item_value_per_item)
= m_connectivity_computer
. computeLocalItemNumberInChildItem<SubItemValuePerItemType::item_t,
SubItemValuePerItemType::sub_item_t>(*this);
}
return sub_item_value_per_item;
}
friend class ConnectivityComputer;
PASTIS_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:
PASTIS_INLINE
const CellValue<const CellType>& cellType() const
{
return m_cell_type;
};
PASTIS_INLINE
const CellValue<const int>& cellNumber() const
{
return m_cell_number;
};
PASTIS_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>
PASTIS_INLINE
const auto getItemToItemMatrix() const
{
return ItemToItemMatrix<source_item_type,
target_item_type>(_getMatrix(source_item_type,
target_item_type));
}
PASTIS_INLINE
const auto cellToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
}
PASTIS_INLINE
const auto cellToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
}
PASTIS_INLINE
const auto cellToNodeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
}
PASTIS_INLINE
const auto faceToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
}
PASTIS_INLINE
const auto faceToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
}
PASTIS_INLINE
const auto faceToNodeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
}
PASTIS_INLINE
const auto edgeToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
}
PASTIS_INLINE
const auto edgeToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
}
PASTIS_INLINE
const auto edgeToNodeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
}
PASTIS_INLINE
const auto nodeToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
}
PASTIS_INLINE
const auto nodeToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
}
PASTIS_INLINE
const auto nodeToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>();
}
PASTIS_INLINE
const auto& cellFaceIsReversed() const
{ static_assert(dimension>1, "reversed faces makes no sense in dimension 1");
return m_cell_face_is_reversed;
}
PASTIS_INLINE
const auto& cellLocalNumbersInTheirNodes() const
{
return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes);
}
PASTIS_INLINE
const auto& cellLocalNumbersInTheirEdges() const
{
if constexpr (dimension>2) {
return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges);
} else {
return cellLocalNumbersInTheirFaces();
}
}
PASTIS_INLINE
const auto& cellLocalNumbersInTheirFaces() const
{
if constexpr (dimension>1) {
return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces);
} else {
return cellLocalNumbersInTheirNodes();
}
}
PASTIS_INLINE
const auto& faceLocalNumbersInTheirCells() const
{
if constexpr(dimension>1) {
return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells);
} else {
return nodeLocalNumbersInTheirCells();
}
}
PASTIS_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);
}
PASTIS_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);
}
PASTIS_INLINE
const auto& edgeLocalNumbersInTheirCells() const
{
if constexpr (dimension>2) {
return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells);
} else {
return faceLocalNumbersInTheirCells();
}
}
PASTIS_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);
}
PASTIS_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);
}
PASTIS_INLINE
const auto& nodeLocalNumbersInTheirCells() const
{
return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells);
}
PASTIS_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);
}
PASTIS_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);
}
void addRefFaceList(const RefFaceList& ref_face_list)
{
m_ref_face_list.push_back(ref_face_list);
}
size_t numberOfRefFaceList() const
{
return m_ref_face_list.size();
}
const RefFaceList& refFaceList(const size_t& i) const
{
return m_ref_face_list[i];
}
void addRefNodeList(const RefNodeList& ref_node_list)
{
m_ref_node_list.push_back(ref_node_list);
}
size_t numberOfRefNodeList() const
{
return m_ref_node_list.size();
}
const RefNodeList& refNodeList(const size_t& i) const
{
return m_ref_node_list[i];
}
PASTIS_INLINE
CSRGraph cellToCellGraph() const
{
std::vector<std::set<int>> cell_cells(this->numberOfCells());
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() == 2) {
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);
}
}
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);
}
PASTIS_INLINE
size_t numberOfNodes() const final
{
const auto& node_to_cell_matrix
= this->_getMatrix(ItemType::node,ItemType::cell);
return node_to_cell_matrix.numRows();
}
PASTIS_INLINE
size_t numberOfEdges() const final
{
const auto& edge_to_node_matrix
= this->_getMatrix(ItemType::edge,ItemType::node);
return edge_to_node_matrix.numRows();
}
PASTIS_INLINE
size_t numberOfFaces() const final
{
const auto& face_to_node_matrix
= this->_getMatrix(ItemType::face,ItemType::cell);
return face_to_node_matrix.numRows();
}
PASTIS_INLINE
size_t numberOfCells() const final
{
const auto& cell_to_node_matrix
= this->_getMatrix(ItemType::cell,ItemType::node);
return cell_to_node_matrix.numRows();
}
const CellValue<const double>& invCellNbNodes() const
{
#warning add calculation on demand when variables will be defined
return m_inv_cell_nb_nodes;
}
unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const
{
const Face face(face_nodes);
auto i_face = m_face_number_map.find(face);
if (i_face == m_face_number_map.end()) {
perr() << "Face " << face << " not found!\n";
throw std::exception();
std::exit(0); }
return i_face->second;
}
Connectivity(const Connectivity&) = delete;
Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
const std::vector<CellType>& cell_type_vector,
const std::vector<int>& cell_number_vector);
~Connectivity()
{
;
}
};
using Connectivity3D = Connectivity<3>;
using Connectivity2D = Connectivity<2>;
using Connectivity1D = Connectivity<1>;
#endif // CONNECTIVITY_HPP