v0.14.0
field_evaluator.cpp
Go to the documentation of this file.
1 /** \file field_evaluator.cpp
2  * \example field_evaluator.cpp
3  * \brief test segments distance
4  *
5  * \ingroup field_evaluator
6  */
7 
8 
9 
10 #include <MoFEM.hpp>
11 
12 using namespace MoFEM;
13 
14 static char help[] = "testing field evaluator\n\n";
15 
16 /**
17  * @brief Operator used to check consistency between local coordinates
18  * and global cooridnates for integrated points and evaluated points
19  *
20  */
21 template <typename OP> struct MyOp : public OP {
22 
24 
25  MyOp(std::array<double, 12> &eval_points)
26  : OP("FIELD1", OP::OPROW),
27  evalPoints(
28  getMatrixAdaptor(eval_points.data(), eval_points.size() / 3, 3)) {}
29 
30  MoFEMErrorCode doWork(int side, EntityType type,
33  if (type == MBVERTEX) {
34 
35  MOFEM_LOG("SELF", Sev::inform) << "FE " << OP::getFEEntityHandle();
36 
37  MOFEM_LOG("SELF", Sev::inform) << "Integration pts" << std::endl;
38  MOFEM_LOG("SELF", Sev::inform) << OP::getGaussPts() << endl;
39 
40  MOFEM_LOG("SELF", Sev::inform) << "Global coordinates " << endl;
41  MOFEM_LOG("SELF", Sev::inform) << OP::getCoordsAtGaussPts() << std::endl;
42 
43  for (int gg = 0; gg != OP::getCoordsAtGaussPts().size1(); ++gg) {
44  int pt_number = OP::getGaussPts()(OP::getGaussPts().size1() - 1, gg);
45 
46  MOFEM_LOG("SELF", Sev::inform) << "gg " << gg << std::endl;
47  MOFEM_LOG("SELF", Sev::inform) << "pt " << pt_number << std::endl;
48 
49  ublas::matrix_row<MatrixDouble> coord_at_gauss_pt(
50  OP::getCoordsAtGaussPts(), gg);
51  ublas::matrix_row<MatrixShallowArrayAdaptor<double>> eval_coord(
52  evalPoints, pt_number);
53 
54  MOFEM_LOG("SELF", Sev::inform) << "coord_at_gauss_pt ";
55  MOFEM_LOG("SELF", Sev::inform) << coord_at_gauss_pt << std::endl;
56 
57  MOFEM_LOG("SELF", Sev::inform) << "eval_coord ";
58  MOFEM_LOG("SELF", Sev::inform) << eval_coord << std::endl;
59 
60  double error = norm_2(coord_at_gauss_pt - eval_coord);
61 
62  if (error > 1e-12)
63  SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
64  "Difference at %d error = %3.4e", pt_number, error);
65  }
66  }
68  }
69 };
70 
71 // Lambda function is used to set integration rule to -1, that indicates
72 // that finite element instace will use non-standard integration points.
73 auto get_rule = [](int order_row, int order_col, int order_data) {
74  return -1;
75 };
76 
77 int main(int argc, char *argv[]) {
78 
79  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
80 
81  try {
82 
83  // Read mesh and create moab and mofem data structures
84  moab::Core mb_instance;
85  moab::Interface &moab = mb_instance;
86 
87  MoFEM::Core core(moab, PETSC_COMM_WORLD);
88  MoFEM::Interface &m_field = core;
89 
90  // Register DM Manager
91  DMType dm_name = "DMMOFEM";
92  CHKERR DMRegister_MoFEM(dm_name);
93 
94  // Simple interface
95  Simple *simple_interface;
96  CHKERR m_field.getInterface(simple_interface);
97  {
98 
99  // get options from command line
100  CHKERR simple_interface->getOptions();
101  // load mesh file
102  CHKERR simple_interface->loadFile();
103  // add fields
104  CHKERR simple_interface->addDomainField("FIELD1", H1,
106  // set fields order
107  CHKERR simple_interface->setFieldOrder("FIELD1", 1);
108  // setup problem
109  CHKERR simple_interface->setUp();
110 
111  FieldEvaluatorInterface *field_eval_ptr;
112  CHKERR m_field.getInterface(field_eval_ptr);
113 
114 
115  auto dm = simple_interface->getDM();
116  const MoFEM::Problem *prb_ptr;
117  CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
118 
119  if (simple_interface->getDim() == 3) {
120 
121  std::array<double, 3> point = {0, 0, 0};
122  const double dist = 0.1;
123 
124  std::array<double, 12> eval_points = {0.0, 0.0, 0.0, 0.0, 0.01, 0.0,
125  -0.1, 0.0, 0.0, 0.1, 0.0, 0.0};
126 
129 
130  // Get pointer of FieldEvaluator data. Note finite element and method
131  // set integrating points is destroyed when this pointer is releases
132  auto data = field_eval_ptr->getData<VolEle>();
133 
134  // Set operators and integration rule
135  if (auto fe_method = data->feMethodPtr.lock()) {
136  fe_method->getRuleHook = get_rule;
137  fe_method->getOpPtrVector().push_back(new MyOp<VolOp>(eval_points));
138  } else
139  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
140  "Pointer to element does not exists");
141 
142  // Build tree for particular element
143  CHKERR field_eval_ptr->buildTree3D(data,
144  simple_interface->getDomainFEName());
145  // Set points to set on finite elements
146  data->setEvalPoints(eval_points.data(), eval_points.size() / 3);
147 
148  // Evaluate points on finite elements
149  CHKERR field_eval_ptr->evalFEAtThePoint3D(
150  &point[0], dist, prb_ptr->getName(),
151  simple_interface->getDomainFEName(), data, m_field.get_comm_rank(),
152  m_field.get_comm_rank(), nullptr, MF_EXIST, VERY_NOISY);
153  }
154 
155  if (simple_interface->getDim() == 2) {
156 
157  std::array<double, 3> point = {0, 0, 0};
158  const double dist = 0.1;
159 
160  std::array<double, 12> eval_points = {
161  0.0, 0.0, 0.0,
162  0.0, 0.01, 0.0,
163  -0.1, 0.0, 0.0,
164  0.1, 0.0, 0.0};
165 
168 
169  // Get pointer of FieldEvaluator data. Note finite element and method
170  // set integrating points is destroyed when this pointer is releases
171  auto data = field_eval_ptr->getData<FaceEle>();
172 
173  // Set operators and integration rule
174  if (auto fe_method = data->feMethodPtr.lock()) {
175  fe_method->getRuleHook = get_rule;
176  fe_method->getOpPtrVector().push_back(new MyOp<FaceOp>(eval_points));
177  } else
178  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
179  "Pointer to element does not exists");
180 
181  // Build tree for particular element
182  CHKERR field_eval_ptr->buildTree2D(data,
183  simple_interface->getDomainFEName());
184  // Set points to set on finite elements
185  data->setEvalPoints(eval_points.data(), eval_points.size() / 3);
186 
187  // Evaluate points on finite elements
188  CHKERR field_eval_ptr->evalFEAtThePoint2D(
189  &point[0], dist, prb_ptr->getName(),
190  simple_interface->getDomainFEName(), data, m_field.get_comm_rank(),
191  m_field.get_comm_rank(), nullptr, MF_EXIST, VERY_NOISY);
192  }
193 
194  }
195  }
196  CATCH_ERRORS;
197 
199 
200  return 0;
201 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
help
static char help[]
Definition: field_evaluator.cpp:14
MoFEM::Types::MatrixShallowArrayAdaptor
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:121
MoFEM::ForcesAndSourcesCore::getRuleHook
RuleHookFun getRuleHook
Hook to get rule.
Definition: ForcesAndSourcesCore.hpp:42
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VERY_NOISY
@ VERY_NOISY
Definition: definitions.h:225
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::FieldEvaluatorInterface::buildTree3D
MoFEMErrorCode buildTree3D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
Definition: FieldEvaluator.cpp:57
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::FieldEvaluatorInterface::buildTree2D
MoFEMErrorCode buildTree2D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
Definition: FieldEvaluator.cpp:63
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
MyOp::MyOp
MyOp(std::array< double, 12 > &eval_points)
Definition: field_evaluator.cpp:25
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
convert.type
type
Definition: convert.py:64
MoFEM::FieldEvaluatorInterface::evalFEAtThePoint3D
MoFEMErrorCode evalFEAtThePoint3D(const double *const point, const double distance, const std::string problem, const std::string finite_element, boost::shared_ptr< SetPtsData > data_ptr, int lower_rank, int upper_rank, boost::shared_ptr< CacheTuple > cache_ptr, MoFEMTypes bh=MF_EXIST, VERBOSITY_LEVELS verb=QUIET)
Evaluate field at artbitray position.
Definition: FieldEvaluator.cpp:364
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::getMatrixAdaptor
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:57
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
MyOp::evalPoints
MatrixShallowArrayAdaptor< double > evalPoints
Definition: field_evaluator.cpp:23
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
MoFEM::FieldEvaluatorInterface::getData
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
Definition: FieldEvaluator.hpp:107
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
MyOp::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: field_evaluator.cpp:30
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MyOp
Operator used to check consistency between local coordinates and global cooridnates for integrated po...
Definition: field_evaluator.cpp:21
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MoFEM::Simple::getDim
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:312
get_rule
auto get_rule
Definition: field_evaluator.cpp:73
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::FieldEvaluatorInterface::evalFEAtThePoint2D
MoFEMErrorCode evalFEAtThePoint2D(const double *const point, const double distance, const std::string problem, const std::string finite_element, boost::shared_ptr< SetPtsData > data_ptr, int lower_rank, int upper_rank, boost::shared_ptr< CacheTuple > cache_ptr, MoFEMTypes bh=MF_EXIST, VERBOSITY_LEVELS verb=QUIET)
Evaluate field at artbitray position.
Definition: FieldEvaluator.cpp:373
main
int main(int argc, char *argv[])
Definition: field_evaluator.cpp:77
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
OP
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359