v0.13.2
Loading...
Searching...
No Matches
hanging_node_approx.cpp
Go to the documentation of this file.
1/**
2 * \example hanging_node_approx.cpp
3 *
4 * Testing approximation with hanging nodes.
5 *
6 */
7
8#include <MoFEM.hpp>
9
10using namespace MoFEM;
11
12static char help[] = "...\n\n";
13
14constexpr char FIELD_NAME[] = "U";
15constexpr int FIELD_DIM = 1;
16constexpr int SPACE_DIM = 2;
17constexpr int nb_ref_levels = 3; ///< Three levels of refinement
18
28
30
31template <int FIELD_DIM> struct ApproxFieldFunction;
32template <int FIELD_DIM> struct ApproxFieldFunctionDerivative;
33
34/**
35 * @brief third order polynomial used for testing
36 *
37 */
38template <> struct ApproxFieldFunction<1> {
39 auto operator()(const double x, const double y, const double z) {
40 return x * x + y * y + x * y * y + x * x * y;
41 }
42};
43
44/**
45 * @brief third order polynomial used for testing
46 *
47 */
48template <> struct ApproxFieldFunctionDerivative<1> {
49 auto operator()(const double x, const double y, const double z) {
50 // x * x + y * y + x * y * y + x * x * y
51
52 return FTensor::Tensor1<double, SPACE_DIM>{2 * x + y * y + 2 * x * y,
53 2 * y + 2 * x * y + x * x};
54 }
55};
56
57/**
58 * @brief evaluate mass matrix
59 *
60 */
63
64/**
65 * @brief evaluate source, i.e. rhs vector
66 *
67 */
70
71/**
72 * @brief set bit
73 *
74 */
75auto bit = [](auto l) { return BitRefLevel().set(l); };
76
77/**
78 * @brief set bit to marker
79 *
80 * Marker is used to mark field entities on skin on which we have hanging nodes
81 */
82auto marker = [](auto l) {
83 return BitRefLevel().set(BITREFLEVEL_SIZE - 1 - l);
84};
85
86/**
87 * @brief set levels of projection operators, which project field data from
88 * parent entities, to child, up to to level, i.e. last mesh refinement.
89 *
90 */
91template <typename PARENT_FE>
93 boost::shared_ptr<FEMethod> &fe_top,
94 ForcesAndSourcesCore::UserDataOperator::OpType op,
95 int verbosity, LogManager::SeverityLevel sev) {
96
97 auto jac_ptr = boost::make_shared<MatrixDouble>();
98 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
99 auto det_ptr = boost::make_shared<VectorDouble>();
100
101 BitRefLevel bit_marker;
102 for (auto l = 1; l <= nb_ref_levels; ++l)
103 bit_marker |= marker(l);
104
105 /**
106 * @brief Collect data from parent elements to child
107 */
108 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
109 add_parent_level =
110 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
111
112 // Evaluate if not last parent element
113 if (level > 0) {
114
115 // Create domain parent FE
116 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
117 new PARENT_FE(m_field));
118 if (op == DomainEleOp::OPSPACE) {
119 // Push base function
120 if (typeid(PARENT_FE) == typeid(DomainParentEle))
122 fe_ptr_current->getOpPtrVector(), {H1});
123 }
124
125 // Call next level
126 add_parent_level(
127 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
128 fe_ptr_current),
129 level - 1);
130
131 // Add data to curent fe level
132 if (op == DomainEleOp::OPSPACE) {
133
134 // Only base
135 parent_fe_pt->getOpPtrVector().push_back(
136
138
139 H1, op, fe_ptr_current,
140
141 BitRefLevel().set(), bit(0).flip(),
142
143 bit_marker, BitRefLevel().set(),
144
145 verbosity, sev));
146
147 } else {
148
149 // Filed data
150 parent_fe_pt->getOpPtrVector().push_back(
151
153
154 FIELD_NAME, op, fe_ptr_current,
155
156 BitRefLevel().set(), bit(0).flip(),
157
158 bit_marker, BitRefLevel().set(),
159
160 verbosity, sev));
161 }
162 }
163 };
164
165 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
167};
168
169/**
170 * @brief lambda function used to select elements on which finite element
171 * pipelines are executed.
172 *
173 * @note childs elements on pipeline, retrieve data from parents using operators
174 * pushed by \ref set_parent_dofs
175 *
176 */
177auto test_bit_child = [](FEMethod *fe_ptr) {
178 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
180};
181
182struct AtomTest {
183
184 AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
185
187
188private:
191
194
195 /**
196 * @brief red mesh and randomly refine three times
197 *
198 * @return MoFEMErrorCode
199 */
201
202 /**
203 * @brief add field, and set up problem
204 *
205 * @return MoFEMErrorCode
206 */
212 struct CommonData {
213 boost::shared_ptr<VectorDouble> approxVals;
214 boost::shared_ptr<MatrixDouble> divApproxVals;
217 };
218
219 template <int FIELD_DIM> struct OpError;
220
221 template <int FIELD_DIM> struct OpErrorSkel;
222};
223
228template <> struct AtomTest::OpError<1> : public DomainEleOp {
229 boost::shared_ptr<CommonData> commonDataPtr;
230 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
231 : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
232 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
234
235 if (const size_t nb_dofs = data.getIndices().size()) {
236
238
239 const int nb_integration_pts = getGaussPts().size2();
240 auto t_w = getFTensor0IntegrationWeight();
241 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
242 auto t_grad_val =
243 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->divApproxVals));
244 auto t_coords = getFTensor1CoordsAtGaussPts();
245
246 VectorDouble nf(nb_dofs, false);
247 nf.clear();
248
249 const double volume = getMeasure();
250
251 auto t_row_base = data.getFTensor0N();
252 double error = 0;
253 for (int gg = 0; gg != nb_integration_pts; ++gg) {
254
255 const double alpha = t_w * volume;
256 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
257 t_coords(2));
258
259 auto t_grad_diff =
260 AtomTest::divApproxFunction(t_coords(0), t_coords(1), t_coords(2));
261 t_grad_diff(i) -= t_grad_val(i);
262
263 error += alpha * (pow(diff, 2) + t_grad_diff(i) * t_grad_diff(i));
264
265 for (size_t r = 0; r != nb_dofs; ++r) {
266 nf[r] += alpha * t_row_base * diff;
267 ++t_row_base;
268 }
269
270 ++t_w;
271 ++t_val;
272 ++t_grad_val;
273 ++t_coords;
274 }
275
276 const int index = 0;
277 CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
278 CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
279 }
280
282 }
283};
284
285template <> struct AtomTest::OpErrorSkel<1> : public BoundaryEleOp {
286 boost::shared_ptr<CommonData> commonDataPtr;
287 OpErrorSkel(boost::shared_ptr<CommonData> &common_data_ptr)
288 : BoundaryEleOp(H1, OPSPACE), commonDataPtr(common_data_ptr) {}
289 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
291
293
294 const int nb_integration_pts = getGaussPts().size2();
295 auto t_w = getFTensor0IntegrationWeight();
296 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
297 auto t_coords = getFTensor1CoordsAtGaussPts();
298
299 const double volume = getMeasure();
300
301 double error2 = 0;
302 for (int gg = 0; gg != nb_integration_pts; ++gg) {
303
304 const double alpha = t_w * volume;
305 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
306 t_coords(2));
307 error2 += alpha * (pow(diff, 2));
308
309 ++t_w;
310 ++t_val;
311 ++t_coords;
312 }
313
314 MOFEM_LOG("SELF", Sev::verbose) << "Boundary error " << sqrt(error2);
315
316 constexpr double eps = 1e-8;
317 if (sqrt(error2) > eps)
318 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
319 "Error on boundary = %6.4e", sqrt(error2));
320
322 }
323};
324
325//! [Run programme]
335}
336//! [Run programme]
337
338//! [Read mesh]
341 ParallelComm *pcomm =
342 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
343 Skinner skin(&mField.get_moab());
345
349
350 MOFEM_LOG("WORLD", Sev::verbose) << "Dim " << simpleInterface->getDim();
351
352 auto &moab = mField.get_moab();
353
354 Range level0_ents;
355 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
356 bit(0), BitRefLevel().set(), SPACE_DIM, level0_ents);
357 Range level0_skin;
358 CHKERR skin.find_skin(0, level0_ents, false, level0_skin);
359 CHKERR pcomm->filter_pstatus(level0_skin,
360 PSTATUS_SHARED | PSTATUS_MULTISHARED,
361 PSTATUS_NOT, -1, nullptr);
362
363 auto refine_mesh = [&](auto l) {
365
366 auto refine = mField.getInterface<MeshRefinement>();
367
368 auto meshset_level0_ptr = get_temp_meshset_ptr(moab);
369 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(l - 1), BitRefLevel().set(),
370 SPACE_DIM, *meshset_level0_ptr);
371
372 // random mesh refinement
373 auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
374
375 Range els;
376 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr, SPACE_DIM, els);
377 CHKERR bit_mng->filterEntitiesByRefLevel(bit(l - 1), bit(l - 1), els);
378
379 Range ele_to_refine;
380
381 if (l == 1) {
382 int ii = 0;
383 for (auto t : els) {
384 if ((ii % 2)) {
385 ele_to_refine.insert(t);
386 std::vector<EntityHandle> adj_edges;
387 CHKERR mField.get_moab().get_adjacencies(&t, 1, SPACE_DIM - 1, false,
388 adj_edges);
389 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
390 adj_edges.size());
391 }
392 ++ii;
393 }
394 } else {
395 Range level_skin;
396 CHKERR skin.find_skin(0, els, false, level_skin);
397 CHKERR pcomm->filter_pstatus(level_skin,
398 PSTATUS_SHARED | PSTATUS_MULTISHARED,
399 PSTATUS_NOT, -1, nullptr);
400 level_skin = subtract(level_skin, level0_skin);
401 Range adj;
402 CHKERR mField.get_moab().get_adjacencies(level_skin, SPACE_DIM, false,
403 adj, moab::Interface::UNION);
404 els = subtract(els, adj);
405 ele_to_refine.merge(els);
406 Range adj_edges;
407 CHKERR mField.get_moab().get_adjacencies(
408 els, SPACE_DIM - 1, false, adj_edges, moab::Interface::UNION);
409 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
410 }
411
412 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit(l),
413 false, VERBOSE);
414 CHKERR refine->refineTrisHangingNodes(*meshset_level0_ptr, bit(l), VERBOSE);
415 CHKERR bit_mng->updateRangeByChildren(level0_skin, level0_skin);
419 true);
420
421 CHKERR bit_mng->writeBitLevelByDim(
422 bit(l), BitRefLevel().set(), SPACE_DIM,
423 (boost::lexical_cast<std::string>(l) + "_ref_mesh.vtk").c_str(), "VTK",
424 "");
425 CHKERR bit_mng->writeBitLevelByDim(
426 bit(l), bit(l), MBTRI,
427 (boost::lexical_cast<std::string>(l) + "_only_ref_mesh.vtk").c_str(),
428 "VTK", "");
429
431 };
432
433 auto mark_skins = [&](auto l, auto m) {
435 Range ents;
437 ents);
438 Range level_skin;
439 CHKERR skin.find_skin(0, ents, false, level_skin);
440 CHKERR pcomm->filter_pstatus(level_skin,
441 PSTATUS_SHARED | PSTATUS_MULTISHARED,
442 PSTATUS_NOT, -1, nullptr);
443 level_skin = subtract(level_skin, level0_skin);
444 CHKERR mField.get_moab().get_adjacencies(level_skin, 0, false, level_skin,
445 moab::Interface::UNION);
446 CHKERR bit_mng->addBitRefLevel(level_skin, marker(m));
448 };
449
450 BitRefLevel bit_sum;
451 for (auto l = 0; l != nb_ref_levels; ++l) {
452 CHKERR refine_mesh(l + 1);
453 CHKERR mark_skins(l, l + 1);
454 CHKERR mark_skins(l + 1, l + 1);
455 bit_sum |= bit(l);
456 }
457 bit_sum |= bit(nb_ref_levels);
458
459 simpleInterface->getBitRefLevel() = bit_sum;
461
463}
464//! [Read mesh]
465
466//! [Set up problem]
469 // Add field
474
475 constexpr int order = 3;
477
478 // Simple interface will resolve adjacency to DOFs of parent of the element.
479 // Using that information MAtrixManager allocate appropriately size of
480 // matrix.
482
484
487
488 // remove obsolete DOFs from problem
489
490 for (int l = 0; l != nb_ref_levels; ++l) {
491 CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
492 FIELD_NAME, bit(l), bit(l));
493 CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
494 FIELD_NAME, marker(l + 1),
495 bit(l).flip());
496 }
498}
499//! [Set up problem]
500
501//! [Push operators to pipeline]
505
506 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
507
508 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
509 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
510
511 pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
512 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
513
514 auto beta = [](const double, const double, const double) { return 1; };
515 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
516 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
517 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
518 DomainEleOp::OPROW, QUIET, Sev::noisy);
519 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
520 DomainEleOp::OPCOL, QUIET, Sev::noisy);
521 pipeline_mng->getOpDomainLhsPipeline().push_back(
523
524 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
526 field_op_row->doWorkRhsHook = [](DataOperator *op_ptr, int side,
530 if (type == MBENTITYSET) {
531
532 MOFEM_LOG("SELF", Sev::verbose)
533 << "ROW: side/type: " << side << "/" << CN::EntityTypeName(type)
534 << " op space/base: " << FieldSpaceNames[data.getSpace()] << "/"
535 << ApproximationBaseNames[data.getBase()] << " DOFs "
536 << data.getIndices() << " nb base functions " << data.getN().size2()
537 << " nb base functions integration points " << data.getN().size1();
538
539 for (auto &field_ent : data.getFieldEntities()) {
540 MOFEM_LOG("SELF", Sev::verbose)
541 << "\t" << CN::EntityTypeName(field_ent->getEntType());
542 }
543 }
545 };
546
547 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
548 DomainEleOp::OPSPACE, VERBOSE, Sev::verbose);
549 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
550 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
551 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
552
553 pipeline_mng->getOpDomainRhsPipeline().push_back(
555
557}
558//! [Push operators to pipeline]
559
560//! [Solve]
564
565 MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
566
567 auto solver = pipeline_mng->createKSP();
568 CHKERR KSPSetFromOptions(solver);
569 CHKERR KSPSetUp(solver);
570
571 auto dm = simpleInterface->getDM();
572 auto D = smartCreateDMVector(dm);
573 auto F = smartVectorDuplicate(D);
574
575 CHKERR KSPSolve(solver, F, D);
576 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
577 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
578 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
580}
581
582//! [Check results]
586 pipeline_mng->getDomainLhsFE().reset();
587 pipeline_mng->getDomainRhsFE().reset();
588 pipeline_mng->getBoundaryRhsFE().reset();
589
590 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
591 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
592 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
593 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
594 pipeline_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
595
596 auto common_data_ptr = boost::make_shared<CommonData>();
597 common_data_ptr->resVec = smartCreateDMVector(simpleInterface->getDM());
598 common_data_ptr->L2Vec = createSmartVectorMPI(
599 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
600 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
601 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
602
604 pipeline_mng->getOpDomainRhsPipeline(), {H1});
605 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
606 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
607 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
608 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
609
610 pipeline_mng->getOpDomainRhsPipeline().push_back(
612 FIELD_NAME, common_data_ptr->divApproxVals));
613 pipeline_mng->getOpDomainRhsPipeline().push_back(
615 common_data_ptr->approxVals));
616
617 pipeline_mng->getOpDomainRhsPipeline().push_back(
618 new OpError<FIELD_DIM>(common_data_ptr));
619
620 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
621 BoundaryEleOp::OPSPACE, QUIET, Sev::noisy);
622 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
623 BoundaryEleOp::OPROW, VERBOSE, Sev::noisy);
624 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
626 common_data_ptr->approxVals));
627 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
628 new OpErrorSkel<FIELD_DIM>(common_data_ptr));
629
630 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
631 CHKERR VecZeroEntries(common_data_ptr->resVec);
632
633 CHKERR pipeline_mng->loopFiniteElements();
634
635 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
636 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
637 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
638 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
639 double nrm2;
640 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
641 const double *array;
642 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
643 MOFEM_LOG_C("WORLD", Sev::inform, "Error %6.4e Vec norm %6.4e\n",
644 std::sqrt(array[0]), nrm2);
645
646 constexpr double eps = 1e-8;
647 if (nrm2 > eps)
648 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
649 "Not converged solution err = %6.4e", nrm2);
650 if (std::sqrt(array[0]) > eps)
651 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
652 "Error in approximation err = %6.4e", std::sqrt(array[0]));
653
654 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
655
657}
658//! [Check results]
659
663 pipeline_mng->getDomainLhsFE().reset();
664 pipeline_mng->getDomainRhsFE().reset();
665
666 auto rule = [](int, int, int p) -> int { return -1; };
667 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
668
669 static_cast<ForcesAndSourcesCore *>(pipeline_mng->getDomainRhsFE().get())
670 ->setRuleHook = [](ForcesAndSourcesCore *fe_raw_ptr, int order_row,
671 int order_col, int order_data) -> MoFEMErrorCode {
673 fe_raw_ptr->gaussPts.resize(3, 3);
674 fe_raw_ptr->gaussPts(0, 0) = 0;
675 fe_raw_ptr->gaussPts(1, 0) = 0;
676 fe_raw_ptr->gaussPts(2, 0) = 0;
677 fe_raw_ptr->gaussPts(0, 1) = 1;
678 fe_raw_ptr->gaussPts(1, 1) = 0;
679 fe_raw_ptr->gaussPts(2, 1) = 0;
680 fe_raw_ptr->gaussPts(0, 2) = 0;
681 fe_raw_ptr->gaussPts(1, 2) = 1;
682 fe_raw_ptr->gaussPts(2, 2) = 0;
684 };
685
686 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
688
689 auto approx_vals = boost::make_shared<VectorDouble>();
690
691 auto &moab = mField.get_moab();
692 Tag th;
693 double def_val[] = {0};
694 CHKERR moab.tag_get_handle("FIELD", 1, MB_TYPE_DOUBLE, th,
695 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
696
697 field_op_row->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
698 EntityType type,
701 if (type == MBVERTEX) {
702 auto op_ptr =
704 base_op_ptr);
705 auto t_field = getFTensor0FromVec(*approx_vals);
706 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
707 if (nb_gauss_pts != 3)
708 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
709 "Should be three guass pts.");
710 auto conn = op_ptr->getConn();
711 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
712 const double v = t_field;
713 CHKERR moab.tag_set_data(th, &conn[gg], 1, &v);
714 ++t_field;
715 }
716 }
718 };
719
720 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
721 DomainEleOp::OPSPACE, VERBOSE, Sev::noisy);
722 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
723 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
724 pipeline_mng->getOpDomainRhsPipeline().push_back(
725 new OpCalculateScalarFieldValues(FIELD_NAME, approx_vals));
726 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
727 CHKERR pipeline_mng->loopFiniteElements();
728
729 CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByType(
730 bit(nb_ref_levels), BitRefLevel().set(), MBTRI, "out.vtk", "VTK", "");
731
733}
734
735int main(int argc, char *argv[]) {
736
737 // Initialisation of MoFEM/PETSc and MOAB data structures
738 MoFEM::Core::Initialize(&argc, &argv, NULL, help);
739
740 try {
741
742 //! [Register MoFEM discrete manager in PETSc]
743 DMType dm_name = "DMMOFEM";
744 CHKERR DMRegister_MoFEM(dm_name);
745 //! [Register MoFEM discrete manager in PETSc
746
747 //! [Create MoAB]
748 moab::Core mb_instance; ///< mesh database
749 moab::Interface &moab = mb_instance; ///< mesh database interface
750 //! [Create MoAB]
751
752 //! [Create MoFEM]
753 MoFEM::Core core(moab); ///< finite element database
754 MoFEM::Interface &m_field = core; ///< finite element database insterface
755 //! [Create MoFEM]
756
757 //! [AtomTest]
758 AtomTest ex(m_field);
759 CHKERR ex.runProblem();
760 //! [AtomTest]
761 }
763
765}
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
int main()
Definition: adol-c_atom.cpp:46
static const double eps
ElementsAndOps< SPACE_DIM >::DomianParentEle DomainParentEle
ElementsAndOps< SPACE_DIM >::BoundaryParentEle BoundaryParentEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ QUIET
Definition: definitions.h:208
@ VERBOSE
Definition: definitions.h:209
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
static const char *const FieldSpaceNames[]
Definition: definitions.h:92
#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
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#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
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'm', SPACE_DIM > m
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:511
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
MoFEMErrorCode updateMeshsetByEntitiesChildren(const EntityHandle parent, const BitRefLevel &parent_bit, const BitRefLevel &parent_mask, const BitRefLevel &child_bit, const BitRefLevel &child_mask, const EntityHandle child, EntityType child_type, const bool recursive=false, int verb=0)
Get child entities form meshset containing parent entities.
MoFEMErrorCode addBitRefLevel(const Range &ents, const BitRefLevel &bit, int verb=QUIET) const
add bit ref level to ref entity
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childrens.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
auto marker
set bit to marker
constexpr int nb_ref_levels
Three levels of refinement.
static char help[]
constexpr char FIELD_NAME[]
auto test_bit_child
lambda function used to select elements on which finite element pipelines are executed.
constexpr int SPACE_DIM
constexpr int FIELD_DIM
auto bit
set bit
auto set_parent_dofs(MoFEM::Interface &m_field, boost::shared_ptr< FEMethod > &fe_top, ForcesAndSourcesCore::UserDataOperator::OpType op, int verbosity, LogManager::SeverityLevel sev)
set levels of projection operators, which project field data from parent entities,...
FTensor::Index< 'i', SPACE_DIM > i
constexpr char FIELD_NAME[]
constexpr int SPACE_DIM
constexpr int FIELD_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< 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
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1626
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
constexpr double t
plate stiffness
Definition: plate.cpp:59
auto operator()(const double x, const double y, const double z)
auto operator()(const double x, const double y, const double z)
Collected data use d by operator to evaluate errors for the test.
boost::shared_ptr< VectorDouble > approxVals
SmartPetscObj< Vec > L2Vec
SmartPetscObj< Vec > resVec
boost::shared_ptr< MatrixDouble > divApproxVals
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Operator to evaluate errors.
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpErrorSkel(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
Simple * simpleInterface
MoFEMErrorCode checkResults()
[Check results]
MoFEMErrorCode solveSystem()
static ApproxFieldFunction< FIELD_DIM > approxFunction
AtomTest(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
add field, and set up problem
MoFEMErrorCode readMesh()
red mesh and randomly refine three times
MoFEMErrorCode assembleSystem()
static ApproxFieldFunctionDerivative< FIELD_DIM > divApproxFunction
MoFEM::Interface & mField
MoFEMErrorCode printResults()
[Check results]
MoFEMErrorCode runProblem()
Add operators pushing bases from local to physical configuration.
Managing BitRefLevels.
MoFEMErrorCode writeBitLevelByDim(const BitRefLevel bit, const BitRefLevel mask, const int dim, const char *file_name, const char *file_type, const char *options, const bool check_for_empty=true) const
Write bit ref level to file.
MoFEMErrorCode filterEntitiesByRefLevel(const BitRefLevel bit, const BitRefLevel mask, Range &ents, int verb=QUIET) const
filter entities by bit ref level
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure to get information form mofem into EntitiesFieldData
MatrixDouble gaussPts
Matrix of integration points.
Mesh refinement interface.
Operator to project base functions from parent entity to child.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:271
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
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:671
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:476
MoFEMErrorCode addBoundaryField(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 boundary.
Definition: Simple.cpp:282
BitRefLevel & getBitRefLevelMask()
Get the BitRefLevel.
Definition: Simple.hpp:320
EntityHandle & getBoundaryMeshSet()
Get the BoundaryMeshSet object.
Definition: Simple.hpp:299
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:614
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:348
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:313
bool & getParentAdjacencies()
Get the addParentAdjacencies.
Definition: Simple.hpp:432
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.