v0.13.1
Loading...
Searching...
No Matches
hanging_node_approx.cpp
Go to the documentation of this file.
1/**
2 * \example hanging_node_approx.cpp
3 *
4 * Tetsing approximation with hanging nodes.
5 *
6 */
7
8
9
10#include <MoFEM.hpp>
11
12using namespace MoFEM;
13
14static char help[] = "...\n\n";
15
16constexpr char FIELD_NAME[] = "U";
17constexpr int FIELD_DIM = 1;
18constexpr int SPACE_DIM = 2;
19constexpr int nb_ref_levels = 3; ///< Three levels of refinement
20
21template <int DIM> struct ElementsAndOps {};
22
23template <> struct ElementsAndOps<2> {
30};
31
32template <> struct ElementsAndOps<3> {
35};
36
43
45
46template <int FIELD_DIM> struct ApproxFieldFunction;
47template <int FIELD_DIM> struct ApproxFieldFunctionDerivative;
48
49/**
50 * @brief third order polynomial used for testing
51 *
52 */
53template <> struct ApproxFieldFunction<1> {
54 auto operator()(const double x, const double y, const double z) {
55 return x * x + y * y + x * y * y + x * x * y;
56 }
57};
58
59/**
60 * @brief third order polynomial used for testing
61 *
62 */
63template <> struct ApproxFieldFunctionDerivative<1> {
64 auto operator()(const double x, const double y, const double z) {
65 // x * x + y * y + x * y * y + x * x * y
66
67 return FTensor::Tensor1<double, SPACE_DIM>{2 * x + y * y + 2 * x * y,
68 2 * y + 2 * x * y + x * x};
69 }
70};
71
72/**
73 * @brief evaluate mass matrix
74 *
75 */
78
79/**
80 * @brief evaluate source, i.e. rhs vector
81 *
82 */
85
86/**
87 * @brief set bit
88 *
89 */
90auto bit = [](auto l) { return BitRefLevel().set(l); };
91
92/**
93 * @brief set bit to marker
94 *
95 * Marker is used to mark field entities on skin on which we have hanging nodes
96 */
97auto marker = [](auto l) {
98 return BitRefLevel().set(BITREFLEVEL_SIZE - 1 - l);
99};
100
101/**
102 * @brief set levels of projection operators, which project field data from
103 * parent entities, to child, up to to level, i.e. last mesh refinement.
104 *
105 */
106template <typename PARENT_FE>
108 boost::shared_ptr<FEMethod> &fe_top,
109 ForcesAndSourcesCore::UserDataOperator::OpType op,
110 int verbosity, LogManager::SeverityLevel sev) {
111
112 auto jac_ptr = boost::make_shared<MatrixDouble>();
113 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
114 auto det_ptr = boost::make_shared<VectorDouble>();
115
116 BitRefLevel bit_marker;
117 for (auto l = 1; l <= nb_ref_levels; ++l)
118 bit_marker |= marker(l);
119
120 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
121 add_parent_level =
122 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
123 if (level > 0) {
124
125 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
126 new PARENT_FE(m_field));
127 if (op == DomainEleOp::OPSPACE) {
128 if (typeid(PARENT_FE) == typeid(DomainParentEle)) {
129 fe_ptr_current->getOpPtrVector().push_back(
130 new OpCalculateHOJac<2>(jac_ptr));
131 fe_ptr_current->getOpPtrVector().push_back(
132 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
133 fe_ptr_current->getOpPtrVector().push_back(
134 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
135 }
136 }
137
138 add_parent_level(
139 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
140 fe_ptr_current),
141 level - 1);
142
143 if (op == DomainEleOp::OPSPACE) {
144
145 parent_fe_pt->getOpPtrVector().push_back(
146
148
149 H1, op, fe_ptr_current,
150
151 BitRefLevel().set(), bit(0).flip(),
152
153 bit_marker, BitRefLevel().set(),
154
155 verbosity, sev));
156
157 } else {
158
159 parent_fe_pt->getOpPtrVector().push_back(
160
162
163 FIELD_NAME, op, fe_ptr_current,
164
165 BitRefLevel().set(), bit(0).flip(),
166
167 bit_marker, BitRefLevel().set(),
168
169 verbosity, sev));
170 }
171 }
172 };
173
174 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
176};
177
178/**
179 * @brief lambda function used to select elements on which finite element
180 * pipelines are executed.
181 *
182 * @note childs elements on pipeline, retrive data from parents using operators
183 * pushed by \ref set_parent_dofs
184 *
185 */
186auto test_bit_child = [](FEMethod *fe_ptr) {
187 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
189};
190
191struct AtomTest {
192
193 AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
194
196
197private:
200
203
204 /**
205 * @brief red mesh and randomly refine three times
206 *
207 * @return MoFEMErrorCode
208 */
210
211 /**
212 * @brief add field, and set up problem
213 *
214 * @return MoFEMErrorCode
215 */
221 struct CommonData {
222 boost::shared_ptr<VectorDouble> approxVals;
223 boost::shared_ptr<MatrixDouble> divApproxVals;
226 };
227
228 template <int FIELD_DIM> struct OpError;
229
230 template <int FIELD_DIM> struct OpErrorSkel;
231};
232
237template <> struct AtomTest::OpError<1> : public DomainEleOp {
238 boost::shared_ptr<CommonData> commonDataPtr;
239 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
240 : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
241 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
243
244 if (const size_t nb_dofs = data.getIndices().size()) {
245
247
248 const int nb_integration_pts = getGaussPts().size2();
249 auto t_w = getFTensor0IntegrationWeight();
250 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
251 auto t_grad_val =
252 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->divApproxVals));
253 auto t_coords = getFTensor1CoordsAtGaussPts();
254
255 VectorDouble nf(nb_dofs, false);
256 nf.clear();
257
258 const double volume = getMeasure();
259
260 auto t_row_base = data.getFTensor0N();
261 double error = 0;
262 for (int gg = 0; gg != nb_integration_pts; ++gg) {
263
264 const double alpha = t_w * volume;
265 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
266 t_coords(2));
267
268 auto t_grad_diff =
269 AtomTest::divApproxFunction(t_coords(0), t_coords(1), t_coords(2));
270 t_grad_diff(i) -= t_grad_val(i);
271
272 error += alpha * (pow(diff, 2) + t_grad_diff(i) * t_grad_diff(i));
273
274 for (size_t r = 0; r != nb_dofs; ++r) {
275 nf[r] += alpha * t_row_base * diff;
276 ++t_row_base;
277 }
278
279 ++t_w;
280 ++t_val;
281 ++t_grad_val;
282 ++t_coords;
283 }
284
285 const int index = 0;
286 CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
287 CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
288 }
289
291 }
292};
293
294template <> struct AtomTest::OpErrorSkel<1> : public BoundaryEleOp {
295 boost::shared_ptr<CommonData> commonDataPtr;
296 OpErrorSkel(boost::shared_ptr<CommonData> &common_data_ptr)
297 : BoundaryEleOp(H1, OPSPACE), commonDataPtr(common_data_ptr) {}
298 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
300
302
303 const int nb_integration_pts = getGaussPts().size2();
304 auto t_w = getFTensor0IntegrationWeight();
305 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
306 auto t_coords = getFTensor1CoordsAtGaussPts();
307
308 const double volume = getMeasure();
309
310 double error2 = 0;
311 for (int gg = 0; gg != nb_integration_pts; ++gg) {
312
313 const double alpha = t_w * volume;
314 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
315 t_coords(2));
316 error2 += alpha * (pow(diff, 2));
317
318 ++t_w;
319 ++t_val;
320 ++t_coords;
321 }
322
323 MOFEM_LOG("SELF", Sev::verbose) << "Boundary error " << sqrt(error2);
324
325 constexpr double eps = 1e-8;
326 if (sqrt(error2) > eps)
327 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
328 "Error on boundary = %6.4e", sqrt(error2));
329
330
332 }
333};
334
335//! [Run programme]
345}
346//! [Run programme]
347
348//! [Read mesh]
351 ParallelComm *pcomm =
352 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
353 Skinner skin(&mField.get_moab());
355
359
360 MOFEM_LOG("WORLD", Sev::verbose) << "Dim " << simpleInterface->getDim();
361
362 auto &moab = mField.get_moab();
363
364 Range level0_ents;
365 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
366 bit(0), BitRefLevel().set(), SPACE_DIM, level0_ents);
367 Range level0_skin;
368 CHKERR skin.find_skin(0, level0_ents, false, level0_skin);
369 CHKERR pcomm->filter_pstatus(level0_skin,
370 PSTATUS_SHARED | PSTATUS_MULTISHARED,
371 PSTATUS_NOT, -1, nullptr);
372
373 auto refine_mesh = [&](auto l) {
375
376 auto refine = mField.getInterface<MeshRefinement>();
377
378 auto meshset_level0_ptr = get_temp_meshset_ptr(moab);
379 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(l - 1), BitRefLevel().set(),
380 SPACE_DIM, *meshset_level0_ptr);
381
382 // random mesh refinement
383 auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
384
385 Range els;
386 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr, SPACE_DIM, els);
387 CHKERR bit_mng->filterEntitiesByRefLevel(bit(l - 1), bit(l - 1), els);
388
389 Range ele_to_refine;
390
391 if (l == 1) {
392 int ii = 0;
393 for (auto t : els) {
394 if ((ii % 2)) {
395 ele_to_refine.insert(t);
396 std::vector<EntityHandle> adj_edges;
397 CHKERR mField.get_moab().get_adjacencies(&t, 1, SPACE_DIM - 1, false,
398 adj_edges);
399 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
400 adj_edges.size());
401 }
402 ++ii;
403 }
404 } else {
405 Range level_skin;
406 CHKERR skin.find_skin(0, els, false, level_skin);
407 CHKERR pcomm->filter_pstatus(level_skin,
408 PSTATUS_SHARED | PSTATUS_MULTISHARED,
409 PSTATUS_NOT, -1, nullptr);
410 level_skin = subtract(level_skin, level0_skin);
411 Range adj;
412 CHKERR mField.get_moab().get_adjacencies(level_skin, SPACE_DIM, false,
413 adj, moab::Interface::UNION);
414 els = subtract(els, adj);
415 ele_to_refine.merge(els);
416 Range adj_edges;
417 CHKERR mField.get_moab().get_adjacencies(
418 els, SPACE_DIM - 1, false, adj_edges, moab::Interface::UNION);
419 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
420 }
421
422 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit(l),
423 false, VERBOSE);
424 CHKERR refine->refineTrisHangingNodes(*meshset_level0_ptr, bit(l), VERBOSE);
425 CHKERR bit_mng->updateRangeByChildren(level0_skin, level0_skin);
429 true);
430
431 CHKERR bit_mng->writeBitLevelByDim(
432 bit(l), BitRefLevel().set(), SPACE_DIM,
433 (boost::lexical_cast<std::string>(l) + "_ref_mesh.vtk").c_str(), "VTK",
434 "");
435 CHKERR bit_mng->writeBitLevelByDim(
436 bit(l), bit(l), MBTRI,
437 (boost::lexical_cast<std::string>(l) + "_only_ref_mesh.vtk").c_str(),
438 "VTK", "");
439
441 };
442
443 auto mark_skins = [&](auto l, auto m) {
445 Range ents;
447 ents);
448 Range level_skin;
449 CHKERR skin.find_skin(0, ents, false, level_skin);
450 CHKERR pcomm->filter_pstatus(level_skin,
451 PSTATUS_SHARED | PSTATUS_MULTISHARED,
452 PSTATUS_NOT, -1, nullptr);
453 level_skin = subtract(level_skin, level0_skin);
454 CHKERR mField.get_moab().get_adjacencies(level_skin, 0, false, level_skin,
455 moab::Interface::UNION);
456 CHKERR bit_mng->addBitRefLevel(level_skin, marker(m));
458 };
459
460 BitRefLevel bit_sum;
461 for (auto l = 0; l != nb_ref_levels; ++l) {
462 CHKERR refine_mesh(l + 1);
463 CHKERR mark_skins(l, l + 1);
464 CHKERR mark_skins(l + 1, l + 1);
465 bit_sum |= bit(l);
466 }
467 bit_sum |= bit(nb_ref_levels);
468
469 simpleInterface->getBitRefLevel() = bit_sum;
471
473}
474//! [Read mesh]
475
476//! [Set up problem]
479 // Add field
484
485 constexpr int order = 3;
487
488 // Simple interface will resolve adjacency to DOFs of parent of the element.
489 // Using that information MAtrixManager allocate appropriately size of
490 // matrix.
492 BitRefLevel bit_marker;
493 for (auto l = 1; l <= nb_ref_levels; ++l)
494 bit_marker |= marker(l);
495 simpleInterface->getBitAdjEnt() = bit_marker;
496
498
501
502 // remove obsolete DOFs from problem
503
504 for (int l = 0; l != nb_ref_levels; ++l) {
505 CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
506 FIELD_NAME, bit(l), bit(l));
507 CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
508 FIELD_NAME, marker(l + 1),
509 bit(l).flip());
510 }
512}
513//! [Set up problem]
514
515//! [Push operators to pipeline]
519
520 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
521
522 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
523 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
524
525 pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
526 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
527
528 auto beta = [](const double, const double, const double) { return 1; };
529 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
530 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
531 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
532 DomainEleOp::OPROW, QUIET, Sev::noisy);
533 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
534 DomainEleOp::OPCOL, QUIET, Sev::noisy);
535 pipeline_mng->getOpDomainLhsPipeline().push_back(
537
538 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
540 field_op_row->doWorkRhsHook = [](DataOperator *op_ptr, int side,
544 if (type == MBENTITYSET) {
545
546 MOFEM_LOG("SELF", Sev::verbose)
547 << "ROW: side/type: " << side << "/" << CN::EntityTypeName(type)
548 << " op space/base: " << FieldSpaceNames[data.getSpace()] << "/"
549 << ApproximationBaseNames[data.getBase()] << " DOFs "
550 << data.getIndices() << " nb base functions " << data.getN().size2()
551 << " nb base functions integration points " << data.getN().size1();
552
553 auto get_bool = [](auto fe, auto bit) {
554 return (bit & fe->getBitRefLevel()).any();
555 };
556
557 for (auto &field_ent : data.getFieldEntities()) {
558 MOFEM_LOG("SELF", Sev::verbose)
559 << "\t" << CN::EntityTypeName(field_ent->getEntType());
560 }
561 }
563 };
564
565 set_parent_dofs<DomainParentEle>(
567 Sev::verbose);
568 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
569 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
570 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
571
572 pipeline_mng->getOpDomainRhsPipeline().push_back(
574
576}
577//! [Push operators to pipeline]
578
579//! [Solve]
583
584 MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
585
586 auto solver = pipeline_mng->createKSP();
587 CHKERR KSPSetFromOptions(solver);
588 CHKERR KSPSetUp(solver);
589
590 auto dm = simpleInterface->getDM();
591 auto D = smartCreateDMVector(dm);
592 auto F = smartVectorDuplicate(D);
593
594 CHKERR KSPSolve(solver, F, D);
595 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
596 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
597 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
599}
600
601//! [Check results]
605 pipeline_mng->getDomainLhsFE().reset();
606 pipeline_mng->getDomainRhsFE().reset();
607 pipeline_mng->getBoundaryRhsFE().reset();
608
609 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
610 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
611 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
612 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
613 pipeline_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
614
615 auto common_data_ptr = boost::make_shared<CommonData>();
616 common_data_ptr->resVec = smartCreateDMVector(simpleInterface->getDM());
617 common_data_ptr->L2Vec = createSmartVectorMPI(
618 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
619 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
620 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
621
622 auto jac_ptr = boost::make_shared<MatrixDouble>();
623 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
624 auto det_ptr = boost::make_shared<VectorDouble>();
625
626 pipeline_mng->getOpDomainRhsPipeline().push_back(
627 new OpCalculateHOJac<2>(jac_ptr));
628 pipeline_mng->getOpDomainRhsPipeline().push_back(
629 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
630 pipeline_mng->getOpDomainRhsPipeline().push_back(
631 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
632
633 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
634 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
635 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
636 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
637
638 pipeline_mng->getOpDomainRhsPipeline().push_back(
640 FIELD_NAME, common_data_ptr->divApproxVals));
641 pipeline_mng->getOpDomainRhsPipeline().push_back(
643 common_data_ptr->approxVals));
644
645 pipeline_mng->getOpDomainRhsPipeline().push_back(
646 new OpError<FIELD_DIM>(common_data_ptr));
647
648 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
649 BoundaryEleOp::OPSPACE, QUIET, Sev::noisy);
650 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
651 BoundaryEleOp::OPROW, VERBOSE, Sev::noisy);
652 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
654 common_data_ptr->approxVals));
655 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
656 new OpErrorSkel<FIELD_DIM>(common_data_ptr));
657
658 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
659 CHKERR VecZeroEntries(common_data_ptr->resVec);
660
661 CHKERR pipeline_mng->loopFiniteElements();
662
663 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
664 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
665 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
666 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
667 double nrm2;
668 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
669 const double *array;
670 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
671 MOFEM_LOG_C("WORLD", Sev::inform, "Error %6.4e Vec norm %6.4e\n",
672 std::sqrt(array[0]), nrm2);
673
674 constexpr double eps = 1e-8;
675 if (nrm2 > eps)
676 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
677 "Not converged solution err = %6.4e", nrm2);
678 if (std::sqrt(array[0]) > eps)
679 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
680 "Error in approximation err = %6.4e", std::sqrt(array[0]));
681
682 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
683
685}
686//! [Check results]
687
691 pipeline_mng->getDomainLhsFE().reset();
692 pipeline_mng->getDomainRhsFE().reset();
693
694 auto rule = [](int, int, int p) -> int { return -1; };
695 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
696
697 static_cast<ForcesAndSourcesCore *>(pipeline_mng->getDomainRhsFE().get())
698 ->setRuleHook = [](ForcesAndSourcesCore *fe_raw_ptr, int order_row,
699 int order_col, int order_data) -> MoFEMErrorCode {
701 fe_raw_ptr->gaussPts.resize(3, 3);
702 fe_raw_ptr->gaussPts(0, 0) = 0;
703 fe_raw_ptr->gaussPts(1, 0) = 0;
704 fe_raw_ptr->gaussPts(2, 0) = 0;
705 fe_raw_ptr->gaussPts(0, 1) = 1;
706 fe_raw_ptr->gaussPts(1, 1) = 0;
707 fe_raw_ptr->gaussPts(2, 1) = 0;
708 fe_raw_ptr->gaussPts(0, 2) = 0;
709 fe_raw_ptr->gaussPts(1, 2) = 1;
710 fe_raw_ptr->gaussPts(2, 2) = 0;
712 };
713
714 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
716
717 auto approx_vals = boost::make_shared<VectorDouble>();
718
719 auto &moab = mField.get_moab();
720 Tag th;
721 double def_val[] = {0};
722 CHKERR moab.tag_get_handle("FIELD", 1, MB_TYPE_DOUBLE, th,
723 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
724
725 field_op_row->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
726 EntityType type,
729 if (type == MBVERTEX) {
730 auto op_ptr =
732 base_op_ptr);
733 auto t_field = getFTensor0FromVec(*approx_vals);
734 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
735 if (nb_gauss_pts != 3)
736 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
737 "Should be three guass pts.");
738 auto conn = op_ptr->getConn();
739 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
740 const double v = t_field;
741 CHKERR moab.tag_set_data(th, &conn[gg], 1, &v);
742 ++t_field;
743 }
744 }
746 };
747
748 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
749 DomainEleOp::OPSPACE, VERBOSE, Sev::noisy);
750 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
751 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
752 pipeline_mng->getOpDomainRhsPipeline().push_back(
753 new OpCalculateScalarFieldValues(FIELD_NAME, approx_vals));
754 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
755 CHKERR pipeline_mng->loopFiniteElements();
756
757 CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByType(
758 bit(nb_ref_levels), BitRefLevel().set(), MBTRI, "out.vtk", "VTK", "");
759
761}
762
763int main(int argc, char *argv[]) {
764
765 // Initialisation of MoFEM/PETSc and MOAB data structures
766 MoFEM::Core::Initialize(&argc, &argv, NULL, help);
767
768 try {
769
770 //! [Register MoFEM discrete manager in PETSc]
771 DMType dm_name = "DMMOFEM";
772 CHKERR DMRegister_MoFEM(dm_name);
773 //! [Register MoFEM discrete manager in PETSc
774
775 //! [Create MoAB]
776 moab::Core mb_instance; ///< mesh database
777 moab::Interface &moab = mb_instance; ///< mesh database interface
778 //! [Create MoAB]
779
780 //! [Create MoFEM]
781 MoFEM::Core core(moab); ///< finite element database
782 MoFEM::Interface &m_field = core; ///< finite element database insterface
783 //! [Create MoFEM]
784
785 //! [AtomTest]
786 AtomTest ex(m_field);
787 CHKERR ex.runProblem();
788 //! [AtomTest]
789 }
791
793}
static Index< 'p', 3 > p
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
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: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< 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 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 addBitRefLevel(const Range &ents, const BitRefLevel bit, int verb=QUIET) const
add bit ref level to ref entity
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childresn.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
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 temprary meshset.
Definition: Templates.hpp:1584
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()
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.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
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:26
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:264
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:373
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:289
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:780
auto & getBitAdjEnt()
bit ref level for parent
Definition: Simple.hpp:439
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:303
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:588
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:391
BitRefLevel & getBitRefLevelMask()
Get the BitRefLevel.
Definition: Simple.hpp:306
EntityHandle & getBoundaryMeshSet()
Get the BoundaryMeshSet object.
Definition: Simple.hpp:285
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:723
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:299
bool & getParentAdjacencies()
Get the addParentAdjacencies.
Definition: Simple.hpp:418
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.