Select Git revision
Polynomial1D.hpp
Connectivity.hpp 19.60 KiB
#ifndef CONNECTIVITY_HPP
#define CONNECTIVITY_HPP
#include <algebra/TinyVector.hpp>
#include <mesh/CellType.hpp>
#include <mesh/ConnectivityComputer.hpp>
#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/IConnectivity.hpp>
#include <mesh/ItemToItemMatrix.hpp>
#include <mesh/ItemType.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/RefId.hpp>
#include <mesh/RefItemList.hpp>
#include <mesh/SubItemValuePerItem.hpp>
#include <utils/CRSGraph.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <utils/PugsTraits.hpp>
#include <utils/PugsUtils.hpp>
#include <algorithm>
#include <iostream>
#include <set>
#include <vector>
class ConnectivityDescriptor;
template <size_t Dim>
class Connectivity final : public IConnectivity
{
private:
std::ostream& _write(std::ostream&) const final;
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:
size_t m_number_of_cells;
size_t m_number_of_faces;
size_t m_number_of_edges;
size_t m_number_of_nodes;
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;
WeakFaceValue<const bool> m_is_boundary_face;
WeakEdgeValue<const int> m_edge_owner;
WeakEdgeValue<const bool> m_edge_is_owned;
WeakEdgeValue<const bool> m_is_boundary_edge;
WeakNodeValue<const int> m_node_owner;
WeakNodeValue<const bool> m_node_is_owned;
WeakNodeValue<const bool> m_is_boundary_node;
WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;
WeakFaceValuePerCell<const uint16_t> m_cell_local_numbers_in_their_faces;
WeakEdgeValuePerCell<const uint16_t> m_cell_local_numbers_in_their_edges;
WeakNodeValuePerCell<const uint16_t> m_cell_local_numbers_in_their_nodes;
WeakCellValuePerFace<const uint16_t> m_face_local_numbers_in_their_cells;
WeakEdgeValuePerFace<const uint16_t> m_face_local_numbers_in_their_edges;
WeakNodeValuePerFace<const uint16_t> m_face_local_numbers_in_their_nodes;
WeakCellValuePerEdge<const uint16_t> m_edge_local_numbers_in_their_cells;
WeakFaceValuePerEdge<const uint16_t> m_edge_local_numbers_in_their_faces;
WeakNodeValuePerEdge<const uint16_t> m_edge_local_numbers_in_their_nodes;
WeakCellValuePerNode<const uint16_t> m_node_local_numbers_in_their_cells;
WeakFaceValuePerNode<const uint16_t> m_node_local_numbers_in_their_faces;
WeakEdgeValuePerNode<const uint16_t> m_node_local_numbers_in_their_edges;
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;
void _computeCellFaceAndFaceNodeConnectivities();
void _buildIsBoundaryFace() const;
void _buildIsBoundaryEdge() const;
void _buildIsBoundaryNode() const;
friend class ConnectivityComputer;
public:
PUGS_INLINE
const ConnectivityMatrix&
getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const final
{
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.computeInverseConnectivityMatrix(*this, item_type_0, item_type_1);
}
return connectivity_matrix;
}
PUGS_INLINE
CellValue<const CellType>
cellType() const
{
return m_cell_type;
}
PUGS_INLINE
CellValue<const int>
cellNumber() const
{ return m_cell_number;
}
PUGS_INLINE
FaceValue<const int>
faceNumber() const
{
return m_face_number;
}
PUGS_INLINE
EdgeValue<const int>
edgeNumber() const
{
return m_edge_number;
}
PUGS_INLINE
NodeValue<const int>
nodeNumber() const
{
return m_node_number;
}
template <ItemType item_type>
PUGS_INLINE 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
CellValue<const int>
cellOwner() const
{
return m_cell_owner;
}
PUGS_INLINE
FaceValue<const int>
faceOwner() const
{
return m_face_owner;
}
PUGS_INLINE
EdgeValue<const int>
edgeOwner() const
{
return m_edge_owner;
}
PUGS_INLINE
NodeValue<const int>
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
CellValue<const bool>
cellIsOwned() const
{
return m_cell_is_owned;
}
PUGS_INLINE
FaceValue<const bool>
faceIsOwned() const
{
return m_face_is_owned;
}
PUGS_INLINE
EdgeValue<const bool>
edgeIsOwned() const
{
return m_edge_is_owned;
}
PUGS_INLINE
NodeValue<const bool>
nodeIsOwned() const
{
return m_node_is_owned;
}
template <ItemType item_type>
PUGS_INLINE 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
FaceValue<const bool>
isBoundaryFace() const
{
if (not m_is_boundary_face.isBuilt()) {
this->_buildIsBoundaryFace();
}
return m_is_boundary_face;
}
PUGS_INLINE
EdgeValue<const bool>
isBoundaryEdge() const
{
if (not m_is_boundary_edge.isBuilt()) {
this->_buildIsBoundaryEdge();
}
return m_is_boundary_edge;
}
PUGS_INLINE
NodeValue<const bool>
isBoundaryNode() const
{
if (not m_is_boundary_node.isBuilt()) {
this->_buildIsBoundaryNode();
}
return m_is_boundary_node;
}
template <ItemType item_type>
PUGS_INLINE auto
isBoundary() const
{
if constexpr (item_type == ItemType::face) {
return isBoundaryFace();
} else if constexpr (item_type == ItemType::edge) {
return isBoundaryEdge();
} else if constexpr (item_type == ItemType::node) {
return isBoundaryNode();
} else {
static_assert(item_type != ItemType::cell, "cell boundary makes no sense");
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
PUGS_INLINE
const bool&
isConnectivityMatrixBuilt(ItemType item_type_0, 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 ItemToItemMatrix<source_item_type, target_item_type>
getItemToItemMatrix() const
{
return ItemToItemMatrix<source_item_type, target_item_type>(getMatrix(source_item_type, target_item_type));
}
PUGS_INLINE
ItemToItemMatrix<ItemType::cell, ItemType::face>
cellToFaceMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::cell, ItemType::edge>
cellToEdgeMatrix() const
{
return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::cell, ItemType::node>
cellToNodeMatrix() const
{ return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::face, ItemType::cell>
faceToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::face, ItemType::edge>
faceToEdgeMatrix() const
{
static_assert(Dimension > 2, "face to edge matrix makes sense in dimension > 2");
return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::face, ItemType::node>
faceToNodeMatrix() const
{
static_assert(Dimension > 1, "face to node matrix makes sense in dimension > 1");
return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::edge, ItemType::cell>
edgeToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::edge, ItemType::face>
edgeToFaceMatrix() const
{
static_assert(Dimension > 2, "edge to face matrix makes sense in dimension > 2");
return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::edge, ItemType::node>
edgeToNodeMatrix() const
{
static_assert(Dimension > 1, "edge to node matrix makes sense in dimension > 1");
return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::node, ItemType::cell>
nodeToCellMatrix() const
{
return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::node, ItemType::face>
nodeToFaceMatrix() const
{
static_assert(Dimension > 1, "node to face matrix makes sense in dimension > 1");
return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
}
PUGS_INLINE
ItemToItemMatrix<ItemType::node, ItemType::edge>
nodeToEdgeMatrix() const
{
static_assert(Dimension > 1, "node to edge matrix makes sense in dimension > 1");
return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>(); }
PUGS_INLINE
FaceValuePerCell<const bool>
cellFaceIsReversed() const
{
static_assert(Dimension > 1, "reversed faces makes no sense in dimension 1");
return m_cell_face_is_reversed;
}
PUGS_INLINE
EdgeValuePerFace<const bool>
faceEdgeIsReversed() const
{
static_assert(Dimension > 2, "reversed edges makes no sense in dimension 1 or 2");
return m_face_edge_is_reversed;
}
PUGS_INLINE
FaceValuePerCell<const uint16_t>
cellLocalNumbersInTheirFaces() const
{
if (not m_cell_local_numbers_in_their_faces.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<FaceOfCell>(*this);
}
return m_cell_local_numbers_in_their_faces;
}
PUGS_INLINE
EdgeValuePerCell<const uint16_t>
cellLocalNumbersInTheirEdges() const
{
if (not m_cell_local_numbers_in_their_edges.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<EdgeOfCell>(*this);
}
return m_cell_local_numbers_in_their_edges;
}
PUGS_INLINE
NodeValuePerCell<const uint16_t>
cellLocalNumbersInTheirNodes() const
{
if (not m_cell_local_numbers_in_their_nodes.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<NodeOfCell>(*this);
}
return m_cell_local_numbers_in_their_nodes;
}
PUGS_INLINE
CellValuePerFace<const uint16_t>
faceLocalNumbersInTheirCells() const
{
if (not m_face_local_numbers_in_their_cells.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<CellOfFace>(*this);
}
return m_face_local_numbers_in_their_cells;
}
PUGS_INLINE
EdgeValuePerFace<const uint16_t>
faceLocalNumbersInTheirEdges() const
{
static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
if (not m_face_local_numbers_in_their_edges.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<EdgeOfFace>(*this);
}
return m_face_local_numbers_in_their_edges;
}
PUGS_INLINE NodeValuePerFace<const uint16_t>
faceLocalNumbersInTheirNodes() const
{
static_assert(Dimension > 1, "this function has no meaning in 1d");
if (not m_face_local_numbers_in_their_nodes.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<NodeOfFace>(*this);
}
return m_face_local_numbers_in_their_nodes;
}
PUGS_INLINE
CellValuePerEdge<const uint16_t>
edgeLocalNumbersInTheirCells() const
{
if (not m_edge_local_numbers_in_their_cells.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<CellOfEdge>(*this);
}
return m_edge_local_numbers_in_their_cells;
}
PUGS_INLINE
FaceValuePerEdge<const uint16_t>
edgeLocalNumbersInTheirFaces() const
{
static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
if (not m_edge_local_numbers_in_their_faces.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<FaceOfEdge>(*this);
}
return m_edge_local_numbers_in_their_faces;
}
PUGS_INLINE
NodeValuePerEdge<const uint16_t>
edgeLocalNumbersInTheirNodes() const
{
static_assert(Dimension > 1, "this function has no meaning in 1d");
if (not m_edge_local_numbers_in_their_nodes.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<NodeOfEdge>(*this);
}
return m_edge_local_numbers_in_their_nodes;
}
PUGS_INLINE
CellValuePerNode<const uint16_t>
nodeLocalNumbersInTheirCells() const
{
if (not m_node_local_numbers_in_their_cells.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<CellOfNode>(*this);
}
return m_node_local_numbers_in_their_cells;
}
PUGS_INLINE
FaceValuePerNode<const uint16_t>
nodeLocalNumbersInTheirFaces() const
{
static_assert(Dimension > 1, "this function has no meaning in 1d");
if (not m_node_local_numbers_in_their_faces.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<FaceOfNode>(*this);
}
return m_node_local_numbers_in_their_faces;
}
PUGS_INLINE
EdgeValuePerNode<const uint16_t>
nodeLocalNumbersInTheirEdges() const
{
static_assert(Dimension > 1, "this function has no meaning in 1d or 2d");
if (not m_node_local_numbers_in_their_edges.isBuilt()) {
m_connectivity_computer.computeLocalItemNumberInChildItem<EdgeOfNode>(*this); }
return m_node_local_numbers_in_their_edges;
}
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(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
CRSGraph
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() > 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); }
}
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 CRSGraph(entries, neighbors);
}
PUGS_INLINE
size_t
numberOfNodes() const final
{
return m_number_of_nodes;
}
PUGS_INLINE
size_t
numberOfEdges() const final
{
return m_number_of_edges;
}
PUGS_INLINE
size_t
numberOfFaces() const final
{
return m_number_of_faces;
}
PUGS_INLINE
size_t
numberOfCells() const final
{
return m_number_of_cells;
}
template <ItemType item_type>
PUGS_INLINE size_t
numberOf() const
{
if constexpr (item_type == ItemType::cell) {
return this->numberOfCells();
} else if constexpr (item_type == ItemType::face) {
return this->numberOfFaces();
} else if constexpr (item_type == ItemType::edge) {
return this->numberOfEdges();
} else if constexpr (item_type == ItemType::node) {
return this->numberOfNodes();
} else {
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
Connectivity(const Connectivity&) = delete;
private:
Connectivity(); void _buildFrom(const ConnectivityDescriptor& descriptor);
public:
~Connectivity() = default;
};
template <size_t Dimension>
PUGS_INLINE std::shared_ptr<Connectivity<Dimension>>
Connectivity<Dimension>::build(const ConnectivityDescriptor& descriptor)
{
// Cannot use std::make_shared in this context since the default constructor is private
std::shared_ptr<Connectivity<Dimension>> connectivity_ptr(new Connectivity<Dimension>);
connectivity_ptr->_buildFrom(descriptor);
return connectivity_ptr;
}
using Connectivity3D = Connectivity<3>;
using Connectivity2D = Connectivity<2>;
using Connectivity1D = Connectivity<1>;
#endif // CONNECTIVITY_HPP