v0.13.1
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/* This file is part of MoFEM.
9 * MoFEM is free software: you can redistribute it and/or modify it under
10 * the terms of the GNU Lesser General Public License as published by the
11 * Free Software Foundation, either version 3 of the License, or (at your
12 * option) any later version.
13 *
14 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 * License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21
22#include <MoFEM.hpp>
23
24using namespace MoFEM;
25
26static char help[] = "...\n\n";
27
28constexpr char FIELD_NAME[] = "U";
29constexpr int FIELD_DIM = 1;
30constexpr int SPACE_DIM = 2;
31constexpr int nb_ref_levels = 3; ///< Three levels of refinement
32
33template <int DIM> struct ElementsAndOps {};
34
35template <> struct ElementsAndOps<2> {
42};
43
44template <> struct ElementsAndOps<3> {
47};
48
55
57
58template <int FIELD_DIM> struct ApproxFieldFunction;
59template <int FIELD_DIM> struct ApproxFieldFunctionDerivative;
60
61/**
62 * @brief third order polynomial used for testing
63 *
64 */
65template <> struct ApproxFieldFunction<1> {
66 auto operator()(const double x, const double y, const double z) {
67 return x * x + y * y + x * y * y + x * x * y;
68 }
69};
70
71/**
72 * @brief third order polynomial used for testing
73 *
74 */
75template <> struct ApproxFieldFunctionDerivative<1> {
76 auto operator()(const double x, const double y, const double z) {
77 // x * x + y * y + x * y * y + x * x * y
78
79 return FTensor::Tensor1<double, SPACE_DIM>{2 * x + y * y + 2 * x * y,
80 2 * y + 2 * x * y + x * x};
81 }
82};
83
84/**
85 * @brief evaluate mass matrix
86 *
87 */
90
91/**
92 * @brief evaluate source, i.e. rhs vector
93 *
94 */
97
98/**
99 * @brief set bit
100 *
101 */
102auto bit = [](auto l) { return BitRefLevel().set(l); };
103
104/**
105 * @brief set bit to marker
106 *
107 * Marker is used to mark field entities on skin on which we have hanging nodes
108 */
109auto marker = [](auto l) {
110 return BitRefLevel().set(BITREFLEVEL_SIZE - 1 - l);
111};
112
113/**
114 * @brief set levels of projection operators, which project field data from
115 * parent entities, to child, up to to level, i.e. last mesh refinement.
116 *
117 */
118template <typename PARENT_FE>
120 boost::shared_ptr<FEMethod> &fe_top,
121 ForcesAndSourcesCore::UserDataOperator::OpType op,
122 int verbosity, LogManager::SeverityLevel sev) {
123
124 auto jac_ptr = boost::make_shared<MatrixDouble>();
125 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
126 auto det_ptr = boost::make_shared<VectorDouble>();
127
128 BitRefLevel bit_marker;
129 for (auto l = 1; l <= nb_ref_levels; ++l)
130 bit_marker |= marker(l);
131
132 boost::function<void(boost::shared_ptr<ForcesAndSourcesCore>, int)>
133 add_parent_level =
134 [&](boost::shared_ptr<ForcesAndSourcesCore> parent_fe_pt, int level) {
135 if (level > 0) {
136
137 auto fe_ptr_current = boost::shared_ptr<ForcesAndSourcesCore>(
138 new PARENT_FE(m_field));
139 if (op == DomainEleOp::OPSPACE) {
140 if (typeid(PARENT_FE) == typeid(DomainParentEle)) {
141 fe_ptr_current->getOpPtrVector().push_back(
142 new OpCalculateHOJac<2>(jac_ptr));
143 fe_ptr_current->getOpPtrVector().push_back(
144 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
145 fe_ptr_current->getOpPtrVector().push_back(
146 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
147 }
148 }
149
150 add_parent_level(
151 boost::dynamic_pointer_cast<ForcesAndSourcesCore>(
152 fe_ptr_current),
153 level - 1);
154
155 if (op == DomainEleOp::OPSPACE) {
156
157 parent_fe_pt->getOpPtrVector().push_back(
158
160
161 H1, op, fe_ptr_current,
162
163 BitRefLevel().set(), bit(0).flip(),
164
165 bit_marker, BitRefLevel().set(),
166
167 verbosity, sev));
168
169 } else {
170
171 parent_fe_pt->getOpPtrVector().push_back(
172
174
175 FIELD_NAME, op, fe_ptr_current,
176
177 BitRefLevel().set(), bit(0).flip(),
178
179 bit_marker, BitRefLevel().set(),
180
181 verbosity, sev));
182 }
183 }
184 };
185
186 add_parent_level(boost::dynamic_pointer_cast<ForcesAndSourcesCore>(fe_top),
188};
189
190/**
191 * @brief lambda function used to select elements on which finite element
192 * pipelines are executed.
193 *
194 * @note childs elements on pipeline, retrive data from parents using operators
195 * pushed by \ref set_parent_dofs
196 *
197 */
198auto test_bit_child = [](FEMethod *fe_ptr) {
199 return fe_ptr->numeredEntFiniteElementPtr->getBitRefLevel().test(
201};
202
203struct AtomTest {
204
205 AtomTest(MoFEM::Interface &m_field) : mField(m_field) {}
206
208
209private:
212
215
216 /**
217 * @brief red mesh and randomly refine three times
218 *
219 * @return MoFEMErrorCode
220 */
222
223 /**
224 * @brief add field, and set up problem
225 *
226 * @return MoFEMErrorCode
227 */
233 struct CommonData {
234 boost::shared_ptr<VectorDouble> approxVals;
235 boost::shared_ptr<MatrixDouble> divApproxVals;
238 };
239
240 template <int FIELD_DIM> struct OpError;
241
242 template <int FIELD_DIM> struct OpErrorSkel;
243};
244
249template <> struct AtomTest::OpError<1> : public DomainEleOp {
250 boost::shared_ptr<CommonData> commonDataPtr;
251 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
252 : DomainEleOp(FIELD_NAME, OPROW), commonDataPtr(common_data_ptr) {}
255
256 if (const size_t nb_dofs = data.getIndices().size()) {
257
259
260 const int nb_integration_pts = getGaussPts().size2();
261 auto t_w = getFTensor0IntegrationWeight();
262 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
263 auto t_grad_val =
264 getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->divApproxVals));
265 auto t_coords = getFTensor1CoordsAtGaussPts();
266
267 VectorDouble nf(nb_dofs, false);
268 nf.clear();
269
270 const double volume = getMeasure();
271
272 auto t_row_base = data.getFTensor0N();
273 double error = 0;
274 for (int gg = 0; gg != nb_integration_pts; ++gg) {
275
276 const double alpha = t_w * volume;
277 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
278 t_coords(2));
279
280 auto t_grad_diff =
281 AtomTest::divApproxFunction(t_coords(0), t_coords(1), t_coords(2));
282 t_grad_diff(i) -= t_grad_val(i);
283
284 error += alpha * (pow(diff, 2) + t_grad_diff(i) * t_grad_diff(i));
285
286 for (size_t r = 0; r != nb_dofs; ++r) {
287 nf[r] += alpha * t_row_base * diff;
288 ++t_row_base;
289 }
290
291 ++t_w;
292 ++t_val;
293 ++t_grad_val;
294 ++t_coords;
295 }
296
297 const int index = 0;
298 CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
299 CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
300 }
301
303 }
304};
305
306template <> struct AtomTest::OpErrorSkel<1> : public BoundaryEleOp {
307 boost::shared_ptr<CommonData> commonDataPtr;
308 OpErrorSkel(boost::shared_ptr<CommonData> &common_data_ptr)
309 : BoundaryEleOp(H1, OPSPACE), commonDataPtr(common_data_ptr) {}
312
314
315 const int nb_integration_pts = getGaussPts().size2();
316 auto t_w = getFTensor0IntegrationWeight();
317 auto t_val = getFTensor0FromVec(*(commonDataPtr->approxVals));
318 auto t_coords = getFTensor1CoordsAtGaussPts();
319
320 const double volume = getMeasure();
321
322 double error2 = 0;
323 for (int gg = 0; gg != nb_integration_pts; ++gg) {
324
325 const double alpha = t_w * volume;
326 double diff = t_val - AtomTest::approxFunction(t_coords(0), t_coords(1),
327 t_coords(2));
328 error2 += alpha * (pow(diff, 2));
329
330 ++t_w;
331 ++t_val;
332 ++t_coords;
333 }
334
335 MOFEM_LOG("SELF", Sev::verbose) << "Boundary error " << sqrt(error2);
336
337 constexpr double eps = 1e-8;
338 if (sqrt(error2) > eps)
339 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
340 "Error on boundary = %6.4e", sqrt(error2));
341
342
344 }
345};
346
347//! [Run programme]
357}
358//! [Run programme]
359
360//! [Read mesh]
363 ParallelComm *pcomm =
364 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
365 Skinner skin(&mField.get_moab());
367
371
372 MOFEM_LOG("WORLD", Sev::verbose) << "Dim " << simpleInterface->getDim();
373
374 auto &moab = mField.get_moab();
375
376 Range level0_ents;
377 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
378 bit(0), BitRefLevel().set(), SPACE_DIM, level0_ents);
379 Range level0_skin;
380 CHKERR skin.find_skin(0, level0_ents, false, level0_skin);
381 CHKERR pcomm->filter_pstatus(level0_skin,
382 PSTATUS_SHARED | PSTATUS_MULTISHARED,
383 PSTATUS_NOT, -1, nullptr);
384
385 auto refine_mesh = [&](auto l) {
387
388 auto refine = mField.getInterface<MeshRefinement>();
389
390 auto meshset_level0_ptr = get_temp_meshset_ptr(moab);
391 CHKERR bit_mng->getEntitiesByDimAndRefLevel(bit(l - 1), BitRefLevel().set(),
392 SPACE_DIM, *meshset_level0_ptr);
393
394 // random mesh refinement
395 auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
396
397 Range els;
398 CHKERR moab.get_entities_by_dimension(*meshset_level0_ptr, SPACE_DIM, els);
399 CHKERR bit_mng->filterEntitiesByRefLevel(bit(l - 1), bit(l - 1), els);
400
401 Range ele_to_refine;
402
403 if (l == 1) {
404 int ii = 0;
405 for (auto t : els) {
406 if ((ii % 2)) {
407 ele_to_refine.insert(t);
408 std::vector<EntityHandle> adj_edges;
409 CHKERR mField.get_moab().get_adjacencies(&t, 1, SPACE_DIM - 1, false,
410 adj_edges);
411 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_edges.begin(),
412 adj_edges.size());
413 }
414 ++ii;
415 }
416 } else {
417 Range level_skin;
418 CHKERR skin.find_skin(0, els, false, level_skin);
419 CHKERR pcomm->filter_pstatus(level_skin,
420 PSTATUS_SHARED | PSTATUS_MULTISHARED,
421 PSTATUS_NOT, -1, nullptr);
422 level_skin = subtract(level_skin, level0_skin);
423 Range adj;
424 CHKERR mField.get_moab().get_adjacencies(level_skin, SPACE_DIM, false,
425 adj, moab::Interface::UNION);
426 els = subtract(els, adj);
427 ele_to_refine.merge(els);
428 Range adj_edges;
429 CHKERR mField.get_moab().get_adjacencies(
430 els, SPACE_DIM - 1, false, adj_edges, moab::Interface::UNION);
431 CHKERR moab.add_entities(*meshset_ref_edges_ptr, adj_edges);
432 }
433
434 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit(l),
435 false, VERBOSE);
436 CHKERR refine->refineTrisHangingNodes(*meshset_level0_ptr, bit(l), VERBOSE);
437 CHKERR bit_mng->updateRangeByChildren(level0_skin, level0_skin);
441 true);
442
443 CHKERR bit_mng->writeBitLevelByDim(
444 bit(l), BitRefLevel().set(), SPACE_DIM,
445 (boost::lexical_cast<std::string>(l) + "_ref_mesh.vtk").c_str(), "VTK",
446 "");
447 CHKERR bit_mng->writeBitLevelByDim(
448 bit(l), bit(l), MBTRI,
449 (boost::lexical_cast<std::string>(l) + "_only_ref_mesh.vtk").c_str(),
450 "VTK", "");
451
453 };
454
455 auto mark_skins = [&](auto l, auto m) {
457 Range ents;
459 ents);
460 Range level_skin;
461 CHKERR skin.find_skin(0, ents, false, level_skin);
462 CHKERR pcomm->filter_pstatus(level_skin,
463 PSTATUS_SHARED | PSTATUS_MULTISHARED,
464 PSTATUS_NOT, -1, nullptr);
465 level_skin = subtract(level_skin, level0_skin);
466 CHKERR mField.get_moab().get_adjacencies(level_skin, 0, false, level_skin,
467 moab::Interface::UNION);
468 CHKERR bit_mng->addBitRefLevel(level_skin, marker(m));
470 };
471
472 BitRefLevel bit_sum;
473 for (auto l = 0; l != nb_ref_levels; ++l) {
474 CHKERR refine_mesh(l + 1);
475 CHKERR mark_skins(l, l + 1);
476 CHKERR mark_skins(l + 1, l + 1);
477 bit_sum |= bit(l);
478 }
479 bit_sum |= bit(nb_ref_levels);
480
481 simpleInterface->getBitRefLevel() = bit_sum;
483
485}
486//! [Read mesh]
487
488//! [Set up problem]
491 // Add field
496
497 constexpr int order = 3;
499
500 // Simple interface will resolve adjacency to DOFs of parent of the element.
501 // Using that information MAtrixManager allocate appropriately size of
502 // matrix.
504 BitRefLevel bit_marker;
505 for (auto l = 1; l <= nb_ref_levels; ++l)
506 bit_marker |= marker(l);
507 simpleInterface->getBitAdjEnt() = bit_marker;
508
510
513
514 // remove obsolete DOFs from problem
515
516 for (int l = 0; l != nb_ref_levels; ++l) {
517 CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
518 FIELD_NAME, bit(l), bit(l));
519 CHKERR prb_mng->removeDofsOnEntities(simpleInterface->getProblemName(),
520 FIELD_NAME, marker(l + 1),
521 bit(l).flip());
522 }
524}
525//! [Set up problem]
526
527//! [Push operators to pipeline]
531
532 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
533
534 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
535 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
536
537 pipeline_mng->getDomainLhsFE()->exeTestHook = test_bit_child;
538 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
539
540 auto beta = [](const double, const double, const double) { return 1; };
541 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
542 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
543 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
544 DomainEleOp::OPROW, QUIET, Sev::noisy);
545 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainLhsFE(),
546 DomainEleOp::OPCOL, QUIET, Sev::noisy);
547 pipeline_mng->getOpDomainLhsPipeline().push_back(
549
550 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
552 field_op_row->doWorkRhsHook = [](DataOperator *op_ptr, int side,
556 if (type == MBENTITYSET) {
557
558 MOFEM_LOG("SELF", Sev::verbose)
559 << "ROW: side/type: " << side << "/" << CN::EntityTypeName(type)
560 << " op space/base: " << FieldSpaceNames[data.getSpace()] << "/"
561 << ApproximationBaseNames[data.getBase()] << " DOFs "
562 << data.getIndices() << " nb base functions " << data.getN().size2()
563 << " nb base functions integration points " << data.getN().size1();
564
565 auto get_bool = [](auto fe, auto bit) {
566 return (bit & fe->getBitRefLevel()).any();
567 };
568
569 for (auto &field_ent : data.getFieldEntities()) {
570 MOFEM_LOG("SELF", Sev::verbose)
571 << "\t" << CN::EntityTypeName(field_ent->getEntType());
572 }
573 }
575 };
576
577 set_parent_dofs<DomainParentEle>(
579 Sev::verbose);
580 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
581 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
582 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
583
584 pipeline_mng->getOpDomainRhsPipeline().push_back(
586
588}
589//! [Push operators to pipeline]
590
591//! [Solve]
595
596 MOFEM_LOG("WORLD", Sev::inform) << "Solve problem";
597
598 auto solver = pipeline_mng->createKSP();
599 CHKERR KSPSetFromOptions(solver);
600 CHKERR KSPSetUp(solver);
601
602 auto dm = simpleInterface->getDM();
603 auto D = smartCreateDMVector(dm);
604 auto F = smartVectorDuplicate(D);
605
606 CHKERR KSPSolve(solver, F, D);
607 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
608 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
609 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
611}
612
613//! [Check results]
617 pipeline_mng->getDomainLhsFE().reset();
618 pipeline_mng->getDomainRhsFE().reset();
619 pipeline_mng->getBoundaryRhsFE().reset();
620
621 auto rule = [](int, int, int p) -> int { return 2 * p + 1; };
622 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
623 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(rule);
624 pipeline_mng->getDomainRhsFE()->exeTestHook = test_bit_child;
625 pipeline_mng->getBoundaryRhsFE()->exeTestHook = test_bit_child;
626
627 auto common_data_ptr = boost::make_shared<CommonData>();
628 common_data_ptr->resVec = smartCreateDMVector(simpleInterface->getDM());
629 common_data_ptr->L2Vec = createSmartVectorMPI(
630 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
631 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
632 common_data_ptr->divApproxVals = boost::make_shared<MatrixDouble>();
633
634 auto jac_ptr = boost::make_shared<MatrixDouble>();
635 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
636 auto det_ptr = boost::make_shared<VectorDouble>();
637
638 pipeline_mng->getOpDomainRhsPipeline().push_back(
639 new OpCalculateHOJac<2>(jac_ptr));
640 pipeline_mng->getOpDomainRhsPipeline().push_back(
641 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
642 pipeline_mng->getOpDomainRhsPipeline().push_back(
643 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
644
645 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
646 DomainEleOp::OPSPACE, QUIET, Sev::noisy);
647 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
648 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
649
650 pipeline_mng->getOpDomainRhsPipeline().push_back(
652 FIELD_NAME, common_data_ptr->divApproxVals));
653 pipeline_mng->getOpDomainRhsPipeline().push_back(
655 common_data_ptr->approxVals));
656
657 pipeline_mng->getOpDomainRhsPipeline().push_back(
658 new OpError<FIELD_DIM>(common_data_ptr));
659
660 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
661 BoundaryEleOp::OPSPACE, QUIET, Sev::noisy);
662 set_parent_dofs<BoundaryParentEle>(mField, pipeline_mng->getBoundaryRhsFE(),
663 BoundaryEleOp::OPROW, VERBOSE, Sev::noisy);
664 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
666 common_data_ptr->approxVals));
667 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
668 new OpErrorSkel<FIELD_DIM>(common_data_ptr));
669
670 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
671 CHKERR VecZeroEntries(common_data_ptr->resVec);
672
673 CHKERR pipeline_mng->loopFiniteElements();
674
675 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
676 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
677 CHKERR VecAssemblyBegin(common_data_ptr->resVec);
678 CHKERR VecAssemblyEnd(common_data_ptr->resVec);
679 double nrm2;
680 CHKERR VecNorm(common_data_ptr->resVec, NORM_2, &nrm2);
681 const double *array;
682 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
683 MOFEM_LOG_C("WORLD", Sev::inform, "Error %6.4e Vec norm %6.4e\n",
684 std::sqrt(array[0]), nrm2);
685
686 constexpr double eps = 1e-8;
687 if (nrm2 > eps)
688 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
689 "Not converged solution err = %6.4e", nrm2);
690 if (std::sqrt(array[0]) > eps)
691 SETERRQ1(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
692 "Error in approximation err = %6.4e", std::sqrt(array[0]));
693
694 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
695
697}
698//! [Check results]
699
703 pipeline_mng->getDomainLhsFE().reset();
704 pipeline_mng->getDomainRhsFE().reset();
705
706 auto rule = [](int, int, int p) -> int { return -1; };
707 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
708
709 static_cast<ForcesAndSourcesCore *>(pipeline_mng->getDomainRhsFE().get())
710 ->setRuleHook = [](ForcesAndSourcesCore *fe_raw_ptr, int order_row,
711 int order_col, int order_data) -> MoFEMErrorCode {
713 fe_raw_ptr->gaussPts.resize(3, 3);
714 fe_raw_ptr->gaussPts(0, 0) = 0;
715 fe_raw_ptr->gaussPts(1, 0) = 0;
716 fe_raw_ptr->gaussPts(2, 0) = 0;
717 fe_raw_ptr->gaussPts(0, 1) = 1;
718 fe_raw_ptr->gaussPts(1, 1) = 0;
719 fe_raw_ptr->gaussPts(2, 1) = 0;
720 fe_raw_ptr->gaussPts(0, 2) = 0;
721 fe_raw_ptr->gaussPts(1, 2) = 1;
722 fe_raw_ptr->gaussPts(2, 2) = 0;
724 };
725
726 auto field_op_row = new ForcesAndSourcesCore::UserDataOperator(
728
729 auto approx_vals = boost::make_shared<VectorDouble>();
730
731 auto &moab = mField.get_moab();
732 Tag th;
733 double def_val[] = {0};
734 CHKERR moab.tag_get_handle("FIELD", 1, MB_TYPE_DOUBLE, th,
735 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
736
737 field_op_row->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
741 if (type == MBVERTEX) {
742 auto op_ptr =
744 base_op_ptr);
745 auto t_field = getFTensor0FromVec(*approx_vals);
746 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
747 if (nb_gauss_pts != 3)
748 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
749 "Should be three guass pts.");
750 auto conn = op_ptr->getConn();
751 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
752 const double v = t_field;
753 CHKERR moab.tag_set_data(th, &conn[gg], 1, &v);
754 ++t_field;
755 }
756 }
758 };
759
760 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
761 DomainEleOp::OPSPACE, VERBOSE, Sev::noisy);
762 set_parent_dofs<DomainParentEle>(mField, pipeline_mng->getDomainRhsFE(),
763 DomainEleOp::OPROW, VERBOSE, Sev::noisy);
764 pipeline_mng->getOpDomainRhsPipeline().push_back(
765 new OpCalculateScalarFieldValues(FIELD_NAME, approx_vals));
766 pipeline_mng->getOpDomainRhsPipeline().push_back(field_op_row);
767 CHKERR pipeline_mng->loopFiniteElements();
768
769 CHKERR mField.getInterface<BitRefManager>()->writeBitLevelByType(
770 bit(nb_ref_levels), BitRefLevel().set(), MBTRI, "out.vtk", "VTK", "");
771
773}
774
775int main(int argc, char *argv[]) {
776
777 // Initialisation of MoFEM/PETSc and MOAB data structures
778 MoFEM::Core::Initialize(&argc, &argv, NULL, help);
779
780 try {
781
782 //! [Register MoFEM discrete manager in PETSc]
783 DMType dm_name = "DMMOFEM";
784 CHKERR DMRegister_MoFEM(dm_name);
785 //! [Register MoFEM discrete manager in PETSc
786
787 //! [Create MoAB]
788 moab::Core mb_instance; ///< mesh database
789 moab::Interface &moab = mb_instance; ///< mesh database interface
790 //! [Create MoAB]
791
792 //! [Create MoFEM]
793 MoFEM::Core core(moab); ///< finite element database
794 MoFEM::Interface &m_field = core; ///< finite element database insterface
795 //! [Create MoFEM]
796
797 //! [AtomTest]
798 AtomTest ex(m_field);
799 CHKERR ex.runProblem();
800 //! [AtomTest]
801 }
803
805}
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ QUIET
Definition: definitions.h:221
@ VERBOSE
Definition: definitions.h:222
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
static const char *const FieldSpaceNames[]
Definition: definitions.h:105
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
static const char *const ApproximationBaseNames[]
Definition: definitions.h:85
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
constexpr double beta
constexpr double alpha
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:482
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
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:311
SeverityLevel
Severity levels.
Definition: LogManager.hpp:43
auto marker
set bit to marker
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
evaluate source, i.e. rhs vector
constexpr int SPACE_DIM
constexpr int FIELD_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
evaluate mass matrix
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
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temprary meshset.
Definition: Templates.hpp:1479
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
const double r
rate factor
constexpr double t
plate stiffness
Definition: plate.cpp:76
FTensor::Index< 'm', 3 > m
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:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
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.
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:35
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:272
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:378
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:294
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:810
auto & getBitAdjEnt()
bit ref level for parent
Definition: Simple.hpp:447
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:308
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:593
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:396
BitRefLevel & getBitRefLevelMask()
Get the BitRefLevel.
Definition: Simple.hpp:314
EntityHandle & getBoundaryMeshSet()
Get the BoundaryMeshSet object.
Definition: Simple.hpp:293
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:753
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:342
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition: Simple.hpp:307
bool & getParentAdjacencies()
Get the addParentAdjacencies.
Definition: Simple.hpp:426
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.