field_evaluator.cpp File Reference

#include <MoFEM.hpp>

Go to the source code of this file.

Classes
struct	MyOp< OP >
	Operator used to check consistency between local coordinates and global cooridnates for integrated points and evaluated points. More...

Functions
int	main (int argc, char *argv[])

Variables
static char	help [] = "testing field evaluator\n\n"
auto	get_rule

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 77 of file field_evaluator.cpp.

                                 {
 
  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;
}

Variable Documentation

get_rule

auto get_rule

Initial value:

= [](int order_row, int order_col, int order_data) {
  return -1;
}

Examples

mofem/atom_tests/field_evaluator.cpp, and mofem/tutorials/max-1_lorentz_force//lorentz_force.cpp.

Definition at line 73 of file field_evaluator.cpp.

                                                                 {
  return -1;
};

help

char help[] = "testing field evaluator\n\n"

static

Definition at line 14 of file field_evaluator.cpp.