v0.13.1
Simple.cpp
Go to the documentation of this file.
1/** \file Simple.cpp
2 * \brief Implementation of simple interface
3 * \ingroup mofem_simple_interface
4 */
5
6
7namespace MoFEM {
8
9MoFEMErrorCode Simple::query_interface(boost::typeindex::type_index type_index,
10 UnknownInterface **iface) const {
11 *iface = const_cast<Simple *>(this);
12 return 0;
13}
14
15template <int DIM>
17 static_assert(DIM == 2 || DIM == 3, "not implemented");
19}
20
21template <>
22MoFEMErrorCode Simple::setSkeletonAdjacency<2>(std::string fe_name) {
23 Interface &m_field = cOre;
25
26 auto defaultSkeletonEdge =
27 [&](moab::Interface &moab, const Field &field, const EntFiniteElement &fe,
28 std::vector<EntityHandle> &adjacency) -> MoFEMErrorCode {
30
31 CHKERR DefaultElementAdjacency::defaultEdge(moab, field, fe, adjacency);
32
33 if (std::find(domainFields.begin(), domainFields.end(), field.getName()) !=
34 domainFields.end()) {
35
36 const EntityHandle fe_ent = fe.getEnt();
37 std::vector<EntityHandle> bride_adjacency_edge;
38 CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, bride_adjacency_edge);
39
40 switch (field.getSpace()) {
41 case H1:
42 CHKERR moab.get_connectivity(&*bride_adjacency_edge.begin(),
43 bride_adjacency_edge.size(), adjacency,
44 true);
45 case HCURL:
46 case HDIV:
47 CHKERR moab.get_adjacencies(&*bride_adjacency_edge.begin(),
48 bride_adjacency_edge.size(), 1, false,
49 adjacency, moab::Interface::UNION);
50 case L2:
51 adjacency.insert(adjacency.end(), bride_adjacency_edge.begin(),
52 bride_adjacency_edge.end());
53 break;
54 case NOFIELD:
55 break;
56 default:
57 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
58 "this field is not implemented for TRI finite element");
59 }
60
61 std::sort(adjacency.begin(), adjacency.end());
62 auto it = std::unique(adjacency.begin(), adjacency.end());
63
64 std::vector<EntityHandle> new_adjacency(
65 std::distance(adjacency.begin(), it));
66 std::copy(adjacency.begin(), it, new_adjacency.begin());
67
68 for (auto e : new_adjacency) {
69 auto side_table = fe.getSideNumberTable();
70 if (side_table.find(e) == side_table.end())
71 const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
72 .insert(
73 boost::shared_ptr<SideNumber>(new SideNumber(e, -1, 0, 0)));
74 }
75
76 adjacency.swap(new_adjacency);
77 }
78
80 };
81
82 CHKERR m_field.modify_finite_element_adjacency_table(fe_name, MBEDGE,
83 defaultSkeletonEdge);
84
86}
87
88template <>
89MoFEMErrorCode Simple::setSkeletonAdjacency<3>(std::string fe_name) {
90 Interface &m_field = cOre;
92
93 auto defaultSkeletonEdge =
94 [&](moab::Interface &moab, const Field &field, const EntFiniteElement &fe,
95 std::vector<EntityHandle> &adjacency) -> MoFEMErrorCode {
97
98 CHKERR DefaultElementAdjacency::defaultFace(moab, field, fe, adjacency);
99
100 if (std::find(domainFields.begin(), domainFields.end(), field.getName()) !=
101 domainFields.end()) {
102
103 const EntityHandle fe_ent = fe.getEnt();
104 std::vector<EntityHandle> bride_adjacency_edge;
105 CHKERR moab.get_adjacencies(&fe_ent, 1, 2, false, bride_adjacency_edge);
106
107 switch (field.getSpace()) {
108 case H1:
109 CHKERR moab.get_connectivity(&*bride_adjacency_edge.begin(),
110 bride_adjacency_edge.size(), adjacency,
111 true);
112 case HCURL:
113 CHKERR moab.get_adjacencies(&*bride_adjacency_edge.begin(),
114 bride_adjacency_edge.size(), 1, false,
115 adjacency, moab::Interface::UNION);
116 case HDIV:
117 CHKERR moab.get_adjacencies(&*bride_adjacency_edge.begin(),
118 bride_adjacency_edge.size(), 2, false,
119 adjacency, moab::Interface::UNION);
120 case L2:
121 adjacency.insert(adjacency.end(), bride_adjacency_edge.begin(),
122 bride_adjacency_edge.end());
123 break;
124 case NOFIELD:
125 break;
126 default:
127 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
128 "this field is not implemented for TRI finite element");
129 }
130
131 std::sort(adjacency.begin(), adjacency.end());
132 auto it = std::unique(adjacency.begin(), adjacency.end());
133
134 std::vector<EntityHandle> new_adjacency(
135 std::distance(adjacency.begin(), it));
136 std::copy(adjacency.begin(), it, new_adjacency.begin());
137
138 for (auto e : new_adjacency) {
139 auto side_table = fe.getSideNumberTable();
140 if (side_table.find(e) == side_table.end())
141 const_cast<SideNumber_multiIndex &>(fe.getSideNumberTable())
142 .insert(
143 boost::shared_ptr<SideNumber>(new SideNumber(e, -1, 0, 0)));
144 }
145
146 adjacency.swap(new_adjacency);
147 }
148
150 };
151
152 CHKERR m_field.modify_finite_element_adjacency_table(fe_name, MBTRI,
153 defaultSkeletonEdge);
154 CHKERR m_field.modify_finite_element_adjacency_table(fe_name, MBQUAD,
155 defaultSkeletonEdge);
156
158}
159
160template <>
163
164 switch (getDim()) {
165 case 1:
166 THROW_MESSAGE("Not implemented");
167 case 2:
168 return setSkeletonAdjacency<2>(fe_name);
169 case 3:
170 return setSkeletonAdjacency<3>(fe_name);
171 default:
172 THROW_MESSAGE("Not implemented");
173 }
175}
176
178 static_assert(DIM == 2 || DIM == 3, "not implemented");
180}
181
182template <> MoFEMErrorCode Simple::setParentAdjacency<3>() {
183 Interface &m_field = cOre;
185
187 boost::make_shared<ParentFiniteElementAdjacencyFunction<3>>(
190 boost::make_shared<ParentFiniteElementAdjacencyFunction<2>>(
192
197 if (addBoundaryFE || !boundaryFields.empty()) {
202 }
203 if (addSkeletonFE || !skeletonFields.empty()) {
208 }
209
211}
212
213template <> MoFEMErrorCode Simple::setParentAdjacency<2>() {
214 Interface &m_field = cOre;
216
218 boost::make_shared<ParentFiniteElementAdjacencyFunction<2>>(
221 boost::make_shared<ParentFiniteElementAdjacencyFunction<1>>(
223
228 if (addBoundaryFE || !boundaryFields.empty())
231 if (addSkeletonFE || !skeletonFields.empty())
234
236}
237
240 switch (getDim()) {
241 case 1:
242 THROW_MESSAGE("Not implemented");
243 case 2:
244 return setParentAdjacency<2>();
245 case 3:
246 return setParentAdjacency<3>();
247 default:
248 THROW_MESSAGE("Not implemented");
249 }
251}
252
255 if (dim == -1)
256 dim = getDim();
257
258 if (fe_name.empty())
259 fe_name = skeletonFE;
260
261 switch (dim) {
262 case 2:
263 return setSkeletonAdjacency<2>(fe_name);
264 case 3:
265 return setSkeletonAdjacency<3>(fe_name);
266 default:
267 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Not implemented");
268 }
270}
271
273 : cOre(const_cast<Core &>(core)), bitLevel(BitRefLevel().set(0)),
274 bitLevelMask(BitRefLevel().set()), meshSet(0), boundaryMeshset(0),
275 skeletonMeshset(0), nameOfProblem("SimpleProblem"), domainFE("dFE"),
276 boundaryFE("bFE"), skeletonFE("sFE"), dIm(-1), addSkeletonFE(false),
277 addBoundaryFE(false), addParentAdjacencies(false),
278 bitAdjParent(BitRefLevel().set()), bitAdjParentMask(BitRefLevel().set()),
279 bitAdjEnt(BitRefLevel().set()), bitAdjEntMask(BitRefLevel().set()) {
280 PetscLogEventRegister("SimpleSetUp", 0, &MOFEM_EVENT_SimpleSetUP);
281 PetscLogEventRegister("SimpleLoadMesh", 0, &MOFEM_EVENT_SimpleLoadMesh);
282 PetscLogEventRegister("SimpleBuildFields", 0, &MOFEM_EVENT_SimpleBuildFields);
283 PetscLogEventRegister("SimpleBuildFiniteElements", 0,
285 PetscLogEventRegister("SimpleSetUp", 0, &MOFEM_EVENT_SimpleBuildProblem);
286 PetscLogEventRegister("SimpleKSPSolve", 0, &MOFEM_EVENT_SimpleKSPSolve);
287 strcpy(meshFileName, "mesh.h5m");
288}
289
291 PetscBool flg = PETSC_TRUE;
293 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Simple interface options",
294 "none");
295 CHKERRG(ierr);
296 ierr = PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
297 meshFileName, 255, &flg);
298 CHKERRG(ierr);
299 ierr = PetscOptionsEnd();
300 CHKERRG(ierr);
302}
303
304MoFEMErrorCode Simple::loadFile(const std::string options,
305 const std::string mesh_file_name) {
306 Interface &m_field = cOre;
308 PetscLogEventBegin(MOFEM_EVENT_SimpleLoadMesh, 0, 0, 0, 0);
309
310 if (!mesh_file_name.empty())
311 strcpy(meshFileName, mesh_file_name.c_str());
312
313 // This is a case of distributed mesh and algebra. In that case each processor
314 // keep only part of the problem.
315 CHKERR m_field.get_moab().load_file(meshFileName, 0, options.c_str());
316 CHKERR m_field.rebuild_database();
317
318 // determine problem dimension
319 if (dIm == -1) {
320 int nb_ents_3d;
321 CHKERR m_field.get_moab().get_number_entities_by_dimension(
322 meshSet, 3, nb_ents_3d, true);
323 if (nb_ents_3d > 0) {
324 dIm = 3;
325 } else {
326 int nb_ents_2d;
327 CHKERR m_field.get_moab().get_number_entities_by_dimension(
328 meshSet, 2, nb_ents_2d, true);
329 if (nb_ents_2d > 0) {
330 dIm = 2;
331 } else {
332 dIm = 1;
333 }
334 }
335 }
336
337 if (!boundaryMeshset)
339 if (!skeletonMeshset)
341 if (addSkeletonFE)
343
344 if (bitLevel.any()) {
345 Range ents;
346 CHKERR m_field.get_moab().get_entities_by_dimension(meshSet, dIm, ents,
347 true);
348 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents, bitLevel,
349 false);
350 } else {
351 MOFEM_LOG("WORLD", Sev::warning) << "BitRefLevel is none and not set";
352 CHKERR m_field.getInterface<BitRefManager>()->addToDatabaseBitRefLevelByDim(
353 dIm, BitRefLevel().set(), BitRefLevel().set());
354 }
355
356 PetscLogEventEnd(MOFEM_EVENT_SimpleLoadMesh, 0, 0, 0, 0);
358}
359
362 Interface &m_field = cOre;
363 if (m_field.get_comm_size() == 1)
365 else
366 CHKERR loadFile("PARALLEL=READ_PART;"
367 "PARALLEL_RESOLVE_SHARED_ENTS;"
368 "PARTITION=PARALLEL_PARTITION;",
371}
372
374Simple::addDomainField(const std::string &name, const FieldSpace space,
375 const FieldApproximationBase base,
376 const FieldCoefficientsNumber nb_of_cooficients,
377 const TagType tag_type, const enum MoFEMTypes bh,
378 int verb) {
379
380 Interface &m_field = cOre;
382 CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
383 verb);
384 if (space == NOFIELD)
385 noFieldFields.push_back(name);
386 else
387 domainFields.push_back(name);
389}
390
392Simple::addBoundaryField(const std::string &name, const FieldSpace space,
393 const FieldApproximationBase base,
394 const FieldCoefficientsNumber nb_of_cooficients,
395 const TagType tag_type, const enum MoFEMTypes bh,
396 int verb) {
397 Interface &m_field = cOre;
399 CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
400 verb);
401 boundaryFields.push_back(name);
402 if (space == NOFIELD)
403 SETERRQ(
404 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
405 "NOFIELD space for boundary filed not implemented in Simple interface");
407}
408
410Simple::addSkeletonField(const std::string &name, const FieldSpace space,
411 const FieldApproximationBase base,
412 const FieldCoefficientsNumber nb_of_cooficients,
413 const TagType tag_type, const enum MoFEMTypes bh,
414 int verb) {
415
416 Interface &m_field = cOre;
418 CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
419 verb);
420 skeletonFields.push_back(name);
421 if (space == NOFIELD)
422 SETERRQ(
423 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
424 "NOFIELD space for boundary filed not implemented in Simple interface");
425
427}
428
430Simple::addDataField(const std::string &name, const FieldSpace space,
431 const FieldApproximationBase base,
432 const FieldCoefficientsNumber nb_of_cooficients,
433 const TagType tag_type, const enum MoFEMTypes bh,
434 int verb) {
435
436 Interface &m_field = cOre;
438 CHKERR m_field.add_field(name, space, base, nb_of_cooficients, tag_type, bh,
439 verb);
440 if (space == NOFIELD)
441 noFieldDataFields.push_back(name);
442 else
443 dataFields.push_back(name);
445}
446
448 Interface &m_field = cOre;
450
451 auto remove_field_from_list = [&](auto &vec) {
452 auto it = std::find(vec.begin(), vec.end(), name);
453 if (it != vec.end())
454 vec.erase(it);
455 };
456
457 remove_field_from_list(noFieldFields);
458 remove_field_from_list(domainFields);
459
461}
462
464 Interface &m_field = cOre;
466
467 auto remove_field_from_list = [&](auto &vec) {
468 auto it = std::find(vec.begin(), vec.end(), name);
469 if (it != vec.end())
470 vec.erase(it);
471 };
472
473 remove_field_from_list(boundaryFields);
474
476}
477
479 Interface &m_field = cOre;
481
482 auto remove_field_from_list = [&](auto &vec) {
483 auto it = std::find(vec.begin(), vec.end(), name);
484 if (it != vec.end())
485 vec.erase(it);
486 };
487
488 remove_field_from_list(skeletonFields);
489
491}
492
494 Interface &m_field = cOre;
496
497 auto clear_rows_and_cols = [&](auto &fe_name) {
499 auto fe_ptr = m_field.get_finite_elements();
500 auto &fe_by_name = const_cast<FiniteElement_multiIndex *>(fe_ptr)
501 ->get<FiniteElement_name_mi_tag>();
502 auto it_fe = fe_by_name.find(fe_name);
503 if (it_fe != fe_by_name.end()) {
504
505 if (!fe_by_name.modify(it_fe, FiniteElement_row_change_bit_reset()))
506 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
507 "modification unsuccessful");
508
509 if (!fe_by_name.modify(it_fe, FiniteElement_col_change_bit_reset()))
510 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
511 "modification unsuccessful");
512 }
514 };
515 CHKERR clear_rows_and_cols(domainFE);
516 CHKERR clear_rows_and_cols(boundaryFE);
517 CHKERR clear_rows_and_cols(skeletonFE);
518
519 // Define finite elements
521
522 auto add_fields = [&](auto &fe_name, auto &fields) {
524 for (auto &field : fields) {
525 CHKERR m_field.modify_finite_element_add_field_row(fe_name, field);
526 CHKERR m_field.modify_finite_element_add_field_col(fe_name, field);
527 CHKERR m_field.modify_finite_element_add_field_data(fe_name, field);
528 }
530 };
531
532 auto add_data_fields = [&](auto &fe_name, auto &fields) {
534 for (auto &field : fields)
535 CHKERR m_field.modify_finite_element_add_field_data(fe_name, field);
537 };
538
539 CHKERR add_fields(domainFE, domainFields);
540 CHKERR add_data_fields(domainFE, dataFields);
541 CHKERR add_fields(domainFE, noFieldFields);
542 CHKERR add_data_fields(domainFE, noFieldDataFields);
543
544 if (addBoundaryFE || !boundaryFields.empty()) {
546 CHKERR add_fields(boundaryFE, domainFields);
547 if (!boundaryFields.empty())
548 CHKERR add_fields(boundaryFE, boundaryFields);
549 CHKERR add_data_fields(boundaryFE, dataFields);
550 CHKERR add_data_fields(boundaryFE, noFieldDataFields);
551 CHKERR add_fields(boundaryFE, noFieldFields);
552 }
553 if (addSkeletonFE || !skeletonFields.empty()) {
555 CHKERR add_fields(skeletonFE, domainFields);
556 if (!skeletonFields.empty())
557 CHKERR add_fields(skeletonFE, skeletonFields);
558 CHKERR add_data_fields(skeletonFE, dataFields);
559 CHKERR add_data_fields(skeletonFE, noFieldDataFields);
560 CHKERR add_fields(skeletonFE, noFieldFields);
561 }
563}
564
565MoFEMErrorCode Simple::defineProblem(const PetscBool is_partitioned) {
566 Interface &m_field = cOre;
568 // Create dm instance
569 dM = createSmartDM(m_field.get_comm(), "DMMOFEM");
570 // set dm data structure which created mofem data structures
573 CHKERR DMSetFromOptions(dM);
575 if (addBoundaryFE || !boundaryFields.empty()) {
577 }
578 if (addSkeletonFE || !skeletonFields.empty()) {
580 }
581 for (std::vector<std::string>::iterator fit = otherFEs.begin();
582 fit != otherFEs.end(); ++fit) {
583 CHKERR DMMoFEMAddElement(dM, fit->c_str());
584 }
585 CHKERR DMMoFEMSetIsPartitioned(dM, is_partitioned);
587}
588
590 const int order, const Range *ents) {
592 fieldsOrder.insert(
593 {field_name, {order, ents == NULL ? Range() : Range(*ents)}});
595}
596
598 Interface &m_field = cOre;
600 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildFields, 0, 0, 0, 0);
601
602 auto comm_interface_ptr = m_field.getInterface<CommInterface>();
603 auto bit_ref_ptr = m_field.getInterface<BitRefManager>();
604
605 // Add entities to the fields
606 auto add_ents_to_field = [&](auto meshset, auto dim, auto &fields) {
608 for (auto &field : fields) {
609 CHKERR m_field.add_ents_to_field_by_dim(meshset, dim, field);
610 CHKERR comm_interface_ptr->synchroniseFieldEntities(field, 0);
611 }
613 };
614
615 auto make_no_field_ents = [&](auto &fields) {
617 for (auto &field : fields) {
618 std::array<double, 6> coords = {0, 0, 0, 0, 0, 0};
619 CHKERR m_field.create_vertices_and_add_to_field(field, coords.data(), 1);
620 CHKERR comm_interface_ptr->makeFieldEntitiesMultishared(field, 0);
621 CHKERR bit_ref_ptr->setFieldEntitiesBitRefLevel(field, bitLevel);
622 }
624 };
625
626 CHKERR add_ents_to_field(meshSet, dIm, domainFields);
627 CHKERR add_ents_to_field(meshSet, dIm, dataFields);
628 CHKERR add_ents_to_field(boundaryMeshset, dIm - 1, boundaryFields);
629 CHKERR add_ents_to_field(skeletonMeshset, dIm - 1, skeletonFields);
630
631 std::set<std::string> nofield_fields;
632 for (auto &f : noFieldFields)
633 nofield_fields.insert(f);
634 for (auto &f : noFieldDataFields)
635 nofield_fields.insert(f);
636
637 CHKERR make_no_field_ents(nofield_fields);
638
639 // Set order
640 auto set_order = [&](auto meshset, auto dim, auto &fields) {
642 for (auto &field : fields) {
643 if (fieldsOrder.find(field) == fieldsOrder.end()) {
644 SETERRQ1(PETSC_COMM_WORLD, MOFEM_INVALID_DATA,
645 "Order for field not set %s", field.c_str());
646 }
647 int dds = 0;
648 const Field *field_ptr = m_field.get_field_structure(field);
649 switch (field_ptr->getSpace()) {
650 case L2:
651 dds = dim;
652 break;
653 case HDIV:
654 dds = 2;
655 break;
656 case HCURL:
657 dds = 1;
658 break;
659 case H1:
660 dds = 1;
661 break;
662 default:
663 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
664 "Glasgow we have a problem");
665 }
666
667 auto set_order = [&](auto field, auto &ents) {
669
670 auto range = fieldsOrder.equal_range(field);
671 for (auto o = range.first; o != range.second; ++o) {
672 if (!o->second.second.empty())
673 ents = intersect(ents, o->second.second);
674 CHKERR m_field.set_field_order(ents, field, o->second.first);
675 }
676
678 };
679
680 if (field_ptr->getSpace() == H1) {
682 Range ents;
683 CHKERR m_field.get_field_entities_by_dimension(field, 0, ents);
684 CHKERR set_order(field, ents);
685 } else {
686 CHKERR m_field.set_field_order(meshSet, MBVERTEX, field, 1);
687 }
688 }
689 for (int dd = dds; dd <= dim; dd++) {
690 Range ents;
691 CHKERR m_field.get_field_entities_by_dimension(field, dd, ents);
692 CHKERR set_order(field, ents);
693 }
694 }
696 };
697
698 CHKERR set_order(meshSet, dIm, domainFields);
699 CHKERR set_order(meshSet, dIm, dataFields);
700 CHKERR set_order(meshSet, dIm - 1, boundaryFields);
701 CHKERR set_order(meshSet, dIm - 1, skeletonFields);
702
703 // Build fields
704 CHKERR m_field.build_fields();
705 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildFields, 0, 0, 0, 0);
707}
708
710 Interface &m_field = cOre;
712 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildFiniteElements, 0, 0, 0, 0);
713 // Add finite elements
715 true);
717 if (addBoundaryFE || !boundaryFields.empty()) {
719 boundaryFE, true);
721 }
722 if (addSkeletonFE || !skeletonFields.empty()) {
724 skeletonFE, true);
726 }
727 for (std::vector<std::string>::iterator fit = otherFEs.begin();
728 fit != otherFEs.end(); ++fit) {
729 CHKERR m_field.build_finite_elements(*fit);
730 }
731 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildFiniteElements, 0, 0, 0, 0);
733}
734
736 Interface &m_field = cOre;
738 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
740 // Set problem by the DOFs on the fields rather that by adding DOFs on the
741 // elements
742 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
744 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
745 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
747}
748
749MoFEMErrorCode Simple::setUp(const PetscBool is_partitioned) {
750 Interface &m_field = cOre;
752
753 PetscLogEventBegin(MOFEM_EVENT_SimpleSetUP, 0, 0, 0, 0);
754
756
757 if (addSkeletonFE || !skeletonFields.empty())
759
762
763 CHKERR defineProblem(is_partitioned);
767
768 PetscLogEventEnd(MOFEM_EVENT_SimpleSetUP, 0, 0, 0, 0);
770}
771
773 Interface &m_field = cOre;
775 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
776
777 if (!only_dm) {
780 if (addSkeletonFE || !skeletonFields.empty())
785 }
786
788
789 const Problem *problem_ptr;
790 CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
791 const auto problem_name = problem_ptr->getName();
792 CHKERR m_field.modify_problem_ref_level_set_bit(problem_name, bitLevel);
795
796 // Set problem by the DOFs on the fields rather that by adding DOFs on the
797 // elements
798 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
800 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
801
802 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
804}
805
808 CHKERR PetscObjectReference(getPetscObject(dM.get()));
809 *dm = dM.get();
811}
812
813/**
814 * @brief Delete dm
815 *
816 * @return MoFEMErrorCode
817 */
820 dM.reset();
822}
823
824/**
825 * @brief Delete finite elements
826 *
827 * @return MoFEMErrorCode
828 */
830 Interface &m_field = cOre;
832 for (auto fe : {domainFE, boundaryFE, skeletonFE}) {
833 if (m_field.check_finite_element(fe)) {
834 CHKERR m_field.delete_finite_element(fe);
835 }
836 }
838}
839
841Simple::addFieldToEmptyFieldBlocks(const std::string row_field,
842 const std::string col_field) const {
843 Interface &m_field = cOre;
846 getProblemName(), row_field, col_field);
848}
849
851 Interface &m_field = cOre;
853 ParallelComm *pcomm =
854 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
855
856 auto get_skin = [&](auto meshset) {
857 // filter not owned entities, those are not on boundary
858
859 Range domain_ents;
860 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, dIm,
861 domain_ents, true);
862 CHKERR pcomm->filter_pstatus(domain_ents,
863 PSTATUS_SHARED | PSTATUS_MULTISHARED,
864 PSTATUS_NOT, -1, nullptr);
865
866 Skinner skin(&m_field.get_moab());
867 Range domain_skin;
868 CHKERR skin.find_skin(0, domain_ents, false, domain_skin);
869 CHKERR pcomm->filter_pstatus(domain_skin,
870 PSTATUS_SHARED | PSTATUS_MULTISHARED,
871 PSTATUS_NOT, -1, nullptr);
872 return domain_skin;
873 };
874
875 auto create_boundary_meshset = [&](auto &&domain_skin) {
877 // create boundary meshset
878 if (boundaryMeshset != 0) {
880 }
881 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, boundaryMeshset);
882 CHKERR m_field.get_moab().add_entities(boundaryMeshset, domain_skin);
883 for (int dd = 0; dd != dIm - 1; dd++) {
884 Range adj;
885 CHKERR m_field.get_moab().get_adjacencies(domain_skin, dd, false, adj,
886 moab::Interface::UNION);
887 CHKERR m_field.get_moab().add_entities(boundaryMeshset, adj);
888 }
890 };
891
892 CHKERR create_boundary_meshset(get_skin(meshSet));
893
895}
896
898 Interface &m_field = cOre;
900
901 ParallelComm *pcomm =
902 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
903
904 auto create_skeleton_meshset = [&](auto meshset) {
906 // create boundary meshset
907 if (skeletonMeshset != 0) {
909 }
910 Range boundary_ents, skeleton_ents;
911 CHKERR m_field.get_moab().get_entities_by_dimension(boundaryMeshset,
912 dIm - 1, boundary_ents);
913 Range domain_ents;
914 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, dIm,
915 domain_ents, true);
916 CHKERR m_field.get_moab().get_adjacencies(
917 domain_ents, dIm - 1, false, skeleton_ents, moab::Interface::UNION);
918 skeleton_ents = subtract(skeleton_ents, boundary_ents);
919 CHKERR pcomm->filter_pstatus(skeleton_ents, PSTATUS_NOT_OWNED, PSTATUS_NOT,
920 -1, nullptr);
921 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, skeletonMeshset);
922 CHKERR m_field.get_moab().add_entities(skeletonMeshset, skeleton_ents);
924 };
925
926 CHKERR create_skeleton_meshset(meshSet);
927
929}
930
932 Interface &m_field = cOre;
934 MOFEM_LOG("WORLD", Sev::verbose) << "Exchange ghost cells";
935
936 ParallelComm *pcomm =
937 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
938 if (pcomm == NULL)
939 pcomm = new ParallelComm(&m_field.get_moab(), m_field.get_comm());
940
941 Range verts;
942 CHKERR m_field.get_moab().get_entities_by_type(0, MBVERTEX, verts);
943 CHKERR pcomm->exchange_ghost_cells(getDim(), getDim() - 1, 1,
944 3 /**get all adjacent ghosted entities */,
945 true, false, meshSet ? &meshSet : nullptr);
946
947 Range shared;
948 CHKERR m_field.get_moab().get_entities_by_dimension(0, dIm, shared);
949 for (auto d = dIm - 1; d >= 1; --d) {
950 CHKERR m_field.get_moab().get_adjacencies(shared, d, false, shared,
951 moab::Interface::UNION);
952 }
953 CHKERR pcomm->filter_pstatus(shared, PSTATUS_SHARED | PSTATUS_MULTISHARED,
954 PSTATUS_OR, -1, &shared);
955 Tag part_tag = pcomm->part_tag();
956 CHKERR pcomm->exchange_tags(part_tag, shared);
957
959}
960
961} // namespace MoFEM
static Index< 'o', 3 > o
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:97
@ MF_ZERO
Definition: definitions.h:98
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
const int dim
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:1070
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:118
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:450
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition: DMMMoFEM.cpp:1228
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
multi_index_container< boost::shared_ptr< FiniteElement >, indexed_by< hashed_unique< tag< FiniteElement_Meshset_mi_tag >, member< FiniteElement, EntityHandle, &FiniteElement::meshset > >, hashed_unique< tag< BitFEId_mi_tag >, const_mem_fun< FiniteElement, BitFEId, &FiniteElement::getId >, HashBit< BitFEId >, EqBit< BitFEId > >, ordered_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< FiniteElement, boost::string_ref, &FiniteElement::getNameRef > > > > FiniteElement_multiIndex
MultiIndex for entities for FiniteElement.
virtual const FiniteElement_multiIndex * get_finite_elements() const =0
Get the finite elements object.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual Field * get_field_structure(const std::string &name)=0
get field structure
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode get_field_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the field by dimension
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
virtual MoFEMErrorCode modify_problem_mask_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set dof mask ref level for problem
virtual MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)=0
set ref level for problem
char mesh_file_name[255]
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
int FieldCoefficientsNumber
Number of field coefficients.
Definition: Types.hpp:27
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: MoFEM.hpp:24
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
PetscObject getPetscObject(T obj)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
constexpr auto field_name
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MPI_Comm & get_comm() const =0
virtual MoFEMErrorCode delete_finite_element(const std::string name, int verb=DEFAULT_VERBOSITY)=0
delete finite element from mofem database
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MoFEMErrorCode rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
virtual MoFEMErrorCode create_vertices_and_add_to_field(const std::string name, const double coords[], int size, int verb=DEFAULT_VERBOSITY)=0
Create a vertices and add to field object.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode defaultEdge(Interface &moab, const Field &field, const EntFiniteElement &fe, std::vector< EntityHandle > &adjacency)
static MoFEMErrorCode defaultFace(Interface &moab, const Field &field, const EntFiniteElement &fe, std::vector< EntityHandle > &adjacency)
Deprecated interface functions.
Finite element data for entity.
Provide data structure for (tensor) field approximation.
FieldApproximationBase getApproxBase() const
Get approximation base.
FieldSpace getSpace() const
Get field approximation space.
keeps basic data about problem
Problem manager is used to build and partition problems.
keeps information about side number for the finite element
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
boost::shared_ptr< ParentFiniteElementAdjacencyFunction< 2 > > parentAdjFunctionDim2
Definition: Simple.hpp:524
EntityHandle boundaryMeshset
meshset with boundary
Definition: Simple.hpp:477
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:735
std::vector< std::string > dataFields
Data fields.
Definition: Simple.hpp:492
MoFEMErrorCode setParentAdjacency()
Definition: Simple.cpp:177
std::vector< std::string > domainFields
domain fields
Definition: Simple.hpp:489
MoFEMErrorCode removeDomainField(const std::string &name)
Remove field form domain.
Definition: Simple.cpp:447
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:493
MoFEM::Core & cOre
Definition: Simple.hpp:464
BitRefLevel bitAdjEnt
bit ref level for parent
Definition: Simple.hpp:486
char meshFileName[255]
mesh file name
Definition: Simple.hpp:506
EntityHandle meshSet
domain meshset
Definition: Simple.hpp:476
int getDim() const
Get the problem dimension.
Definition: Simple.hpp:264
std::vector< std::string > otherFEs
Other finite elements.
Definition: Simple.hpp:504
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:709
MoFEMErrorCode addDataField(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:430
MoFEMErrorCode exchangeGhostCells()
Definition: Simple.cpp:931
std::multimap< std::string, std::pair< int, Range > > fieldsOrder
fields order
Definition: Simple.hpp:497
PetscLogEvent MOFEM_EVENT_SimpleBuildFields
Definition: Simple.hpp:471
PetscLogEvent MOFEM_EVENT_SimpleLoadMesh
Definition: Simple.hpp:470
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:374
MoFEMErrorCode reSetUp(bool only_dm=false)
Rebuild internal MoFEM data structures.
Definition: Simple.cpp:772
std::vector< std::string > noFieldFields
NOFIELD field name.
Definition: Simple.hpp:493
bool addSkeletonFE
Add skeleton FE.
Definition: Simple.hpp:480
BitRefLevel bitAdjParent
bit ref level for parent
Definition: Simple.hpp:484
bool addBoundaryFE
Add boundary FE.
Definition: Simple.hpp:481
MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string row_field, const std::string col_field) const
add empty block to problem
Definition: Simple.cpp:841
std::string skeletonFE
skeleton finite element
Definition: Simple.hpp:502
std::vector< std::string > noFieldDataFields
NOFIELD field name.
Definition: Simple.hpp:494
std::string boundaryFE
boundary finite element
Definition: Simple.hpp:501
std::string nameOfProblem
problem name
Definition: Simple.hpp:499
boost::shared_ptr< ParentFiniteElementAdjacencyFunction< 3 > > parentAdjFunctionDim3
Definition: Simple.hpp:522
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
std::string domainFE
domain finite element
Definition: Simple.hpp:500
std::vector< std::string > boundaryFields
boundary fields
Definition: Simple.hpp:490
BitRefLevel bitAdjEntMask
bit ref level for parent marent
Definition: Simple.hpp:487
SmartPetscObj< DM > dM
Discrete manager (interface to PETSc/MoFEM functions)
Definition: Simple.hpp:510
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:597
SmartPetscObj< DM > getDM()
Return smart DM object.
Definition: Simple.hpp:254
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
BitRefLevel bitLevel
BitRefLevel of the probelm.
Definition: Simple.hpp:466
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
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:392
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:565
std::vector< std::string > skeletonFields
fields on the skeleton
Definition: Simple.hpp:491
MoFEMErrorCode addSkeletonField(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 skeleton.
Definition: Simple.cpp:410
MoFEMErrorCode removeSkeletonField(const std::string &name)
Remove field form skeleton.
Definition: Simple.cpp:478
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: Simple.cpp:9
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode removeBoundaryField(const std::string &name)
Remove field form boundary.
Definition: Simple.cpp:463
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
int dIm
dimension of problem
Definition: Simple.hpp:507
MoFEMErrorCode setSkeletonAdjacency(int dim=-1, std::string fe_name="")
Set the skeleton adjacency object.
Definition: Simple.cpp:253
MoFEMErrorCode deleteDM()
Delete dm.
Definition: Simple.cpp:818
PetscLogEvent MOFEM_EVENT_SimpleSetUP
Definition: Simple.hpp:469
PetscLogEvent MOFEM_EVENT_SimpleKSPSolve
Definition: Simple.hpp:474
MoFEMErrorCode createBoundaryMeshset()
Definition: Simple.cpp:850
PetscLogEvent MOFEM_EVENT_SimpleBuildFiniteElements
Definition: Simple.hpp:472
BitRefLevel bitLevelMask
BitRefLevel of the probelm.
Definition: Simple.hpp:467
BitRefLevel bitAdjParentMask
bit ref level for parent marent
Definition: Simple.hpp:485
PetscLogEvent MOFEM_EVENT_SimpleBuildProblem
Definition: Simple.hpp:473
boost::shared_ptr< ParentFiniteElementAdjacencyFunction< 1 > > parentAdjFunctionDim1
Definition: Simple.hpp:526
MoFEMErrorCode deleteFiniteElements()
Delete finite elements.
Definition: Simple.cpp:829
Simple(const MoFEM::Core &core)
Definition: Simple.cpp:272
EntityHandle skeletonMeshset
skeleton meshset with boundary
Definition: Simple.hpp:478
MoFEMErrorCode createSkeletonMeshset()
Definition: Simple.cpp:897
bool addParentAdjacencies
If set to true parent adjacencies are build.
Definition: Simple.hpp:482
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.