base_tensor_base.cpp

Go to the documentation of this file.

/** @file base_tensor_base.cpp
 * @brief Check mixed scalar/vector/tensor basis bilinear contractions.
 */
 
#define FTENSOR_DEBUG
#include <MoFEM.hpp>
using namespace MoFEM;
 
namespace {
 
using BilinearForm =
    FormsIntegrators<VolumeElementForcesAndSourcesCore::UserDataOperator>::
        Assembly<AssemblyType::PETSC>::BiLinearForm<IntegrationType::GAUSS>;
 
template <int BaseDim, int FieldDim, int SpaceDim, int S = -1>
struct OpCheck
    : BilinearForm::OpBaseTensorBase<BaseDim, FieldDim, SpaceDim, S> {
  using Base = BilinearForm::OpBaseTensorBase<BaseDim, FieldDim, SpaceDim, S>;
  using Base::Base;
 
  MoFEMErrorCode check(ForcesAndSourcesCore &fe,
                       EntitiesFieldData::EntData &row_data,
                       EntitiesFieldData::EntData &col_data,
                       const MatrixDouble &tensor, const double alpha,
                       const bool transpose, const bool only_transpose) {
    MoFEMFunctionBegin;
    constexpr bool tensor_basis = BaseDim == 4 || BaseDim == 9;
    constexpr int col_dim = tensor_basis ? 1 : SpaceDim;
    constexpr int extent = BaseDim == 4 ? 2 : 3;
    CHKERR this->setPtrFE(&fe);
    if (this->assembleTranspose != transpose ||
        this->onlyTranspose != only_transpose || this->sYmm)
      SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
              "Incorrect mixed-form assembly flags");
    this->nbRows = row_data.getIndices().size();
    this->nbCols = col_data.getIndices().size();
    this->nbIntegrationPts = row_data.getN().size1();
    this->nbRowBaseFunctions = row_data.getN().size2();
    this->getGaussPts().resize(4, this->nbIntegrationPts, false);
    this->getGaussPts().clear();
    this->getMeasure() = 1.25;
    for (int gg = 0; gg != this->nbIntegrationPts; ++gg)
      this->getGaussPts()(3, gg) = .17 + .11 * gg;
    this->locMat.resize(this->nbRows, this->nbCols, false);
    std::fill(this->locMat.data().begin(), this->locMat.data().end(), .13);
    CHKERR this->iNtegrate(row_data, col_data);
 
    // Independent scalar oracle checks component order and every local block.
    for (int row = 0; row != this->nbRows; ++row) {
      const int rr = row / FieldDim;
      const int a = row % FieldDim;
      for (int col = 0; col != this->nbCols; ++col) {
        const int cc = col / col_dim;
        const int component = col % col_dim;
        double expected = .13;
        for (int gg = 0; gg != this->nbIntegrationPts; ++gg) {
          const double scale = alpha * this->getMeasure() *
                               this->getGaussPts()(3, gg) *
                               row_data.getN()(gg, rr);
          for (int i = 0; i != SpaceDim; ++i)
            for (int j = 0; j != SpaceDim; ++j)
              if (tensor_basis || i == component) {
                const int basis_component = tensor_basis ? extent * i + j : j;
                expected +=
                    scale *
                    tensor(S == 0 ? 0 : gg, (SpaceDim * i + j) * FieldDim + a) *
                    col_data.getN()(gg, BaseDim * cc + basis_component);
              }
        }
        const double value = this->locMat(row, col);
        if (!std::isfinite(value) ||
            std::abs(value - expected) > 2.e-12 * (1. + std::abs(expected)))
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                  "BASE_DIM=%d FIELD_DIM=%d SPACE_DIM=%d S=%d: block (%d,%d) "
                  "has value %g, expected %g",
                  BaseDim, FieldDim, SpaceDim, S, row, col, value, expected);
      }
    }
    MoFEMFunctionReturn(0);
  }
};
 
template <int BaseDim, int FieldDim, int SpaceDim, int S = -1>
MoFEMErrorCode checkForm(MoFEM::Interface &m_field) {
  MoFEMFunctionBegin;
  constexpr int col_dim = BaseDim == 4 || BaseDim == 9 ? 1 : SpaceDim;
  constexpr double alpha = -.7;
  VolumeElementForcesAndSourcesCore fe(m_field);
  EntitiesFieldData::EntData row_data, col_data;
  row_data.getBase() = col_data.getBase() = AINSWORTH_LEGENDRE_BASE;
  row_data.getIndices().resize(2 * FieldDim, false);
  col_data.getIndices().resize(3 * col_dim, false);
  auto tensor = boost::make_shared<MatrixDouble>();
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  using Tensor = GetFTensor3FromMatType<SpaceDim, SpaceDim, FieldDim, S, DL>;
 
  for (const int points : {3, 4}) {
    MatrixSizeHelper<Tensor, DL>::size(*tensor, S == 0 ? 1 : points);
    for (int gg = 0; gg != tensor->size1(); ++gg)
      for (int i = 0; i != SpaceDim; ++i)
        for (int j = 0; j != SpaceDim; ++j)
          for (int a = 0; a != FieldDim; ++a)
            (*tensor)(gg, (SpaceDim * i + j) * FieldDim + a) =
                .1 * (1 + a) + .021 * (1 + i) * (2 + j) +
                .037 * gg * (2 + i - a);
    row_data.getN().resize(points, 3, false);
    col_data.getN().resize(points, 4 * BaseDim, false);
    for (int gg = 0; gg != points; ++gg) {
      for (int rr = 0; rr != 3; ++rr)
        row_data.getN()(gg, rr) = rr == 2 ? 37. : .2 * (gg + 1) - .13 * rr;
      for (int cc = 0; cc != 4; ++cc)
        for (int c = 0; c != BaseDim; ++c)
          col_data.getN()(gg, BaseDim * cc + c) =
              cc == 3 ? 53. : .09 * (gg + 1) * (c + 2) - .17 * (cc + 1);
    }
    for (const int mode : {0, 1, 2}) {
      const bool transpose = mode != 0;
      const bool only_transpose = mode == 2;
      OpCheck<BaseDim, FieldDim, SpaceDim, S> op(
          "row", "col", tensor, [alpha]() { return alpha; }, transpose,
          only_transpose);
      CHKERR op.check(fe, row_data, col_data, *tensor, alpha, transpose,
                      only_transpose);
    }
  }
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr,
                          "Check mixed tensor bilinear integration\n");
  try {
    moab::Core moab;
    MoFEM::Core core(moab);
    CHKERR checkForm<2, 1, 2>(core);
    CHKERR checkForm<2, 3, 2>(core);
    CHKERR checkForm<3, 3, 2>(core);
    CHKERR checkForm<4, 1, 2>(core);
    CHKERR checkForm<4, 3, 2>(core);
    CHKERR checkForm<9, 3, 2>(core);
    CHKERR checkForm<3, 1, 3>(core);
    CHKERR checkForm<3, 5, 3>(core);
    CHKERR checkForm<9, 1, 3>(core);
    CHKERR checkForm<9, 5, 3>(core);
    CHKERR checkForm<3, 5, 3, 0>(core);
    CHKERR checkForm<9, 3, 2, 0>(core);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}