64 {
65
66
68
69 try {
70
71
72 moab::Core moab_core;
73 moab::Interface &moab = moab_core;
74
75
78
79
80 DMType dm_name = "DMMOFEM";
82
83
84 auto core_log = logging::core::get();
85 core_log->add_sink(
89
90
92
93
95
97
98
99
103
104
107
108
110
111
113
115
116
117 auto post_proc = [&](auto dm, auto f_res, auto out_name) {
119 auto post_proc_fe =
120 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(m_field);
121
123
124 auto scal_vec = boost::make_shared<VectorDouble>();
125 auto vec_mat = boost::make_shared<MatrixDouble>();
126
127 post_proc_fe->getOpPtrVector().push_back(
129 post_proc_fe->getOpPtrVector().push_back(
131 f_res));
132
133 post_proc_fe->getOpPtrVector().push_back(
134
136
137 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
138
139 {{"SCALAR", scal_vec}},
140
141 {{"VECTOR", vec_mat}},
142
143 {}, {})
144
145 );
146
148 post_proc_fe);
149 post_proc_fe->writeFile(out_name);
151 };
152
153
154
155
156
157
158 auto TestOpGradGrad = [&]() {
160 MOFEM_LOG(
"OpTester", Sev::verbose) <<
"TestOpGradGrad";
161
162 pip->getOpDomainLhsPipeline().clear();
163 pip->getOpDomainRhsPipeline().clear();
164
165 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
166 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
167
169 "SCALAR", "SCALAR", [](double, double, double) { return 1; }));
171 "VECTOR", "VECTOR", [](double, double, double) { return 1; }));
172
173 auto scl_mat = boost::make_shared<MatrixDouble>();
174 auto vec_mat = boost::make_shared<MatrixDouble>();
175
176 pip->getOpDomainRhsPipeline().push_back(
178 pip->getOpDomainRhsPipeline().push_back(
180 vec_mat));
181
183 "SCALAR", scl_mat, [](double, double, double) { return 1; }));
185 "VECTOR", vec_mat, [](double, double, double) { return 1; }));
186
187 constexpr double eps = 1e-6;
188
189 auto x = opt->setRandomFields(
simple->getDM(),
190 {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
191 auto diff_x = opt->setRandomFields(
192 simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
193
194 auto diff_res = opt->checkCentralFiniteDifference(
195 simple->getDM(),
simple->getDomainFEName(), pip->getDomainRhsFE(),
198
199
200
201 double fnorm;
202 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
203 MOFEM_LOG_C(
"OpTester", Sev::inform,
"TestOpGradGrad %3.4e", fnorm);
204
205 constexpr double err = 1e-9;
206 if (fnorm > err)
208 "Norm of directional derivative too large");
209
211
212
213
215
216 opt->directionalCentralFiniteDifference(
220
221 "out_directional_directive_gradgrad.h5m");
222 }
223
225 };
226
227
228 auto TestOpConvection = [&]() {
230 MOFEM_LOG(
"OpTester", Sev::verbose) <<
"TestOpConvection";
231
232 pip->getOpDomainLhsPipeline().clear();
233 pip->getOpDomainRhsPipeline().clear();
234
235 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
236 pip->setDomainLhsIntegrationRule([](
int,
int,
int o) {
return 2 *
o; });
237
238 auto evaluate_field = [&](auto &pipe) {
239 auto scl_mat = boost::make_shared<MatrixDouble>();
240 auto vec_mat = boost::make_shared<MatrixDouble>();
241 auto dot_vec_vel = boost::make_shared<MatrixDouble>();
242 pipe.push_back(
245 "VECTOR", vec_mat));
247 "VECTOR", dot_vec_vel));
248 return std::make_tuple(scl_mat, vec_mat, dot_vec_vel);
249 };
250
251 auto [rhs_scl_mat, rhs_vec_mat, rhs_dot_vec_vel] =
252 evaluate_field(pip->getOpDomainRhsPipeline());
253 pip->getOpDomainRhsPipeline().push_back(
255 pip->getOpDomainRhsPipeline().push_back(
257 rhs_vec_mat));
258
259 auto [lhs_scl_mat, lhs_vec_mat, lhs_dot_vec_vel] =
260 evaluate_field(pip->getOpDomainLhsPipeline());
261
262
263
264 auto op_convective_term_lhs_du_scalar =
266 op_convective_term_lhs_du_scalar->feScalingFun =
267 [](
const FEMethod *fe_ptr) {
return fe_ptr->ts_a; };
268 pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_scalar);
269 pip->getOpDomainLhsPipeline().push_back(
271
272 auto op_convective_term_lhs_du_vec =
274 op_convective_term_lhs_du_vec->feScalingFun = [](
const FEMethod *fe_ptr) {
275 return fe_ptr->ts_a;
276 };
277 pip->getOpDomainLhsPipeline().push_back(op_convective_term_lhs_du_vec);
278 pip->getOpDomainLhsPipeline().push_back(
280 lhs_dot_vec_vel));
281
282 constexpr double eps = 1e-6;
283
284 auto x = opt->setRandomFields(
simple->getDM(),
285 {{"SCALAR", {-1, 1}}, {"VECTOR", {-1, 1}}});
286 auto x_t = opt->setRandomFields(
287 simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
288 auto diff_x = opt->setRandomFields(
289 simple->getDM(), {{
"SCALAR", {-1, 1}}, {
"VECTOR", {-1, 1}}});
290
291 auto diff_res = opt->checkCentralFiniteDifference(
292 simple->getDM(),
simple->getDomainFEName(), pip->getDomainRhsFE(),
295
296
297
298 double fnorm;
299 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
300 MOFEM_LOG_C(
"OpTester", Sev::inform,
"TestOpConvection %3.4e", fnorm);
301
302 constexpr double err = 1e-9;
303 if (fnorm > err) {
305 "Norm of directional derivative too large");
306 }
307
309
310
311
313
314 opt->directionalCentralFiniteDifference(
316 pip->getDomainRhsFE(), x, x_t,
318
319 "out_directional_directive_convection.h5m");
320 }
321
323 };
324
326 CHKERR TestOpConvection();
327 }
329
330
332
333 return 0;
334}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr int SPACE_DIM
[Define dimension]
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDy< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDy
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpConvectiveTermRhs< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermRhs
FormsIntegrators< DomainEleOp >::Assembly< A >::TriLinearForm< I >::OpConvectiveTermLhsDu< 1, FIELD_DIM, SPACE_DIM > OpConvectiveTermLhsDu
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
structure for User Loop Methods on finite elements
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
PipelineManager interface.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.