v0.14.0
Searching...
No Matches
scalar_check_approximation.cpp
Go to the documentation of this file.
1/**
2 * \file scalar_check_approximation.cpp
3 * \example scalar_check_approximation.cpp
4 *
5 * Approximate polynomial in 2D and check derivatives
6 *
7 */
8
9#include <MoFEM.hpp>
10
11using namespace MoFEM;
12
13static char help[] = "...\n\n";
14
15static int approx_order = 4;
16
17template <int DIM> struct ApproxFunctionsImpl {};
18
19template <int DIM> struct ElementsAndOps {};
20
21constexpr int SPACE_DIM =
22 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
23
28
31
32template <> struct ApproxFunctionsImpl<2> {
33 static double fUn(const double x, const double y, double z) {
34 double r = 1;
35 for (int o = 1; o <= approx_order; ++o) {
36 for (int i = 0; i <= o; ++i) {
37 int j = o - i;
38 if (j >= 0)
39 r += pow(x, i) * pow(y, j);
40 }
41 }
42 return r;
43 }
44
45 static FTensor::Tensor1<double, 2> diffFun(const double x, const double y,
46 double z) {
48 for (int o = 1; o <= approx_order; ++o) {
49 for (int i = 0; i <= o; ++i) {
50 int j = o - i;
51 if (j >= 0) {
52 r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) : 0;
53 r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) : 0;
54 }
55 }
56 }
57 return r;
58 }
59};
60
61template <> struct ApproxFunctionsImpl<3> {
62 static double fUn(const double x, const double y, double z) {
63 double r = 1;
64 for (int o = 1; o <= approx_order; ++o) {
65 for (int i = 0; i <= o; ++i) {
66 for (int j = 0; j <= o - i; j++) {
67 int k = o - i - j;
68 if (k >= 0) {
69 r += pow(x, i) * pow(y, j) * pow(z, k);
70 }
71 }
72 }
73 }
74 return r;
75 }
76
77 static FTensor::Tensor1<double, 3> diffFun(const double x, const double y,
78 double z) {
79 FTensor::Tensor1<double, 3> r{0., 0., 0.};
80 for (int o = 1; o <= approx_order; ++o) {
81 for (int i = 0; i <= o; ++i) {
82 for (int j = 0; j <= o - i; j++) {
83 int k = o - i - j;
84 if (k >= 0) {
85 r(0) += i > 0 ? i * pow(x, i - 1) * pow(y, j) * pow(z, k) : 0;
86 r(1) += j > 0 ? j * pow(x, i) * pow(y, j - 1) * pow(z, k) : 0;
87 r(2) += k > 0 ? k * pow(x, i) * pow(y, j) * pow(z, k - 1) : 0;
88 }
89 }
90 }
91 }
92 return r;
93 }
94};
95
97
98struct OpValsDiffVals : public DomainEleOp {
102 OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
103 : DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
105
107
111 const int nb_gauss_pts = getGaussPts().size2();
112 if (type == MBVERTEX) {
113 vAls.resize(nb_gauss_pts, false);
114 diffVals.resize(SPACE_DIM, nb_gauss_pts, false);
115 vAls.clear();
116 diffVals.clear();
117 }
118
119 const int nb_dofs = data.getIndices().size();
120 if (nb_dofs) {
121
122 MOFEM_LOG("AT", Sev::noisy)
123 << "Type " << moab::CN::EntityTypeName(type) << " side " << side;
124 MOFEM_LOG("AT", Sev::noisy) << data.getN();
125 MOFEM_LOG("AT", Sev::noisy) << data.getDiffN();
126
127 auto t_vals = getFTensor0FromVec(vAls);
128 auto t_base_fun = data.getFTensor0N();
129 for (int gg = 0; gg != nb_gauss_pts; gg++) {
130 auto t_data = data.getFTensor0FieldData();
131 for (int bb = 0; bb != nb_dofs; bb++) {
132 t_vals += t_base_fun * t_data;
133 ++t_base_fun;
134 ++t_data;
135 }
136 const double v = t_vals;
137 if (!std::isnormal(v))
138 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
139 ++t_vals;
140 }
141
143 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
144 auto t_diff_base_fun = data.getFTensor1DiffN<SPACE_DIM>();
145 for (int gg = 0; gg != nb_gauss_pts; gg++) {
146 auto t_data = data.getFTensor0FieldData();
147 for (int bb = 0; bb != nb_dofs; bb++) {
148 t_diff_vals(i) += t_diff_base_fun(i) * t_data;
149 ++t_diff_base_fun;
150 ++t_data;
151 }
152 for (int d = 0; d != SPACE_DIM; ++d)
153 if (!std::isnormal(t_diff_vals(d)))
154 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Not a number");
155 ++t_diff_vals;
156 }
157 }
158 }
160 }
161};
162
163struct OpCheckValsDiffVals : public DomainEleOp {
166 boost::shared_ptr<VectorDouble> ptrVals;
167 boost::shared_ptr<MatrixDouble> ptrDiffVals;
169
171 boost::shared_ptr<VectorDouble> &ptr_vals,
172 boost::shared_ptr<MatrixDouble> &ptr_diff_vals,
174 : DomainEleOp("FIELD1", OPROW), vAls(vals), diffVals(diff_vals),
175 ptrVals(ptr_vals), ptrDiffVals(ptr_diff_vals),
177 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
178 }
179
181
185 const double eps = 1e-6;
186 const int nb_gauss_pts = data.getN().size1();
187
188 auto t_vals = getFTensor0FromVec(vAls);
189
190 auto t_ptr_vals = getFTensor0FromVec(*ptrVals);
191
192 for (int gg = 0; gg != nb_gauss_pts; gg++) {
193
194 double err_val;
195
196 // Check user data operators
197 err_val = std::abs(t_vals - t_ptr_vals);
198 MOFEM_LOG("AT", Sev::noisy) << "Val op error " << err_val;
199
200 if (err_val > eps)
201 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
202 "Wrong value from operator %4.3e", err_val);
203
204 const double x = getCoordsAtGaussPts()(gg, 0);
205 const double y = getCoordsAtGaussPts()(gg, 1);
206 const double z = getCoordsAtGaussPts()(gg, 2);
207
208 // Check approximation
209 const double delta_val = t_vals - ApproxFunctions::fUn(x, y, z);
210
211 err_val = std::fabs(delta_val * delta_val);
212 MOFEM_LOG("AT", Sev::verbose) << err_val << " : " << t_vals;
213 if (err_val > eps)
214 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Wrong value %4.3e",
215 err_val);
216
217 ++t_vals;
218 ++t_ptr_vals;
219 }
220
222
223 auto t_diff_vals = getFTensor1FromMat<SPACE_DIM>(diffVals);
224 auto t_ptr_diff_vals = getFTensor1FromMat<SPACE_DIM>(*ptrDiffVals);
225
226 for (int gg = 0; gg != nb_gauss_pts; gg++) {
227
229 double err_diff_val;
230
231 t_delta_diff_val(i) = t_diff_vals(i) - t_ptr_diff_vals(i);
232 err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
233 MOFEM_LOG("AT", Sev::noisy) << "Diff val op error " << err_diff_val;
234
235 if (err_diff_val > eps)
236 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
237 "Wrong derivatives from operator %4.3e", err_diff_val);
238
239 const double x = getCoordsAtGaussPts()(gg, 0);
240 const double y = getCoordsAtGaussPts()(gg, 1);
241 const double z = getCoordsAtGaussPts()(gg, 2);
242
243 // Check approximation
244 auto t_diff_anal = ApproxFunctions::diffFun(x, y, z);
245 t_delta_diff_val(i) = t_diff_vals(i) - t_diff_anal(i);
246
247 err_diff_val = sqrt(t_delta_diff_val(i) * t_delta_diff_val(i));
248 if (SPACE_DIM == 3)
249 MOFEM_LOG("AT", Sev::noisy)
250 << "Diff val " << err_diff_val << " : "
251 << sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
252 << t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
253 << t_diff_vals(1) << " (" << t_diff_anal(1) << ") "
254 << t_diff_vals(2) << " (" << t_diff_anal(2) << ")";
255 else
256 MOFEM_LOG("AT", Sev::noisy)
257 << "Diff val " << err_diff_val << " : "
258 << sqrt(t_diff_vals(i) * t_diff_vals(i)) << " : "
259 << t_diff_vals(0) << " (" << t_diff_anal(0) << ") "
260 << t_diff_vals(1) << " (" << t_diff_anal(1) << ")";
261
262 MOFEM_LOG("AT", Sev::verbose)
263 << getCoords()(3 * 1 + 0) - getCoords()(3 * 0 + 0);
264 MOFEM_LOG("AT", Sev::verbose)
265 << getCoords()(3 * 1 + 1) - getCoords()(3 * 0 + 1);
266 MOFEM_LOG("AT", Sev::verbose)
267 << getCoords()(3 * 1 + 2) - getCoords()(3 * 0 + 2);
268
269 MOFEM_LOG("AT", Sev::verbose) << "Diff val error " << err_diff_val;
270 if (err_diff_val > eps)
271 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
272 "Wrong derivative of value %4.3e %4.3e", err_diff_val,
273 t_diff_anal.l2());
274
275 ++t_diff_vals;
276 ++t_ptr_diff_vals;
277 }
278 }
279
281 }
282};
283
284int main(int argc, char *argv[]) {
285
286 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
287
288 try {
289
290 DMType dm_name = "DMMOFEM";
291 CHKERR DMRegister_MoFEM(dm_name);
292
293 moab::Core mb_instance;
294 moab::Interface &moab = mb_instance;
295
296 // Add logging channel for example
297 auto core_log = logging::core::get();
300 LogManager::setLog("AT");
301 MOFEM_LOG_TAG("AT", "atom_test");
302
303 // Create MoFEM instance
304 MoFEM::Core core(moab);
305 MoFEM::Interface &m_field = core;
306
307 Simple *simple = m_field.getInterface<Simple>();
308 PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
309 CHKERR simple->getOptions();
310
312
314
315 // Declare elements
316 enum bases {
317 AINSWORTH,
318 AINSWORTH_LOBATTO,
319 DEMKOWICZ,
320 BERNSTEIN,
321 LASBASETOP
322 };
323 const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
324 "bernstein"};
325 PetscBool flg;
326 PetscInt choice_base_value = AINSWORTH;
327 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
328 LASBASETOP, &choice_base_value, &flg);
329
330 if (flg != PETSC_TRUE)
331 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "base not set");
333 if (choice_base_value == AINSWORTH)
335 if (choice_base_value == AINSWORTH_LOBATTO)
337 else if (choice_base_value == DEMKOWICZ)
339 else if (choice_base_value == BERNSTEIN)
341
342 enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
343 const char *list_spaces[] = {"h1", "l2"};
344 PetscInt choice_space_value = H1SPACE;
345 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
346 LASBASETSPACE, &choice_space_value, &flg);
347 if (flg != PETSC_TRUE)
348 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "space not set");
349 FieldSpace space = H1;
350 if (choice_space_value == H1SPACE)
351 space = H1;
352 else if (choice_space_value == L2SPACE)
353 space = L2;
354
355 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
356 PETSC_NULL);
357
358 CHKERR simple->addDomainField("FIELD1", space, base, 1);
359 CHKERR simple->setFieldOrder("FIELD1", approx_order);
360
361 Range edges, faces;
362 CHKERR moab.get_entities_by_dimension(0, 1, edges);
363 CHKERR moab.get_entities_by_dimension(0, 2, faces);
364
365 if (choice_base_value != BERNSTEIN) {
366 Range rise_order;
367
368 int i = 0;
369 for (auto e : edges) {
370 if (!(i % 2)) {
371 rise_order.insert(e);
372 }
373 ++i;
374 }
375
376 for (auto f : faces) {
377 if (!(i % 3)) {
378 rise_order.insert(f);
379 }
380 ++i;
381 }
382
383 Range rise_twice;
384 for (auto e : rise_order) {
385 if (!(i % 2)) {
386 rise_twice.insert(e);
387 }
388 ++i;
389 }
390
391 CHKERR simple->setFieldOrder("FIELD1", approx_order + 1, &rise_order);
392
393 CHKERR simple->setFieldOrder("FIELD1", approx_order + 2, &rise_twice);
394 }
395
396 CHKERR simple->defineFiniteElements();
397
398 auto volume_adj = [](moab::Interface &moab, const Field &field,
399 const EntFiniteElement &fe,
402 EntityHandle fe_ent = fe.getEnt();
403 switch (field.getSpace()) {
404 case H1:
405 CHKERR moab.get_connectivity(&fe_ent, 1, adjacency, true);
406 case HCURL:
408 moab::Interface::UNION);
409 case HDIV:
410 case L2:
412 moab::Interface::UNION);
414 // build side table
415 for (auto ent : adjacency)
416 fe.getSideNumberPtr(ent);
417 break;
418 case NOFIELD: {
420 false);
421 for (auto ent : adjacency) {
422 const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
423 .insert(
424 boost::shared_ptr<SideNumber>(new SideNumber(ent, -1, 0, 0)));
425 }
426 } break;
427 default:
428 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
429 "this field is not implemented for TRI finite element");
430 }
431
433 };
434
439
440 CHKERR simple->defineProblem(PETSC_TRUE);
441 CHKERR simple->buildFields();
442 CHKERR simple->buildFiniteElements();
443 CHKERR simple->buildProblem();
444
445 auto dm = simple->getDM();
446
447 VectorDouble vals;
448 auto jac_ptr = boost::make_shared<MatrixDouble>();
449 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
450 auto det_ptr = boost::make_shared<VectorDouble>();
451 MatrixDouble diff_vals;
452
453 auto assemble_matrices_and_vectors = [&]() {
455
458
460 pipeline_mng->getOpDomainLhsPipeline(), {NOSPACE});
462 pipeline_mng->getOpDomainRhsPipeline(), {NOSPACE});
463
466
467 pipeline_mng->getOpDomainLhsPipeline().push_back(
469 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
470 "FIELD1", "FIELD1", [](double, double, double) { return 1.; }));
471 pipeline_mng->getOpDomainLhsPipeline().push_back(
472 new OpSchurAssembleEnd<SCHUR_DSYSV>({}, {}, {}, {}, {}));
473
474 pipeline_mng->getOpDomainRhsPipeline().push_back(
475 new OpSource("FIELD1", ApproxFunctions::fUn));
476
477 auto integration_rule = [](int, int, int p_data) {
478 return 2 * p_data + 1;
479 };
480 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
481 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
482
484 };
485
486 auto solve_problem = [&] {
488 auto solver = pipeline_mng->createKSP();
489 CHKERR KSPSetFromOptions(solver);
490 CHKERR KSPSetUp(solver);
491
492 auto dm = simple->getDM();
493 auto D = createDMVector(dm);
494 auto F = vectorDuplicate(D);
495
496 CHKERR KSPSolve(solver, F, D);
497 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
498 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
499 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
500
502 };
503
504 auto check_solution = [&] {
506
507 auto ptr_values = boost::make_shared<VectorDouble>();
508 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
509
510 pipeline_mng->getDomainLhsFE().reset();
511 pipeline_mng->getOpDomainRhsPipeline().clear();
512
514 pipeline_mng->getOpDomainRhsPipeline(), {space});
515
516 pipeline_mng->getOpDomainRhsPipeline().push_back(
517 new OpValsDiffVals(vals, diff_vals, true));
518 pipeline_mng->getOpDomainRhsPipeline().push_back(
519 new OpCalculateScalarFieldValues("FIELD1", ptr_values));
520 pipeline_mng->getOpDomainRhsPipeline().push_back(
522 ptr_diff_vals));
523 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
524 vals, diff_vals, ptr_values, ptr_diff_vals, true));
525
526 CHKERR pipeline_mng->loopFiniteElements();
527
529 };
530
531 auto post_proc = [&] {
533
535
536 auto get_domain_post_proc_fe = [&](auto post_proc_mesh_ptr) {
537 auto post_proc_fe =
538 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<DomainEle>>(
539 m_field, post_proc_mesh_ptr);
540
542 post_proc_fe->getOpPtrVector(), {space});
543
544 auto ptr_values = boost::make_shared<VectorDouble>();
545 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
546
547 post_proc_fe->getOpPtrVector().push_back(
548 new OpCalculateScalarFieldValues("FIELD1", ptr_values));
549 post_proc_fe->getOpPtrVector().push_back(
551 ptr_diff_vals));
552
553 post_proc_fe->getOpPtrVector().push_back(
554
555 new OpPPMap(
556
557 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
558
559 {{"FIELD1", ptr_values}},
560
562
563 {}, {})
564
565 );
566 return post_proc_fe;
567 };
568
569 auto get_bdy_post_proc_fe = [&](auto post_proc_mesh_ptr) {
570 auto bdy_post_proc_fe =
571 boost::make_shared<PostProcBrokenMeshInMoabBaseCont<BoundaryEle>>(
572 m_field, post_proc_mesh_ptr);
573
574 auto op_loop_side = new OpLoopSide<EleOnSide>(
575 m_field, simple->getDomainFEName(), SPACE_DIM, Sev::noisy,
576 boost::make_shared<
578
579 auto ptr_values = boost::make_shared<VectorDouble>();
580 auto ptr_diff_vals = boost::make_shared<MatrixDouble>();
581
582 // push operators to side element
584 op_loop_side->getOpPtrVector(), {space});
585 op_loop_side->getOpPtrVector().push_back(
586 new OpCalculateScalarFieldValues("FIELD1", ptr_values));
587 op_loop_side->getOpPtrVector().push_back(
589 ptr_diff_vals));
590 // push op to boundary element
591 bdy_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
592
593 bdy_post_proc_fe->getOpPtrVector().push_back(
594
595 new OpPPMap(
596
597 bdy_post_proc_fe->getPostProcMesh(),
598 bdy_post_proc_fe->getMapGaussPts(),
599
600 {{"FIELD1", ptr_values}},
601
603
604 {}, {})
605
606 );
607
608 return bdy_post_proc_fe;
609 };
610
611 auto post_proc_mesh_ptr = boost::make_shared<moab::Core>();
612 auto post_proc_begin =
613 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
614 m_field, post_proc_mesh_ptr);
615 auto post_proc_end =
616 boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
617 m_field, post_proc_mesh_ptr);
618 auto domain_post_proc_fe = get_domain_post_proc_fe(post_proc_mesh_ptr);
619 auto bdy_post_proc_fe = get_bdy_post_proc_fe(post_proc_mesh_ptr);
620
621 CHKERR DMoFEMPreProcessFiniteElements(dm, post_proc_begin->getFEMethod());
622 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
623 domain_post_proc_fe);
624 CHKERR DMoFEMLoopFiniteElements(dm, simple->getBoundaryFEName(),
625 bdy_post_proc_fe);
626 CHKERR DMoFEMPostProcessFiniteElements(dm, post_proc_end->getFEMethod());
627 CHKERR post_proc_end->writeFile("out_post_proc.h5m");
628
630 };
631
632 CHKERR assemble_matrices_and_vectors();
633 CHKERR solve_problem();
634 CHKERR check_solution();
635 CHKERR post_proc();
636 }
638
640}
static Index< 'o', 3 > o
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
@ F
auto integration_rule
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:542
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
modify finite element table, only for advanced user
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
static int approx_order
static char help[]
constexpr int SPACE_DIM
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
static FTensor::Tensor1< double, 2 > diffFun(const double x, const double y, double z)
static double fUn(const double x, const double y, double z)
static FTensor::Tensor1< double, 3 > diffFun(const double x, const double y, double z)
static double fUn(const double x, const double y, double z)
Add operators pushing bases from local to physical configuration.
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
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Finite element data for entity.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Provide data structure for (tensor) field approximation.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Clear Schur complement internal data.
Definition: Schur.hpp:22
Assemble Schur complement.
Definition: Schur.hpp:104
PipelineManager interface.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FTensor::Index< 'i', 3 > i
OpCheckValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, boost::shared_ptr< VectorDouble > &ptr_vals, boost::shared_ptr< MatrixDouble > &ptr_diff_vals, bool check_grads)
boost::shared_ptr< VectorDouble > ptrVals
FTensor::Index< 'i', SPACE_DIM > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< MatrixDouble > ptrDiffVals
boost::shared_ptr< MatrixDouble > ptrVals
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpValsDiffVals(VectorDouble &vals, MatrixDouble &diff_vals, bool check_grads)
FTensor::Index< 'i', SPACE_DIM > i