v0.13.1
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
12using namespace MoFEM;
13
14static 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 */
21template <typename OP> struct MyOp : public OP {
22
24
25 MyOp(std::array<double, 12> &eval_points)
26 : OP("FIELD1", OP::OPROW),
28 getMatrixAdaptor(eval_points.data(), eval_points.size() / 3, 3)) {}
29
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.
73auto get_rule = [](int order_row, int order_col, int order_data) {
74 return -1;
75};
76
77int 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 }
197
199
200 return 0;
201}
@ VERY_NOISY
Definition: definitions.h:212
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_EXIST
Definition: definitions.h:100
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
int main(int argc, char *argv[])
static char help[]
auto get_rule
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
ublas::matrix< double, ublas::row_major, ublas::shallow_array_adaptor< double > > MatrixShallowArrayAdaptor
Definition: Types.hpp:121
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
CoreTmp< 0 > Core
Definition: Core.hpp:1086
auto getMatrixAdaptor(T1 ptr, const size_t n, const size_t m)
Get Matrix adaptor.
Definition: Templates.hpp:57
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Field evaluator interface.
MoFEMErrorCode buildTree3D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
MoFEMErrorCode buildTree2D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
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.
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.
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.
@ OPROW
operator doWork function is executed on FE rows
RuleHookFun getRuleHook
Hook to get rule.
keeps basic data about problem
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:264
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:313
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator used to check consistency between local coordinates and global cooridnates for integrated po...
MatrixShallowArrayAdaptor< double > evalPoints
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MyOp(std::array< double, 12 > &eval_points)