symmetric_deviator_field_values.cpp

Go to the documentation of this file.

/** @file symmetric_deviator_field_values.cpp
 * @brief Check planar and spatial deviator interpolation and vector sources.
 */
 
#define FTENSOR_DEBUG
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <array>
#include <cmath>
 
namespace {
 
template <int Dim>
struct TestOp : OpCalculateTensor2SymmetricDeviatorFieldValues<Dim> {
  using Base = OpCalculateTensor2SymmetricDeviatorFieldValues<Dim>;
  using Base::Base;
 
  MoFEMErrorCode attach(ForcesAndSourcesCore &fe, const int points) {
    MoFEMFunctionBegin;
    CHKERR this->setPtrFE(&fe);
    this->getGaussPts().resize(4, points, false);
    MoFEMFunctionReturn(0);
  }
};
 
MoFEMErrorCode checkNear(const double value, const double expected) {
  MoFEMFunctionBegin;
  if (!std::isfinite(value) || std::abs(value - expected) > 2.e-13)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Deviator field value %g differs from %g", value, expected);
  MoFEMFunctionReturn(0);
}
 
template <int Dim>
MoFEMErrorCode checkValues(const MatrixDouble &tensor,
                           const MatrixDouble &coordinates,
                           EntitiesFieldData::EntData &data,
                           const VectorDouble &coefficients,
                           const double scale) {
  MoFEMFunctionBegin;
  constexpr int size = Dim == 2 ? 3 : 5;
  const int points = data.getN().size1();
  if (tensor.size1() != points || tensor.size2() != 6 ||
      coordinates.size1() != points || coordinates.size2() != size)
    SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
            "Incorrect deviator output layout");
  for (int gg = 0; gg != points; ++gg) {
    std::array<double, 5> d{};
    // Independent scalar oracle for interpolation and the documented basis.
    for (int c = 0; c != size; ++c) {
      for (int b = 0; b != 2; ++b)
        d[c] += scale * coefficients[b * size + c] * data.getN()(gg, b);
      CHKERR checkNear(coordinates(gg, c), d[c]);
    }
    const std::array<double, 6> expected{
        d[0] / std::sqrt(2.) + d[1] / std::sqrt(6.),
        d[2] / std::sqrt(2.), d[3] / std::sqrt(2.),
        -d[0] / std::sqrt(2.) + d[1] / std::sqrt(6.),
        d[4] / std::sqrt(2.), -2. * d[1] / std::sqrt(6.)};
    for (int c = 0; c != 6; ++c)
      CHKERR checkNear(tensor(gg, c), expected[c]);
    CHKERR checkNear(tensor(gg, 0) + tensor(gg, 3) + tensor(gg, 5), 0.);
    const double norm = tensor(gg, 0) * tensor(gg, 0) +
                        tensor(gg, 3) * tensor(gg, 3) +
                        tensor(gg, 5) * tensor(gg, 5) +
                        2. * (tensor(gg, 1) * tensor(gg, 1) +
                              tensor(gg, 2) * tensor(gg, 2) +
                              tensor(gg, 4) * tensor(gg, 4));
    double coordinate_norm = 0;
    for (const auto value : d)
      coordinate_norm += value * value;
    CHKERR checkNear(norm, coordinate_norm);
  }
  MoFEMFunctionReturn(0);
}
 
template <int Dim> MoFEMErrorCode checkField(MoFEM::Interface &m_field) {
  MoFEMFunctionBegin;
  constexpr int size = Dim == 2 ? 3 : 5;
  VolumeElementForcesAndSourcesCore fe(m_field);
  EntitiesFieldData::EntData data;
  data.getBase() = AINSWORTH_LEGENDRE_BASE;
  data.getFieldData().resize(2 * size, false);
  for (int c = 0; c != 2 * size; ++c)
    data.getFieldData()[c] = .07 * (c + 1);
  const VectorDouble original = data.getFieldData();
  auto tensor = boost::make_shared<MatrixDouble>();
  auto coordinates = boost::make_shared<MatrixDouble>();
  TestOp<Dim> op("D", tensor, SmartPetscObj<Vec>(), MBVERTEX, 0, coordinates);
 
  for (const int points : {3, 4}) {
    CHKERR op.attach(fe, points);
    data.getN().resize(points, 3, false);
    for (int gg = 0; gg != points; ++gg) {
      data.getN()(gg, 0) = .2 * (gg + 1);
      data.getN()(gg, 1) = -.1 * (gg + 2);
      data.getN()(gg, 2) = 17.; // Inactive basis tail must be skipped.
    }
    CHKERR op.doWork(0, MBVERTEX, data);
    CHKERR checkValues<Dim>(*tensor, *coordinates, data, original, 1.);
    CHKERR op.doWork(1, MBVERTEX, data);
    CHKERR checkValues<Dim>(*tensor, *coordinates, data, original, 2.);
    CHKERR op.doWork(0, MBVERTEX, data);
    CHKERR checkValues<Dim>(*tensor, *coordinates, data, original, 1.);
  }
 
  Vec raw_vector;
  CHKERR VecCreateSeq(PETSC_COMM_SELF, 2 * size, &raw_vector);
  SmartPetscObj<Vec> vector(raw_vector);
  PetscScalar *array;
  CHKERR VecGetArray(vector, &array);
  data.getLocalIndices().resize(2 * size, false);
  VectorDouble from_vector(2 * size);
  for (int c = 0; c != 2 * size; ++c)
    array[c] = -.03 * (c + 1);
  for (int c = 0; c != 2 * size; ++c) {
    const int index = c == 1 ? -1 : 2 * size - c - 1;
    data.getLocalIndices()[c] = index;
    from_vector[c] = index == -1 ? 0. : array[index];
  }
  CHKERR VecRestoreArray(vector, &array);
  TestOp<Dim> vector_op("D", tensor, vector, MBVERTEX, 0, coordinates);
  CHKERR vector_op.attach(fe, data.getN().size1());
  CHKERR vector_op.doWork(0, MBVERTEX, data);
  CHKERR checkValues<Dim>(*tensor, *coordinates, data, from_vector, 1.);
  for (int c = 0; c != 2 * size; ++c)
    CHKERR checkNear(data.getFieldData()[c], original[c]);
 
  data.getFieldData().resize(0, false);
  CHKERR vector_op.doWork(0, MBVERTEX, data);
  CHKERR checkValues<Dim>(*tensor, *coordinates, data, original, 0.);
  CHKERR PetscPrintf(PETSC_COMM_SELF, "%dD deviator field evaluation passed\n", Dim);
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, nullptr,
                          "Check symmetric deviator field interpolation\n");
  try {
    moab::Core moab;
    MoFEM::Core core(moab);
    CHKERR checkField<2>(core);
    CHKERR checkField<3>(core);
  }
  CATCH_ERRORS;
  MoFEM::Core::Finalize();
  return 0;
}