test segments distance

/** \file field_evaluator.cpp
* \example field_evaluator.cpp
* \brief test segments distance
* \ingroup field_evaluator
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* License for more details.
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "testing field evaluator\n\n";
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";
// 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,
// set fields order
CHKERR simple_interface->setFieldOrder("FIELD1", 1);
// setup problem
CHKERR simple_interface->setUp();
FieldEvaluatorInterface *field_eval_ptr;
CHKERR m_field.getInterface(field_eval_ptr);
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};
auto dm = simple_interface->getDM();
const MoFEM::Problem *prb_ptr;
* @brief Operator used to check consistency between local coordinates and
* global cooridnates for integrated points and evaluated points
struct MyOp : public VolOp {
MatrixShallowArrayAdaptor<double> evalPoints;
MyOp(std::array<double, 12> &eval_points)
: VolOp("FIELD1", OPROW),
eval_points.size() / 3, 3)) {}
MoFEMErrorCode doWork(int side, EntityType type,
if (type == MBVERTEX) {
std::cout << "Integration pts" << std::endl;
std::cout << getGaussPts() << endl;
std::cout << "Global coordinates " << endl;
std::cout << getCoordsAtGaussPts() << std::endl;
for (int gg = 0; gg != getCoordsAtGaussPts().size1(); ++gg) {
int pt_number = getGaussPts()(3, gg);
std::cout << "gg " << gg << std::endl;
std::cout << "pt " << pt_number << std::endl;
ublas::matrix_row<MatrixDouble> coord_at_gauss_pt(
getCoordsAtGaussPts(), gg);
ublas::matrix_row<MatrixShallowArrayAdaptor<double>> eval_coord(
evalPoints, pt_number);
std::cout << "coord_at_gauss_pt ";
std::cout << coord_at_gauss_pt << std::endl;
std::cout << "eval_coord ";
std::cout << eval_coord << std::endl;
double error = norm_2(coord_at_gauss_pt - eval_coord);
if (error > 1e-12)
"Difference at %d error = %3.4e", pt_number, error);
// 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;
// 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.lock()) {
fe_method->getRuleHook = get_rule;
fe_method->getOpPtrVector().push_back(new MyOp(eval_points));
} else
"Pointer to element does not exists");
// Build tree for particular element
CHKERR field_eval_ptr->buildTree3D(data,
// 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->evalFEAtThePoint3D(
&point[0], dist, prb_ptr->getName(),
simple_interface->getDomainFEName(), data, m_field.get_comm_rank(),
return 0;
Definition: definitions.h:281
Catch errors.
Definition: definitions.h:441
Definition: definitions.h:189
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
Definition: definitions.h:132
Definition: definitions.h:123
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
int main(int argc, char *argv[])
static char help[]
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1129
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:70
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
DataForcesAndSourcesCore::EntData EntData
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Deprecated interface functions.
RuleHookFun getRuleHook
Hook to get rule.
keeps basic data about problem
std::string getName() const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
VolEle::UserDataOperator VolOp