Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
Loading items

Target

Select target project
  • code/pugs
1 result
Select Git revision
Loading items
Show changes
Commits on Source (19)
Showing
with 1206 additions and 707 deletions
......@@ -271,8 +271,9 @@ add_executable(
# Libraries
target_link_libraries(
pastis
PastisMesh
PastisUtils
kokkos
${PARMETIS_LIBRARIES}
${MPI_CXX_LINK_FLAGS} ${MPI_CXX_LIBRARIES}
PastisMesh
PastisUtils)
)
......@@ -296,7 +296,7 @@ template <size_t N, typename T>
PASTIS_INLINE
constexpr T det(const TinyMatrix<N,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "determinent for arbitrary dimension N is defined for floating point types only");
TinyMatrix<N,T> M = A;
......@@ -340,7 +340,7 @@ template <typename T>
PASTIS_INLINE
constexpr T det(const TinyMatrix<1,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0,0);
}
......@@ -348,7 +348,7 @@ template <typename T>
PASTIS_INLINE
constexpr T det(const TinyMatrix<2,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return A(0,0)*A(1,1)-A(1,0)*A(0,1);
}
......@@ -356,7 +356,7 @@ template <typename T>
PASTIS_INLINE
constexpr T det(const TinyMatrix<3,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non-arithmetic types");
return
A(0,0)*(A(1,1)*A(2,2)-A(2,1)*A(1,2))
-A(1,0)*(A(0,1)*A(2,2)-A(2,1)*A(0,2))
......@@ -399,7 +399,7 @@ template <typename T>
PASTIS_INLINE
constexpr TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
TinyMatrix<1,T> A_1(1./A(0,0));
......@@ -412,7 +412,7 @@ constexpr T cofactor(const TinyMatrix<N,T>& A,
const size_t& i,
const size_t& j)
{
static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "cofactor is not defined for non-arithmetic types");
const T sign = ((i+j)%2) ? -1 : 1;
return sign * det(getMinor(A, i, j));
......@@ -422,7 +422,7 @@ template <typename T>
PASTIS_INLINE
constexpr TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
const T determinent = det(A);
......@@ -438,7 +438,7 @@ template <typename T>
PASTIS_INLINE
constexpr TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A)
{
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non-arithmetic types");
static_assert(std::is_floating_point<T>::value, "inverse is defined for floating point types only");
const T determinent = det(A);
......
......@@ -74,9 +74,9 @@ int main(int argc, char *argv[])
case BoundaryConditionDescriptor::Type::symmetry: {
const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
= dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList();
for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefItemList<ItemType::node>();
++i_ref_node_list) {
const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list);
const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
const RefId& ref = ref_node_list.refId();
if (ref == sym_bc_descriptor.boundaryDescriptor()) {
SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
......@@ -191,9 +191,9 @@ int main(int argc, char *argv[])
case BoundaryConditionDescriptor::Type::symmetry: {
const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
= dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefItemList<ItemType::face>();
++i_ref_face_list) {
const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
const RefId& ref = ref_face_list.refId();
if (ref == sym_bc_descriptor.boundaryDescriptor()) {
SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
......@@ -295,9 +295,9 @@ int main(int argc, char *argv[])
case BoundaryConditionDescriptor::Type::symmetry: {
const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
= dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefItemList<ItemType::face>();
++i_ref_face_list) {
const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
const RefId& ref = ref_face_list.refId();
if (ref == sym_bc_descriptor.boundaryDescriptor()) {
SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc
......
......@@ -3,6 +3,8 @@
#include <Messenger.hpp>
#include <ConnectivityDescriptor.hpp>
template<size_t Dimension>
Connectivity<Dimension>::Connectivity() {}
......@@ -10,8 +12,7 @@ template<size_t Dimension>
void Connectivity<Dimension>::
_buildFrom(const ConnectivityDescriptor& descriptor)
{
#warning All of these should be checked by ConnectivityDescriptor
Assert(descriptor.cell_by_node_vector.size() == descriptor.cell_type_vector.size());
Assert(descriptor.cell_to_node_vector.size() == descriptor.cell_type_vector.size());
Assert(descriptor.cell_number_vector.size() == descriptor.cell_type_vector.size());
if constexpr (Dimension>1) {
Assert(descriptor.cell_to_face_vector.size() == descriptor.cell_type_vector.size());
......@@ -21,7 +22,7 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
auto& cell_to_node_matrix
= m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
cell_to_node_matrix = descriptor.cell_by_node_vector;
cell_to_node_matrix = descriptor.cell_to_node_vector;
{
WeakCellValue<CellType> cell_type(*this);
......@@ -45,8 +46,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
{
WeakCellValue<int> cell_global_index(*this);
#warning index must start accounting number of global indices of other procs
#warning must take care of ghost cells
int first_index = 0;
parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
cell_global_index[j] = first_index+j;
......@@ -54,16 +53,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
m_cell_global_index = cell_global_index;
}
{
WeakCellValue<double> inv_cell_nb_nodes(*this);
parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
inv_cell_nb_nodes[j] = 1./cell_nodes.length;
});
m_inv_cell_nb_nodes = inv_cell_nb_nodes;
}
{
WeakCellValue<int> cell_owner(*this);
cell_owner = convert_to_array(descriptor.cell_owner_vector);
......@@ -94,14 +83,13 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
m_node_is_owned = node_is_owned;
}
m_ref_node_list_vector = descriptor.template refItemListVector<ItemType::node>();
m_ref_cell_list_vector = descriptor.template refItemListVector<ItemType::cell>();
if constexpr (Dimension>1) {
auto& face_to_node_matrix
= m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
face_to_node_matrix = descriptor.face_to_node_vector;
m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;
auto& cell_to_face_matrix
= m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
cell_to_face_matrix = descriptor.cell_to_face_vector;
m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;
{
FaceValuePerCell<bool> cell_face_is_reversed(*this);
......@@ -111,14 +99,15 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
cell_face_is_reversed(j,lj) = face_cells_vector[lj];
}
}
m_cell_face_is_reversed = cell_face_is_reversed;
}
{
WeakFaceValue<int> face_number(*this);
face_number = convert_to_array(descriptor.face_number_vector);
m_face_number = face_number;
}
{
WeakFaceValue<int> face_owner(*this);
face_owner = convert_to_array(descriptor.face_owner_vector);
......@@ -134,8 +123,49 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
m_face_is_owned = face_is_owned;
}
m_ref_face_list_vector = descriptor.template refItemListVector<ItemType::face>();
if constexpr (Dimension>2) {
m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;
m_ref_face_list = descriptor.refFaceList();
m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_vector;
m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_vector;
{
EdgeValuePerFace<bool> face_edge_is_reversed(*this);
for (FaceId l=0; l<descriptor.face_edge_is_reversed_vector.size(); ++l) {
const auto& edge_faces_vector = descriptor.face_edge_is_reversed_vector[l];
for (unsigned short el=0; el<edge_faces_vector.size(); ++el) {
face_edge_is_reversed(l,el) = edge_faces_vector[el];
}
}
m_face_edge_is_reversed = face_edge_is_reversed;
}
{
WeakEdgeValue<int> edge_number(*this);
edge_number = convert_to_array(descriptor.edge_number_vector);
m_edge_number = edge_number;
}
{
WeakEdgeValue<int> edge_owner(*this);
edge_owner = convert_to_array(descriptor.edge_owner_vector);
m_edge_owner = edge_owner;
}
{
const int rank = parallel::rank();
WeakEdgeValue<bool> edge_is_owned(*this);
parallel_for(this->numberOfEdges(), PASTIS_LAMBDA(const EdgeId& e) {
edge_is_owned[e] = (m_edge_owner[e] == rank);
});
m_edge_is_owned = edge_is_owned;
}
m_ref_edge_list_vector = descriptor.template refItemListVector<ItemType::edge>();
}
}
}
......
......@@ -7,6 +7,8 @@
#include <PastisOStream.hpp>
#include <PastisUtils.hpp>
#include <PastisTraits.hpp>
#include <TinyVector.hpp>
#include <ItemValue.hpp>
......@@ -28,58 +30,14 @@
#include <RefId.hpp>
#include <ItemType.hpp>
#include <RefNodeList.hpp>
#include <RefFaceList.hpp>
#include <RefItemList.hpp>
#include <SynchronizerManager.hpp>
#include <tuple>
#include <algorithm>
#include <set>
class ConnectivityDescriptor
{
private:
std::vector<RefFaceList> m_ref_face_list;
public:
std::vector<std::vector<unsigned int>> cell_by_node_vector;
#warning should be set in 2d
std::vector<std::vector<unsigned int>> cell_to_face_vector;
std::vector<std::vector<bool>> cell_face_is_reversed_vector;
// std::vector<std::vector<unsigned int>> face_to_cell_vector;
std::vector<std::vector<unsigned int>> face_to_node_vector;
std::vector<CellType> cell_type_vector;
std::vector<int> cell_number_vector;
std::vector<int> face_number_vector;
std::vector<int> node_number_vector;
std::vector<int> cell_owner_vector;
std::vector<int> face_owner_vector;
std::vector<int> node_owner_vector;
const std::vector<RefFaceList>& refFaceList() const
{
return m_ref_face_list;
}
void addRefFaceList(const RefFaceList& ref_face_list)
{
m_ref_face_list.push_back(ref_face_list);
}
ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
ConnectivityDescriptor() = default;
ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
ConnectivityDescriptor(ConnectivityDescriptor&&) = delete;
~ConnectivityDescriptor() = default;
};
class ConnectivityDescriptor;
template <size_t Dim>
class Connectivity final
......@@ -105,12 +63,10 @@ class Connectivity final
ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
WeakCellValue<const CellType> m_cell_type;
#warning is m_cell_global_index really necessary? should it be computed on demand instead?
WeakCellValue<const int> m_cell_global_index;
WeakCellValue<const int> m_cell_number;
WeakFaceValue<const int> m_face_number;
#warning check that m_edge_number is filled correctly
WeakEdgeValue<const int> m_edge_number;
WeakNodeValue<const int> m_node_number;
......@@ -127,6 +83,7 @@ class Connectivity final
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;
......@@ -146,8 +103,10 @@ class Connectivity final
ConnectivityComputer m_connectivity_computer;
std::vector<RefFaceList> m_ref_face_list;
std::vector<RefNodeList> m_ref_node_list;
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;
......@@ -228,8 +187,7 @@ class Connectivity final
} else if constexpr(item_type == ItemType::node) {
return m_node_number;
} else {
static_assert(item_type == ItemType::cell, "unknown ItemType");
return m_cell_number;
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
......@@ -272,8 +230,7 @@ class Connectivity final
} else if constexpr(item_type == ItemType::node) {
return m_node_owner;
} else {
static_assert(item_type == ItemType::cell, "unknown ItemType");
return m_cell_owner;
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
......@@ -316,8 +273,7 @@ class Connectivity final
} else if constexpr(item_type == ItemType::node) {
return m_node_is_owned;
} else {
static_assert(item_type == ItemType::cell, "unknown ItemType");
return m_cell_is_owned;
static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
}
}
......@@ -418,6 +374,13 @@ class Connectivity final
return m_cell_face_is_reversed;
}
PASTIS_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;
}
PASTIS_INLINE
const auto& cellLocalNumbersInTheirNodes() const
{
......@@ -512,29 +475,53 @@ class Connectivity final
return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
}
size_t numberOfRefFaceList() const
template <ItemType item_type>
size_t numberOfRefItemList() const
{
return m_ref_face_list.size();
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");
}
const RefFaceList& refFaceList(const size_t& i) const
{
return m_ref_face_list[i];
}
void addRefNodeList(const RefNodeList& ref_node_list)
template <ItemType item_type>
const RefItemList<item_type>& refItemList(const size_t& i) const
{
m_ref_node_list.push_back(ref_node_list);
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");
}
size_t numberOfRefNodeList() const
{
return m_ref_node_list.size();
}
const RefNodeList& refNodeList(const size_t& i) const
template <ItemType item_type>
void addRefItemList(const RefItemList<item_type>& ref_item_list)
{
return m_ref_node_list[i];
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");
}
}
PASTIS_INLINE
......@@ -624,7 +611,18 @@ class Connectivity final
CellValue<const double> invCellNbNodes() const
{
#warning add calculation on demand when variables will be defined
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(), PASTIS_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;
}
......
#ifndef CONNECTIVITY_DESCRIPTOR_HPP
#define CONNECTIVITY_DESCRIPTOR_HPP
#include <RefItemList.hpp>
#include <PastisTraits.hpp>
#include <ItemOfItemType.hpp>
#include <vector>
class ConnectivityDescriptor
{
private:
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;
public:
std::vector<std::vector<unsigned int>> cell_to_node_vector;
std::vector<std::vector<unsigned int>> cell_to_face_vector;
std::vector<std::vector<unsigned int>> cell_to_edge_vector;
std::vector<std::vector<unsigned int>> face_to_node_vector;
std::vector<std::vector<unsigned int>> face_to_edge_vector;
std::vector<std::vector<unsigned int>> edge_to_node_vector;
template <typename ItemOfItemT>
auto& itemOfItemVector()
{
if constexpr (std::is_same_v<ItemOfItemT,NodeOfCell>) {
return cell_to_node_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,FaceOfCell>) {
return cell_to_face_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfCell>) {
return cell_to_edge_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfFace>) {
return face_to_edge_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,NodeOfFace>) {
return face_to_node_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,NodeOfEdge>) {
return edge_to_node_vector;
} else {
static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type");
}
}
std::vector<Array<bool>> cell_face_is_reversed_vector;
std::vector<Array<bool>> face_edge_is_reversed_vector;
std::vector<CellType> cell_type_vector;
std::vector<int> cell_number_vector;
std::vector<int> face_number_vector;
std::vector<int> edge_number_vector;
std::vector<int> node_number_vector;
template <ItemType item_type>
const std::vector<int>& itemNumberVector() const
{
if constexpr (item_type == ItemType::cell) {
return cell_number_vector;
} else if constexpr (item_type == ItemType::face) {
return face_number_vector;
} else if constexpr (item_type == ItemType::edge) {
return edge_number_vector;
} else if constexpr (item_type == ItemType::node) {
return node_number_vector;
} else {
static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
}
}
std::vector<int> cell_owner_vector;
std::vector<int> face_owner_vector;
std::vector<int> edge_owner_vector;
std::vector<int> node_owner_vector;
template <ItemType item_type>
const std::vector<RefItemList<item_type>>& refItemListVector() const
{
if constexpr (item_type == ItemType::cell) {
return m_ref_cell_list_vector;
} else if constexpr (item_type == ItemType::face) {
return m_ref_face_list_vector;
} else if constexpr (item_type == ItemType::edge) {
return m_ref_edge_list_vector;
} else if constexpr (item_type == ItemType::node) {
return m_ref_node_list_vector;
} 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");
}
}
ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
ConnectivityDescriptor() = default;
ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
ConnectivityDescriptor(ConnectivityDescriptor&&) = delete;
~ConnectivityDescriptor() = default;
};
#endif // CONNECTIVITY_DESCRIPTOR_HPP
This diff is collapsed.
......@@ -6,6 +6,7 @@
#include <ItemValueUtils.hpp>
#include <unordered_map>
#include <ConnectivityDescriptor.hpp>
template <int Dimension>
class ConnectivityDispatcher
......@@ -23,8 +24,8 @@ class ConnectivityDispatcher
{
using ItemId = ItemIdT<item_type>;
ItemValue<const int, item_type> m_new_owner;
Array<const unsigned int> m_list_to_send_size_by_proc;
std::vector<Array<const ItemId>> m_list_to_send_by_proc;
#warning is m_list_to_recv_size_by_proc really necessary?
Array<const unsigned int> m_list_to_recv_size_by_proc;
std::unordered_map<int, int> m_number_to_id_map;
std::vector<Array<const ItemId>> m_recv_id_correspondance_by_proc;
......@@ -65,36 +66,138 @@ class ConnectivityDispatcher
}
}
template <typename ItemToItem>
struct DispatchedItemOfItemInfo
{
std::vector<Array<const int>> m_number_of_sub_item_per_item_to_recv_by_proc;
std::vector<Array<const int>> m_sub_item_numbers_to_recv_by_proc;
};
DispatchedItemOfItemInfo<NodeOfCell> m_dispatched_node_of_cell_info;
DispatchedItemOfItemInfo<EdgeOfCell> m_dispatched_edge_of_cell_info;
DispatchedItemOfItemInfo<FaceOfCell> m_dispatched_face_of_cell_info;
DispatchedItemOfItemInfo<NodeOfEdge> m_dispatched_node_of_edge_info;
DispatchedItemOfItemInfo<FaceOfEdge> m_dispatched_face_of_edge_info;
DispatchedItemOfItemInfo<CellOfEdge> m_dispatched_cell_of_edge_info;
DispatchedItemOfItemInfo<NodeOfFace> m_dispatched_node_of_face_info;
DispatchedItemOfItemInfo<EdgeOfFace> m_dispatched_edge_of_face_info;
DispatchedItemOfItemInfo<CellOfFace> m_dispatched_cell_of_face_info;
DispatchedItemOfItemInfo<EdgeOfNode> m_dispatched_edge_of_node_info;
DispatchedItemOfItemInfo<FaceOfNode> m_dispatched_face_of_node_info;
DispatchedItemOfItemInfo<CellOfNode> m_dispatched_cell_of_node_info;
template <typename ItemOfItem>
PASTIS_INLINE
DispatchedItemOfItemInfo<ItemOfItem>& _dispatchedInfo()
{
if constexpr (std::is_same_v<NodeOfCell, ItemOfItem>) {
return m_dispatched_node_of_cell_info;
} else if constexpr (std::is_same_v<EdgeOfCell, ItemOfItem>) {
return m_dispatched_edge_of_cell_info;
} else if constexpr (std::is_same_v<FaceOfCell, ItemOfItem>) {
return m_dispatched_face_of_cell_info;
} else if constexpr (std::is_same_v<NodeOfEdge, ItemOfItem>) {
return m_dispatched_node_of_edge_info;
} else if constexpr (std::is_same_v<FaceOfEdge, ItemOfItem>) {
return m_dispatched_face_of_edge_info;
} else if constexpr (std::is_same_v<CellOfEdge, ItemOfItem>) {
return m_dispatched_cell_of_edge_info;
} else if constexpr (std::is_same_v<NodeOfFace, ItemOfItem>) {
return m_dispatched_node_of_face_info;
} else if constexpr (std::is_same_v<EdgeOfFace, ItemOfItem>) {
return m_dispatched_edge_of_face_info;
} else if constexpr (std::is_same_v<CellOfFace, ItemOfItem>) {
return m_dispatched_cell_of_face_info;
} else if constexpr (std::is_same_v<EdgeOfNode, ItemOfItem>) {
return m_dispatched_edge_of_node_info;
} else if constexpr (std::is_same_v<FaceOfNode, ItemOfItem>) {
return m_dispatched_face_of_node_info;
} else if constexpr (std::is_same_v<CellOfNode, ItemOfItem>) {
return m_dispatched_cell_of_node_info;
} else {
static_assert(is_false_v<ItemOfItem>, "Unexpected ItemOfItem type");
}
}
template <typename ItemOfItem>
PASTIS_INLINE
const DispatchedItemOfItemInfo<ItemOfItem>& _dispatchedInfo() const
{
if constexpr (std::is_same_v<NodeOfCell, ItemOfItem>) {
return m_dispatched_node_of_cell_info;
} else if constexpr (std::is_same_v<EdgeOfCell, ItemOfItem>) {
return m_dispatched_edge_of_cell_info;
} else if constexpr (std::is_same_v<FaceOfCell, ItemOfItem>) {
return m_dispatched_face_of_cell_info;
} else if constexpr (std::is_same_v<NodeOfEdge, ItemOfItem>) {
return m_dispatched_node_of_edge_info;
} else if constexpr (std::is_same_v<FaceOfEdge, ItemOfItem>) {
return m_dispatched_face_of_edge_info;
} else if constexpr (std::is_same_v<CellOfEdge, ItemOfItem>) {
return m_dispatched_cell_of_edge_info;
} else if constexpr (std::is_same_v<NodeOfFace, ItemOfItem>) {
return m_dispatched_node_of_face_info;
} else if constexpr (std::is_same_v<EdgeOfFace, ItemOfItem>) {
return m_dispatched_edge_of_face_info;
} else if constexpr (std::is_same_v<CellOfFace, ItemOfItem>) {
return m_dispatched_cell_of_face_info;
} else if constexpr (std::is_same_v<EdgeOfNode, ItemOfItem>) {
return m_dispatched_edge_of_node_info;
} else if constexpr (std::is_same_v<FaceOfNode, ItemOfItem>) {
return m_dispatched_face_of_node_info;
} else if constexpr (std::is_same_v<CellOfNode, ItemOfItem>) {
return m_dispatched_cell_of_node_info;
} else {
static_assert(is_false_v<ItemOfItem>, "Unexpected ItemOfItem type");
}
}
template <ItemType item_type>
void _buildNewOwner();
template <ItemType item_type>
void _buildItemListToSend();
Array<const unsigned int> _buildNbCellToSend();
void _buildCellNumberIdMap();
template <typename ItemOfItemT>
void _buildSubItemNumberToIdMap(const std::vector<Array<const int>>& recv_cell_node_number_by_proc);
void _buildSubItemNumberToIdMap();
template <ItemType item_type>
void _buildItemToExchangeLists();
template <ItemType item_type>
void _buildNumberOfItemToExchange();
template <typename ItemOfItemT>
void _buildItemToSubItemDescriptor();
void _dispatchEdges();
void _dispatchFaces();
template<typename DataType, ItemType item_type, typename ConnectivityPtr>
void _gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
std::vector<std::remove_const_t<DataType>>& gathered_vector);
template<typename DataType, typename ItemOfItem, typename ConnectivityPtr>
void _gatherFrom(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
std::vector<Array<std::remove_const_t<DataType>>>& gathered_vector);
template <typename SubItemOfItemT>
std::vector<Array<const int>>
_getRecvNumberOfSubItemPerItemByProc();
void _buildNumberOfSubItemPerItemToRecvByProc();
template <typename SubItemOfItemT>
std::vector<Array<const int>>
_getRecvItemSubItemNumberingByProc(const std::vector<Array<const int>>& recv_number_of_sub_item_per_item_by_proc);
void _buildSubItemNumbersToRecvByProc();
template <ItemType item_type>
void _buildRecvItemIdCorrespondanceByProc();
template <ItemType item_type>
void _buildItemReferenceList();
public:
std::shared_ptr<const ConnectivityType>
dispatchedConnectivity() const
......
This diff is collapsed.
......@@ -262,6 +262,9 @@ private:
template <size_t Dimension>
void __computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
template <size_t Dimension>
void __computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
template <int Dimension>
void _dispatch();
......
......@@ -83,7 +83,7 @@ class ItemToItemMatrix
PASTIS_INLINE
const auto& operator[](const IndexType& source_id) const
{
static_assert(std::is_same<IndexType, SourceItemId>(),
static_assert(std::is_same_v<IndexType, SourceItemId>,
"ItemToItemMatrix must be indexed using correct ItemId");
using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
......
......@@ -118,4 +118,8 @@ struct ItemTypeId<3>
}
};
template <ItemType item_type>
PASTIS_INLINE
constexpr bool is_false_item_type_v = false;
#endif // ITEM_TYPE_HPP
......@@ -95,13 +95,10 @@ class ItemValue
}
template <typename IndexType>
DataType& operator[](const IndexType& i) const noexcept(NO_ASSERT)
DataType& operator[](const IndexType&) const noexcept(NO_ASSERT)
{
static_assert(std::is_same_v<IndexType,ItemId>,
"ItemValue must be indexed by ItemId");
static_assert(not std::is_const_v<DataType>,
"Cannot modify ItemValue of const");
return m_values[i];
}
PASTIS_INLINE
......
......@@ -11,7 +11,6 @@
struct IMesh
{
virtual size_t dimension() const = 0;
// virtual CSRGraph cellToCellGraph() const = 0;
~IMesh() = default;
};
......@@ -66,7 +65,6 @@ public:
return m_xr;
}
[[deprecated("should rework this class: quite ugly")]]
PASTIS_INLINE
NodeValue<Rd> mutableXr() const
{
......
......@@ -178,13 +178,11 @@ class MeshData
const FaceId& l = cell_faces[L];
const auto& face_nodes = face_to_node_matrix[l];
#warning should this lambda be replaced by a precomputed correspondance?
auto local_node_number_in_cell
= [&](const NodeId& node_number) {
for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
if (node_number == cell_nodes[i_node]) {
return i_node;
break;
}
}
return std::numeric_limits<size_t>::max();
......@@ -227,31 +225,37 @@ class MeshData
}
public:
PASTIS_INLINE
const MeshType& mesh() const
{
return m_mesh;
}
PASTIS_INLINE
const NodeValuePerCell<const Rd>& Cjr() const
{
return m_Cjr;
}
PASTIS_INLINE
const NodeValuePerCell<const double>& ljr() const
{
return m_ljr;
}
PASTIS_INLINE
const NodeValuePerCell<const Rd>& njr() const
{
return m_njr;
}
PASTIS_INLINE
const CellValue<const Rd>& xj() const
{
return m_xj;
}
PASTIS_INLINE
const CellValue<const double>& Vj() const
{
return m_Vj;
......
......@@ -7,8 +7,7 @@
#include <Kokkos_Vector.hpp>
#include <TinyVector.hpp>
#include <RefNodeList.hpp>
#include <RefFaceList.hpp>
#include <RefItemList.hpp>
#include <ConnectivityMatrix.hpp>
#include <IConnectivity.hpp>
......@@ -38,7 +37,7 @@ class MeshNodeBoundary
const auto& face_to_cell_matrix
= mesh.connectivity().faceToCellMatrix();
const Array<const FaceId>& face_list = ref_face_list.faceList();
const Array<const FaceId>& face_list = ref_face_list.list();
parallel_for(face_list.size(), PASTIS_LAMBDA(const int& l){
const auto& face_cells = face_to_cell_matrix[face_list[l]];
if (face_cells.size()>1) {
......@@ -74,7 +73,7 @@ class MeshNodeBoundary
template <typename MeshType>
MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
: m_node_list(ref_node_list.nodeList())
: m_node_list(ref_node_list.list())
{
static_assert(Dimension == MeshType::Dimension);
}
......@@ -288,7 +287,6 @@ _getNormal(const MeshType& mesh)
zmax = x;
}
}
#warning re work this part to avoir parallelism dependance
Array<R3> xmin_array = parallel::allGather(xmin);
Array<R3> xmax_array = parallel::allGather(xmax);
Array<R3> ymin_array = parallel::allGather(ymin);
......@@ -321,7 +319,6 @@ _getNormal(const MeshType& mesh)
if (x[2] > zmax[2]) { zmax = x; }
}
const R3 u = xmax-xmin;
const R3 v = ymax-ymin;
const R3 w = zmax-zmin;
......@@ -356,7 +353,6 @@ _getNormal(const MeshType& mesh)
normal *= 1./sqrt(normal_l2);
#warning Add flatness test
// this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
return normal;
......
#ifndef REF_FACE_LIST_HPP
#define REF_FACE_LIST_HPP
#include <Array.hpp>
#include <RefId.hpp>
class RefFaceList
{
private:
RefId m_ref_id;
Array<const FaceId> m_face_id_list;
public:
const RefId& refId() const
{
return m_ref_id;
}
const Array<const FaceId>& faceList() const
{
return m_face_id_list;
}
RefFaceList(const RefId& ref_id,
const Array<const FaceId>& face_id_list)
: m_ref_id(ref_id),
m_face_id_list(face_id_list)
{
;
}
RefFaceList& operator=(const RefFaceList&) = default;
RefFaceList& operator=(RefFaceList&&) = default;
RefFaceList() = default;
RefFaceList(const RefFaceList&) = default;
~RefFaceList() = default;
};
#endif // REF_FACE_LIST_HPP
#ifndef REF_ITEM_LIST_HPP
#define REF_ITEM_LIST_HPP
#include <Array.hpp>
#include <RefId.hpp>
#include <ItemId.hpp>
template <ItemType item_type>
class RefItemList
{
public:
using ItemId = ItemIdT<item_type>;
private:
RefId m_ref_id;
Array<const ItemId> m_item_id_list;
public:
const RefId& refId() const
{
return m_ref_id;
}
const Array<const ItemId>& list() const
{
return m_item_id_list;
}
RefItemList(const RefId& ref_id,
const Array<const ItemId>& item_id_list)
: m_ref_id(ref_id),
m_item_id_list(item_id_list)
{
;
}
RefItemList& operator=(const RefItemList&) = default;
RefItemList& operator=(RefItemList&&) = default;
RefItemList() = default;
RefItemList(const RefItemList&) = default;
~RefItemList() = default;
};
using RefNodeList = RefItemList<ItemType::node>;
using RefEdgeList = RefItemList<ItemType::edge>;
using RefFaceList = RefItemList<ItemType::face>;
using RefCellList = RefItemList<ItemType::cell>;
#endif // REF_ITEM_LIST_HPP
#ifndef REF_NODE_LIST_HPP
#define REF_NODE_LIST_HPP
#include <Array.hpp>
#include <RefId.hpp>
class RefNodeList
{
private:
RefId m_ref_id;
Array<const NodeId> m_node_id_list;
public:
const RefId& refId() const
{
return m_ref_id;
}
const Array<const NodeId>& nodeList() const
{
return m_node_id_list;
}
RefNodeList(const RefId& ref_id,
const Array<const NodeId>& node_id_list)
: m_ref_id(ref_id),
m_node_id_list(node_id_list)
{
;
}
RefNodeList& operator=(const RefNodeList&) = default;
RefNodeList& operator=(RefNodeList&&) = default;
RefNodeList() = default;
RefNodeList(const RefNodeList&) = default;
~RefNodeList() = default;
};
#endif // REF_NODE_LIST_HPP
......@@ -17,20 +17,10 @@
#include <memory>
// template <typename DataType,
// typename ItemOfItem,
// typename ConnectivityPtr = std::shared_ptr<const IConnectivity>,
// typename Allowed=void>
// class SubItemValuePerItem;
template <typename DataType,
typename ItemOfItem,
typename ConnectivityPtr = std::shared_ptr<const IConnectivity>>
class SubItemValuePerItem// <DataType,
// ItemOfItem,
// ConnectivityPtr// ,
// // std::enable_if_t<sub_item_type != item_type>
// >
class SubItemValuePerItem
{
public:
static constexpr ItemType item_type{ItemOfItem::item_type};
......