field_evaluator.cpp

Go to the documentation of this file.

/** \file field_evaluator.cpp
 * \example mofem/atom_tests/field_evaluator.cpp
 * \brief test segments distance
 *
 * \ingroup field_evaluator
 */
 
 
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
static char help[] = "testing field evaluator\n\n";
 
/**
 * @brief Operator used to check consistency between local coordinates
 * and global cooridnates for integrated points and evaluated points
 *
 */
template <typename OP> struct MyOp : public OP {
 
  MatrixShallowArrayAdaptor<double> evalPoints;
 
  MyOp(std::array<double, 12> &eval_points)
      : OP("FIELD1", OP::OPROW),
        evalPoints(
            getMatrixAdaptor(eval_points.data(), eval_points.size() / 3, 3)) {}
 
  MoFEMErrorCode doWork(int side, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBegin;
    if (type == MBVERTEX) {
 
      MOFEM_LOG("SELF", Sev::inform) << "FE " << OP::getFEEntityHandle();
 
      MOFEM_LOG("SELF", Sev::inform) << "Integration pts" << std::endl;
      MOFEM_LOG("SELF", Sev::inform) << OP::getGaussPts() << endl;
 
      MOFEM_LOG("SELF", Sev::inform) << "Global coordinates " << endl;
      MOFEM_LOG("SELF", Sev::inform) << OP::getCoordsAtGaussPts() << std::endl;
 
      for (int gg = 0; gg != OP::getCoordsAtGaussPts().size1(); ++gg) {
        int pt_number = OP::getGaussPts()(OP::getGaussPts().size1() - 1, gg);
 
        MOFEM_LOG("SELF", Sev::inform) << "gg " << gg << std::endl;
        MOFEM_LOG("SELF", Sev::inform) << "pt " << pt_number << std::endl;
 
        ublas::matrix_row<MatrixDouble> coord_at_gauss_pt(
            OP::getCoordsAtGaussPts(), gg);
        ublas::matrix_row<MatrixShallowArrayAdaptor<double>> eval_coord(
            evalPoints, pt_number);
 
        MOFEM_LOG("SELF", Sev::inform) << "coord_at_gauss_pt ";
        MOFEM_LOG("SELF", Sev::inform) << coord_at_gauss_pt << std::endl;
 
        MOFEM_LOG("SELF", Sev::inform) << "eval_coord ";
        MOFEM_LOG("SELF", Sev::inform) << eval_coord << std::endl;
 
        double error = norm_2(coord_at_gauss_pt - eval_coord);
 
        if (error > 1e-12)
          SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
                   "Difference at %d error = %3.4e", pt_number, error);
      }
    }
    MoFEMFunctionReturn(0);
  }
};
 
// Lambda function is used to set integration rule to -1, that indicates
// that finite element instace will use non-standard integration points.
auto get_rule = [](int order_row, int order_col, int order_data) {
  return -1;
};
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    // Read mesh and create moab and mofem data structures
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    MoFEM::Core core(moab, PETSC_COMM_WORLD);
    MoFEM::Interface &m_field = core;
 
    // Register DM Manager
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Simple interface
    Simple *simple_interface;
    CHKERR m_field.getInterface(simple_interface);
    {
 
      // get options from command line
      CHKERR simple_interface->getOptions();
      // load mesh file
      CHKERR simple_interface->loadFile();
      // add fields
      CHKERR simple_interface->addDomainField("FIELD1", H1,
                                              AINSWORTH_LEGENDRE_BASE, 3);
      // set fields order
      CHKERR simple_interface->setFieldOrder("FIELD1", 1);
      // setup problem
      CHKERR simple_interface->setUp();
 
      // Exercise the TestingInterface adapter with a known constant field.
      CHKERR m_field.getInterface<FieldBlas>()->setField(1.25, "FIELD1");
      auto test_point = [&]<int DIM>() {
        MoFEMFunctionBegin;
        auto testing = m_field.getInterface<TestingInterface>();
        auto values = boost::make_shared<MatrixDouble>();
        auto add_ops = [&](auto &pipeline) {
          MoFEMFunctionBegin;
          pipeline.push_back(new OpCalculateVectorFieldValues<3>("FIELD1", values));
          MoFEMFunctionReturn(0);
        };
        CHKERR testing->addFieldEvaluator<DIM>(
            simple_interface->getProblemName(), simple_interface->getDomainFEName(),
            {0., 0.01, 0.}, add_ops);
        using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
        auto t_value = MatrixSizeHelper<GetFTensor1FromMatType<3, -1, DL>,
                                       DL>::get(*values, 1)();
        FTENSOR_INDEX(3, i);
        FTensor::Tensor1<double, 3> t_expected{1.25, 1.25, 1.25};
        FTensor::Tensor1<double, 3> t_error;
        t_error(i) = t_value(i) - t_expected(i);
        CHKERR testing->checkError<TestType::POINTWISE>(
            "constant field evaluation", std::sqrt(t_error(i) * t_error(i)),
            0., {1e-12, 0.});
 
        CHKERR PetscPushErrorHandler(PetscReturnErrorHandler, nullptr);
        const auto missing_error = testing->addFieldEvaluator<DIM>(
            simple_interface->getProblemName(), simple_interface->getDomainFEName(),
            {1e6, 1e6, 1e6}, add_ops);
        CHKERR PetscPopErrorHandler();
        if (missing_error != MOFEM_ATOM_TEST_INVALID)
          SETERRQ(m_field.get_comm(), MOFEM_ATOM_TEST_INVALID,
                  "Missing evaluation point did not fail the test");
 
        CHKERR PetscPushErrorHandler(PetscReturnErrorHandler, nullptr);
        const auto hook_error = testing->addFieldEvaluator<DIM>(
            simple_interface->getProblemName(), simple_interface->getDomainFEName(),
            {0., 0.01, 0.}, [](auto &) -> MoFEMErrorCode {
              return MOFEM_ATOM_TEST_INVALID;
            });
        CHKERR PetscPopErrorHandler();
        if (hook_error != MOFEM_ATOM_TEST_INVALID)
          SETERRQ(m_field.get_comm(), MOFEM_ATOM_TEST_INVALID,
                  "Pipeline hook failure was not propagated");
        MoFEMFunctionReturn(0);
      };
      if (simple_interface->getDim() == 3)
        CHKERR test_point.template operator()<3>();
      else if (simple_interface->getDim() == 2)
        CHKERR test_point.template operator()<2>();
 
      FieldEvaluatorInterface *field_eval_ptr;
      CHKERR m_field.getInterface(field_eval_ptr);
 
 
      auto dm = simple_interface->getDM();
      const MoFEM::Problem *prb_ptr;
      CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
 
      if (simple_interface->getDim() == 3) {
 
        std::array<double, 3> point = {0, 0, 0};
        const double dist = 0.1;
 
        std::array<double, 12> eval_points = {0.0,  0.0, 0.0, 0.0, 0.01, 0.0,
                                              -0.1, 0.0, 0.0, 0.1, 0.0,  0.0};
 
        using VolEle = VolumeElementForcesAndSourcesCore;
        using VolOp = VolEle::UserDataOperator;
 
        // Get pointer of FieldEvaluator data. Note finite element and method
        // set integrating points is destroyed when this pointer is releases
        auto data = field_eval_ptr->getData<VolEle>();
 
        // Set operators and integration rule
        if (auto fe_method = data->feMethodPtr) {
          fe_method->getRuleHook = get_rule;
          fe_method->getOpPtrVector().push_back(new MyOp<VolOp>(eval_points));
        } else
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Pointer to element does not exists");
 
        // Build tree for particular element
        CHKERR field_eval_ptr->buildTree<3>(data,
                                           simple_interface->getDomainFEName());
        // Set points to set on finite elements
        data->setEvalPoints(eval_points.data(), eval_points.size() / 3);
 
        // Evaluate points on finite elements
        CHKERR field_eval_ptr->evalFEAtThePoint<3>(
            &point[0], dist, prb_ptr->getName(),
            simple_interface->getDomainFEName(), data, m_field.get_comm_rank(),
            m_field.get_comm_rank(), nullptr, MF_EXIST, VERY_NOISY);
      }
 
      if (simple_interface->getDim() == 2) {
 
        std::array<double, 3> point = {0, 0, 0};
        const double dist = 0.1;
 
        std::array<double, 12> eval_points = {
          0.0,  0.0, 0.0, 
          0.0, 0.01, 0.0,
          -0.1, 0.0, 0.0, 
          0.1, 0.0,  0.0};
 
        using FaceEle = FaceElementForcesAndSourcesCore;
        using FaceOp = FaceEle::UserDataOperator;
 
        // Get pointer of FieldEvaluator data. Note finite element and method
        // set integrating points is destroyed when this pointer is releases
        auto data = field_eval_ptr->getData<FaceEle>();
 
        // Set operators and integration rule
        if (auto fe_method = data->feMethodPtr) {
          fe_method->getRuleHook = get_rule;
          fe_method->getOpPtrVector().push_back(new MyOp<FaceOp>(eval_points));
        } else
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Pointer to element does not exists");
 
        // Build tree for particular element
        CHKERR field_eval_ptr->buildTree<2>(data,
                                           simple_interface->getDomainFEName());
        // Set points to set on finite elements
        data->setEvalPoints(eval_points.data(), eval_points.size() / 3);
 
        // Evaluate points on finite elements
        CHKERR field_eval_ptr->evalFEAtThePoint<2>(
            &point[0], dist, prb_ptr->getName(),
            simple_interface->getDomainFEName(), data, m_field.get_comm_rank(),
            m_field.get_comm_rank(), nullptr, MF_EXIST, VERY_NOISY);
      }
 
    }
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}