ftensor_matrix_block_views.cpp

Go to the documentation of this file.

/** @file ftensor_matrix_block_views.cpp
 * @brief Check rectangular form blocks, row offsets and column cursor strides.
 */
 
#define FTENSOR_DEBUG
#include <MoFEM.hpp>
#include <type_traits>
using namespace MoFEM;
 
namespace {
 
using FormBase =
    FormsIntegrators<VolumeElementForcesAndSourcesCore::UserDataOperator>::
        Assembly<AssemblyType::PETSC>::OpBase;
 
struct OpCheckBlocks : FormBase {
  OpCheckBlocks() : FormBase("row", "col", OPROWCOL) {}
 
  template <int RowDim, int ColDim> MoFEMErrorCode check() {
    MoFEMFunctionBegin;
    FTensor::Index<'i', RowDim> i;
    FTensor::Index<'j', ColDim> j;
    constexpr int row_offset = 1;
    constexpr int block_count = 3;
    locMat.resize(row_offset + 2 * RowDim + 1, block_count * ColDim + 1, false);
    locMat.clear();
    for (int row = 0; row != 2; ++row) {
      auto t_block = getLocMat<RowDim, ColDim>(row_offset + row * RowDim);
      for (int col = 0; col != block_count; ++col, ++t_block)
        t_block(i, j) = 10. * (row + 1) + col;
    }
    // Scalar inspection is deliberately independent of the tensor adapter.
    for (std::size_t row = 0; row != locMat.size1(); ++row) {
      for (std::size_t col = 0; col != locMat.size2(); ++col) {
        const double expected =
            row < row_offset || row >= row_offset + 2 * RowDim ||
                    col >= block_count * ColDim
                ? 0.
                : 10. * (1 + (row - row_offset) / RowDim) + col / ColDim;
        if (locMat(row, col) != expected)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Incorrect %dx%d form block or cursor stride at (%zu,%zu)",
                  RowDim, ColDim, row, col);
      }
    }
    MoFEMFunctionReturn(0);
  }
};
 
/** Const matrices must produce read-only, zero-copy views with Gauss strides. */
MoFEMErrorCode checkConstMatrixViews() {
  MoFEMFunctionBegin;
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  using VectorView = GetFTensor1FromMatType<3, -1, DL, const MatrixDouble>;
  using TensorView = GetFTensor2FromMatType<3, 3, -1, DL, const MatrixDouble>;
  using SymmetricView =
      GetFTensor2SymmetricFromMatType<3, -1, DL, const MatrixDouble>;
  static_assert(std::is_same_v<
                VectorView,
                FTensor::Tensor1<FTensor::CursorPtr<const double *, 3>, 3>>);
  static_assert(std::is_same_v<
                TensorView,
                FTensor::Tensor2<FTensor::CursorPtr<const double *, 9>, 3, 3>>);
  static_assert(std::is_same_v<SymmetricView,
                               FTensor::Tensor2_symmetric<
                                   FTensor::CursorPtr<const double *, 6>, 3>>);
 
  constexpr int count = 3;
  MatrixDouble vector_values(count, 3), tensor_values(count, 9),
      symmetric_values(count, 6);
  // Independent scalar storage initialisation and indexing are the oracle.
  for (int gg = 0; gg != count; ++gg) {
    for (int col = 0; col != 3; ++col)
      vector_values(gg, col) = 100. * gg + col + 1;
    for (int col = 0; col != 9; ++col)
      tensor_values(gg, col) = 100. * gg + col + 11;
    for (int col = 0; col != 6; ++col)
      symmetric_values(gg, col) = 100. * gg + col + 21;
  }
  const auto &vectors = vector_values;
  const auto &tensors = tensor_values;
  const auto &symmetric = symmetric_values;
  auto t_vector = MatrixSizeHelper<VectorView, DL>::get(vectors, count)();
  auto t_tensor = MatrixSizeHelper<TensorView, DL>::get(tensors, count)();
  auto t_symmetric =
      MatrixSizeHelper<SymmetricView, DL>::get(symmetric, count)();
  static_assert(!std::is_assignable_v<decltype(t_vector(0)), double>);
  static_assert(!std::is_assignable_v<decltype(t_tensor(0, 0)), double>);
  static_assert(!std::is_assignable_v<decltype(t_symmetric(0, 0)), double>);
 
  auto t_vector_mutable = getFTensor1FromMat<3>(vector_values);
  auto t_tensor_mutable = getFTensor2FromMat<3, 3>(tensor_values);
  auto t_symmetric_mutable = getFTensor2SymmetricFromMat<3>(symmetric_values);
  static_assert(std::is_assignable_v<decltype(t_vector_mutable(0)), double>);
  static_assert(std::is_assignable_v<decltype(t_tensor_mutable(0, 0)), double>);
  static_assert(
      std::is_assignable_v<decltype(t_symmetric_mutable(0, 0)), double>);
  // Updates after creating the const views must remain visible through them.
  t_vector_mutable(0) = -31.;
  t_tensor_mutable(0, 0) = -41.;
  t_symmetric_mutable(0, 0) = -51.;
  for (int gg = 0; gg != count; ++gg) {
    for (int row = 0; row != 3; ++row) {
      if (t_vector(row) != vector_values(gg, row))
        SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                "Incorrect const vector view or Gauss stride");
      for (int col = 0; col != 3; ++col)
        if (t_tensor(row, col) != tensor_values(gg, 3 * row + col))
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Incorrect const tensor view or Gauss stride");
      for (int col = row; col != 3; ++col) {
        const int packed = 3 * row - row * (row - 1) / 2 + col - row;
        if (t_symmetric(row, col) != symmetric_values(gg, packed) ||
            t_symmetric(col, row) != symmetric_values(gg, packed))
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "Incorrect const symmetric view or Gauss stride");
      }
    }
    ++t_vector;
    ++t_tensor;
    ++t_symmetric;
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr,
                          "Check FTensor form matrix block views\n");
  try {
    CHKERR checkConstMatrixViews();
    OpCheckBlocks op;
    CHKERR op.check<2, 2>();
    CHKERR op.check<3, 3>();
    CHKERR op.check<5, 5>();
    CHKERR op.check<6, 6>();
    CHKERR op.check<5, 1>();
    CHKERR op.check<1, 5>();
    CHKERR op.check<5, 3>();
    CHKERR op.check<1, 3>();
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}