2477fe2de044255ecdb38d962d73a036731f34fb to 7843720c2e48bfcc1999236713db56570aca3433 · code / pugs

CMakeLists.txt

+3
−2

Original line number
Diff line number
Diff line

@@ -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)

)

src/algebra/TinyMatrix.hpp

+8
−8

Original line number
Diff line number
Diff line

@@ -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);

src/main.cpp

+6
−6

Original line number
Diff line number
Diff line

@@ -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

src/mesh/Connectivity.cpp

+53
−23

Original line number
Diff line number
Diff line

@@ -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>();

    }

  }

}

Compare revisions

Source

Target

Files

CMakeLists.txt

src/algebra/TinyMatrix.hpp

src/main.cpp

src/mesh/Connectivity.cpp