v0.14.0
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 BitRefLevel bit_marker;
98 for (auto l = 1; l <= nb_ref_levels; ++l)
99 bit_marker |= marker(l);
100
101 /**
102 * @brief Collect data from parent elements to child
103 */
104 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
105 add_parent_level =
106 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
107
108 // Evaluate if not last parent element
109 if (level > 0) {
110
111 // Create domain parent FE
112 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
113 new PARENT_FE(m_field));
114 if (op == DomainEleOp::OPSPACE) {
115 // Push base function
116 if (typeid(PARENT_FE) == typeid(DomainParentEle))
118 fe_ptr_current->getOpPtrVector(), {H1});
119 }
120
121 // Call next level
122 add_parent_level(
123 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
124 fe_ptr_current),
125 level - 1);
126
127 // Add data to curent fe level
128 if (op == DomainEleOp::OPSPACE) {
129
130 // Only base
131 parent_fe_pt->getOpPtrVector().push_back(
132
134
135 H1, op, fe_ptr_current,
136
137 BitRefLevel().set(), bit(0).flip(),
138
139 bit_marker, BitRefLevel().set(),
140
141 verbosity, sev));
142
143 } else {
144
145 // Filed data
146 parent_fe_pt->getOpPtrVector().push_back(
147
149
150 FIELD_NAME, op, fe_ptr_current,
151
152 BitRefLevel().set(), bit(0).flip(),
153
154 bit_marker, BitRefLevel().set(),
155
156 verbosity, sev));
157 }
158 }
159 };
160
161 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
163};
164
165/**
166 * @brief lambda function used to select elements on which finite element
167 * pipelines are executed.
168 *
169 * @note childs elements on pipeline, retrieve data from parents using operators
170 * pushed by \ref set_parent_dofs
171 *
172 */
173auto test_bit_child = [](FEMethod *fe_ptr) {
174 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
176};
177
178struct AtomTest {
179
180 AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
181
183
184private:
187
190
191 /**
192 * @brief red mesh and randomly refine three times
193 *
194 * @return MoFEMErrorCode
195 */
197
198 /**
199 * @brief add field, and set up problem
200 *
201 * @return MoFEMErrorCode
202 */
208 struct CommonData {
209 boost::shared_ptr<VectorDouble> approxVals;
210 boost::shared_ptr<MatrixDouble> divApproxVals;
213 };
214
215 template <int FIELD_DIM> struct OpError;
216
217 template <int FIELD_DIM> struct OpErrorSkel;
218};
219
224template <> struct AtomTest::OpError<1> : public DomainEleOp {
225 boost::shared_ptr<CommonData> commonDataPtr;
226 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
227 : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
228 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
230
231 if (const size_t nb_dofs = data.getIndices().size()) {
232
234
235 const int nb_integration_pts = getGaussPts().size2();
236 auto t_w = getFTensor0IntegrationWeight();
237 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
238 auto t_grad_val =
239 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->divApproxVals));
240 auto t_coords = getFTensor1CoordsAtGaussPts();
241
242 VectorDouble nf(nb_dofs, false);
243 nf.clear();
244
245 const double volume = getMeasure();
246
247 auto t_row_base = data.getFTensor0N();
248 double error = 0;
249 for (int gg = 0; gg != nb_integration_pts; ++gg) {
250
251 const double alpha = t_w * volume;
252 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
253 t_coords(2));
254
255 auto t_grad_diff =
256 AtomTest::divApproxFunction(t_coords(0), t_coords(1), t_coords(2));
257 t_grad_diff(i) -= t_grad_val(i);
258
259 error += alpha * (pow(diff, 2) + t_grad_diff(i) * t_grad_diff(i));
260
261 for (size_t r = 0; r != nb_dofs; ++r) {
262 nf[r] += alpha * t_row_base * diff;
263 ++t_row_base;
264 }
265
266 ++t_w;
267 ++t_val;
268 ++t_grad_val;
269 ++t_coords;
270 }
271
272 const int index = 0;
273 CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
274 CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
275 }
276
278 }
279};
280
281template <> struct AtomTest::OpErrorSkel<1> : public BoundaryEleOp {
282 boost::shared_ptr<CommonData> commonDataPtr;
283 OpErrorSkel(boost::shared_ptr<CommonData> &common_data_ptr)
284 : BoundaryEleOp(H1, OPSPACE), commonDataPtr(common_data_ptr) {}
285 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
287
289
290 const int nb_integration_pts = getGaussPts().size2();
291 auto t_w = getFTensor0IntegrationWeight();
292 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
293 auto t_coords = getFTensor1CoordsAtGaussPts();
294
295 const double volume = getMeasure();
296
297 double error2 = 0;
298 for (int gg = 0; gg != nb_integration_pts; ++gg) {
299
300 const double alpha = t_w * volume;
301 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
302 t_coords(2));
303 error2 += alpha * (pow(diff, 2));
304
305 ++t_w;
306 ++t_val;
307 ++t_coords;
308 }
309
310 MOFEM_LOG("SELF", Sev::verbose) << "Boundary error " << sqrt(error2);
311
312 constexpr double eps = 1e-8;
313 if (sqrt(error2) > eps)
314 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
315 "Error on boundary = %6.4e", sqrt(error2));
316
318 }
319};
320
321//! [Run programme]
331}
332//! [Run programme]
333
334//! [Read mesh]
337 ParallelComm *pcomm =
338 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
339 Skinner skin(&mField.get_moab());
341
345
346 MOFEM_LOG("WORLD", Sev::verbose) << "Dim " << simpleInterface->getDim();
347
348 auto &moab = mField.get_moab();
349
350 Range level0_ents;
351 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
352 bit(0), BitRefLevel().set(), SPACE_DIM, level0_ents);
353 Range level0_skin;
354 CHKERR skin.find_skin(0, level0_ents, false, level0_skin);
355 CHKERR pcomm->filter_pstatus(level0_skin,
356 PSTATUS_SHARED | PSTATUS_MULTISHARED,
357 PSTATUS_NOT, -1, nullptr);
358
359 auto refine_mesh = [&](auto l) {
361
362 auto refine = mField.getInterface<MeshRefinement>();
363
364 auto meshset_level0_ptr = get_temp_meshset_ptr(moab);
365 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(l - 1), BitRefLevel().set(),
366 SPACE_DIM, *meshset_level0_ptr);
367
368 // random mesh refinement
369 auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
370
371 Range els;
372 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr, SPACE_DIM, els);
373 CHKERR bit_mng->filterEntitiesByRefLevel(bit(l - 1), bit(l - 1), els);
374
375 Range ele_to_refine;
376
377 if (l == 1) {
378 int ii = 0;
379 for (auto t : els) {
380 if ((ii % 2)) {
381 ele_to_refine.insert(t);
382 std::vector<EntityHandle> adj_edges;
383 CHKERR mField.get_moab().get_adjacencies(&t, 1, SPACE_DIM - 1, false,
384 adj_edges);
385 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
386 adj_edges.size());
387 }
388 ++ii;
389 }
390 } else {
391 Range level_skin;
392 CHKERR skin.find_skin(0, els, false, level_skin);
393 CHKERR pcomm->filter_pstatus(level_skin,
394 PSTATUS_SHARED | PSTATUS_MULTISHARED,
395 PSTATUS_NOT, -1, nullptr);
396 level_skin = subtract(level_skin, level0_skin);
397 Range adj;
398 CHKERR mField.get_moab().get_adjacencies(level_skin, SPACE_DIM, false,
399 adj, moab::Interface::UNION);
400 els = subtract(els, adj);
401 ele_to_refine.merge(els);
402 Range adj_edges;
403 CHKERR mField.get_moab().get_adjacencies(
404 els, SPACE_DIM - 1, false, adj_edges, moab::Interface::UNION);
405 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
406 }
407
408 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit(l),
409 false, VERBOSE);
410 CHKERR refine->refineTrisHangingNodes(*meshset_level0_ptr, bit(l), VERBOSE);
411 CHKERR bit_mng->updateRangeByChildren(level0_skin, level0_skin);
415 true);
416
417 CHKERR bit_mng->writeBitLevelByDim(
418 bit(l), BitRefLevel().set(), SPACE_DIM,
419 (boost::lexical_cast<std::string>(l) + "_ref_mesh.vtk").c_str(), "VTK",
420 "");
421 CHKERR bit_mng->writeBitLevelByDim(
422 bit(l), bit(l), MBTRI,
423 (boost::lexical_cast<std::string>(l) + "_only_ref_mesh.vtk").c_str(),
424 "VTK", "");
425
427 };
428
429 auto mark_skins = [&](auto l, auto m) {
431 Range ents;
433 ents);
434 Range level_skin;
435 CHKERR skin.find_skin(0, ents, false, level_skin);
436 CHKERR pcomm->filter_pstatus(level_skin,
437 PSTATUS_SHARED | PSTATUS_MULTISHARED,
438 PSTATUS_NOT, -1, nullptr);
439 level_skin = subtract(level_skin, level0_skin);
440 CHKERR mField.get_moab().get_adjacencies(level_skin, 0, false, level_skin,
441 moab::Interface::UNION);
442 CHKERR bit_mng->addBitRefLevel(level_skin, marker(m));
444 };
445
446 BitRefLevel bit_sum;
447 for (auto l = 0; l != nb_ref_levels; ++l) {
448 CHKERR refine_mesh(l + 1);
449 CHKERR mark_skins(l, l + 1);
450 CHKERR mark_skins(l + 1, l + 1);
451 bit_sum |= bit(l);
452 }
453 bit_sum |= bit(nb_ref_levels);
454
455 simpleInterface->getBitRefLevel() = bit_sum;
457
459}
460//! [Read mesh]
461
462//! [Set up problem]
465 // Add field
470
471 constexpr int order = 3;
473
474 // Simple interface will resolve adjacency to DOFs of parent of the element.
475 // Using that information MAtrixManager allocate appropriately size of
476 // matrix.
478
480
482
483 // remove obsolete DOFs from problem
484
485 for (int l = 0; l != nb_ref_levels; ++l) {
487 FIELD_NAME, bit(l), bit(l));
489 FIELD_NAME, marker(l + 1),
490 bit(l).flip());
491 }
493}
494//! [Set up problem]
495
496//! [Push operators to pipeline]
500
501 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
502
503 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
504 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
505
506 pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
507 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
508
509 auto beta = [](const double, const double, const double) { return 1; };
510 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
511 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
512 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
513 DomainEleOp::OPROW, QUIET, Sev::noisy);
514 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
515 DomainEleOp::OPCOL, QUIET, Sev::noisy);
516 pipeline_mng->getOpDomainLhsPipeline().push_back(
518
519 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
521 field_op_row->doWorkRhsHook = [](DataOperator *op_ptr, int side,
525 if (type == MBENTITYSET) {
526
527 MOFEM_LOG("SELF", Sev::verbose)
528 << "ROW: side/type: " << side << "/" << CN::EntityTypeName(type)
529 << " op space/base: " << FieldSpaceNames[data.getSpace()] << "/"
530 << ApproximationBaseNames[data.getBase()] << " DOFs "
531 << data.getIndices() << " nb base functions " << data.getN().size2()
532 << " nb base functions integration points " << data.getN().size1();
533
534 for (auto &field_ent : data.getFieldEntities()) {
535 MOFEM_LOG("SELF", Sev::verbose)
536 << "\t" << CN::EntityTypeName(field_ent->getEntType());
537 }
538 }
540 };
541
542 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
543 DomainEleOp::OPSPACE, NOISY, Sev::verbose);
544 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
545 DomainEleOp::OPROW, NOISY, Sev::noisy);
546 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
547
548 pipeline_mng->getOpDomainRhsPipeline().push_back(
550
552}
553//! [Push operators to pipeline]
554
555//! [Solve]
559
560 MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
561
562 auto solver = pipeline_mng->createKSP();
563 CHKERR KSPSetFromOptions(solver);
564 CHKERR KSPSetUp(solver);
565
566 auto dm = simpleInterface->getDM();
567 auto D = createDMVector(dm);
568 auto F = vectorDuplicate(D);
569
570 CHKERR KSPSolve(solver, F, D);
571 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
572 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
573 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
575}
576
577//! [Check results]
581 pipeline_mng->getDomainLhsFE().reset();
582 pipeline_mng->getDomainRhsFE().reset();
583 pipeline_mng->getBoundaryRhsFE().reset();
584
585 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
586 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
587 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
588 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
589 pipeline_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
590
591 auto common_data_ptr = boost::make_shared<CommonData>();
592 common_data_ptr->resVec = createDMVector(simpleInterface->getDM());
593 common_data_ptr->L2Vec = createVectorMPI(
594 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
595 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
596 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
597
599 pipeline_mng->getOpDomainRhsPipeline(), {H1});
600 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
601 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
602 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
603 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
604
605 pipeline_mng->getOpDomainRhsPipeline().push_back(
607 FIELD_NAME, common_data_ptr->divApproxVals));
608 pipeline_mng->getOpDomainRhsPipeline().push_back(
610 common_data_ptr->approxVals));
611
612 pipeline_mng->getOpDomainRhsPipeline().push_back(
613 new OpError<FIELD_DIM>(common_data_ptr));
614
615 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
616 BoundaryEleOp::OPSPACE, QUIET, Sev::noisy);
617 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
618 BoundaryEleOp::OPROW, VERBOSE, Sev::noisy);
619 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
621 common_data_ptr->approxVals));
622 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
623 new OpErrorSkel<FIELD_DIM>(common_data_ptr));
624
625 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
626 CHKERR VecZeroEntries(common_data_ptr->resVec);
627
628 CHKERR pipeline_mng->loopFiniteElements();
629
630 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
631 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
632 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
633 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
634 double nrm2;
635 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
636 const double *array;
637 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
638 MOFEM_LOG_C("WORLD", Sev::inform, "Error %6.4e Vec norm %6.4e\n",
639 std::sqrt(array[0]), nrm2);
640
641 constexpr double eps = 1e-8;
642 if (nrm2 > eps)
643 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
644 "Not converged solution err = %6.4e", nrm2);
645 if (std::sqrt(array[0]) > eps)
646 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
647 "Error in approximation err = %6.4e", std::sqrt(array[0]));
648
649 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
650
652}
653//! [Check results]
654
658 pipeline_mng->getDomainLhsFE().reset();
659 pipeline_mng->getDomainRhsFE().reset();
660
661 auto rule = [](int, int, int p) -> int { return -1; };
662 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
663
664 static_cast<ForcesAndSourcesCore *>(pipeline_mng->getDomainRhsFE().get())
665 ->setRuleHook = [](ForcesAndSourcesCore *fe_raw_ptr, int order_row,
666 int order_col, int order_data) -> MoFEMErrorCode {
668 fe_raw_ptr->gaussPts.resize(3, 3);
669 fe_raw_ptr->gaussPts(0, 0) = 0;
670 fe_raw_ptr->gaussPts(1, 0) = 0;
671 fe_raw_ptr->gaussPts(2, 0) = 0;
672 fe_raw_ptr->gaussPts(0, 1) = 1;
673 fe_raw_ptr->gaussPts(1, 1) = 0;
674 fe_raw_ptr->gaussPts(2, 1) = 0;
675 fe_raw_ptr->gaussPts(0, 2) = 0;
676 fe_raw_ptr->gaussPts(1, 2) = 1;
677 fe_raw_ptr->gaussPts(2, 2) = 0;
679 };
680
681 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
683
684 auto approx_vals = boost::make_shared<VectorDouble>();
685
686 auto &moab = mField.get_moab();
687 Tag th;
688 double def_val[] = {0};
689 CHKERR moab.tag_get_handle("FIELD", 1, MB_TYPE_DOUBLE, th,
690 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
691
692 field_op_row->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
693 EntityType type,
696 if (type == MBVERTEX) {
697 auto op_ptr =
698 static_cast<FaceElementForcesAndSourcesCore::UserDataOperator *>(
699 base_op_ptr);
700 auto t_field = getFTensor0FromVec(*approx_vals);
701 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
702 if (nb_gauss_pts != 3)
703 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
704 "Should be three guass pts.");
705 auto conn = op_ptr->getConn();
706 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
707 const double v = t_field;
708 CHKERR moab.tag_set_data(th, &conn[gg], 1, &v);
709 ++t_field;
710 }
711 }
713 };
714
715 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
716 DomainEleOp::OPSPACE, VERBOSE, Sev::noisy);
717 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
718 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
719 pipeline_mng->getOpDomainRhsPipeline().push_back(
720 new OpCalculateScalarFieldValues(FIELD_NAME, approx_vals));
721 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
722 CHKERR pipeline_mng->loopFiniteElements();
723
724 CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByType(
725 bit(nb_ref_levels), BitRefLevel().set(), MBTRI, "out.vtk", "VTK", "");
726
728}
729
730int main(int argc, char *argv[]) {
731
732 // Initialisation of MoFEM/PETSc and MOAB data structures
733 MoFEM::Core::Initialize(&argc, &argv, NULL, help);
734
735 try {
736
737 //! [Register MoFEM discrete manager in PETSc]
738 DMType dm_name = "DMMOFEM";
739 CHKERR DMRegister_MoFEM(dm_name);
740 //! [Register MoFEM discrete manager in PETSc
741
742 //! [Create MoAB]
743 moab::Core mb_instance; ///< mesh database
744 moab::Interface &moab = mb_instance; ///< mesh database interface
745 //! [Create MoAB]
746
747 //! [Create MoFEM]
748 MoFEM::Core core(moab); ///< finite element database
749 MoFEM::Interface &m_field = core; ///< finite element database insterface
750 //! [Create MoFEM]
751
752 //! [AtomTest]
753 AtomTest ex(m_field);
754 CHKERR ex.runProblem();
755 //! [AtomTest]
756 }
758
760}
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
int main()
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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ QUIET
@ NOISY
@ VERBOSE
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define BITREFLEVEL_SIZE
max number of refinements
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
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 ...
@ 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()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
FTensor::Index< 'm', SPACE_DIM > m
@ F
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
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1003
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.
SeverityLevel
Severity levels.
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
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.
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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
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()
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:285
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:194
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:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:473
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:334
EntityHandle & getBoundaryMeshSet()
Get the BoundaryMeshSet object.
Definition Simple.hpp:313
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:611
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:362
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition Simple.hpp:327
bool & getParentAdjacencies()
Get the addParentAdjacencies.
Definition Simple.hpp:446
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.