v0.15.0
Loading...
Searching...
No Matches
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
6namespace MoFEM {
7
8MoFEMErrorCode Simple::query_interface(boost::typeindex::type_index type_index,
9 UnknownInterface **iface) const {
10 *iface = const_cast<Simple *>(this);
11 return 0;
12}
13
14template <int DIM>
16 static_assert(DIM == 2 || DIM == 3, "not implemented");
18}
19
20template <>
22 Interface &m_field = cOre;
24
26 boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<1>>(
29 fe_name, MBEDGE, *parentAdjSkeletonFunctionDim1);
30
32}
33
34template <>
36 Interface &m_field = cOre;
38
40 boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
42
44 fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
46 fe_name, MBQUAD, *parentAdjSkeletonFunctionDim2);
47
49}
50
51template <>
54
55 switch (getDim()) {
56 case 1:
57 THROW_MESSAGE("Not implemented");
58 case 2:
59 return setSkeletonAdjacency<2>(fe_name);
60 case 3:
61 return setSkeletonAdjacency<3>(fe_name);
62 default:
63 THROW_MESSAGE("Not implemented");
64 }
66}
68 static_assert(DIM == 2 || DIM == 3, "not implemented");
70}
71
102
104 Interface &m_field = cOre;
106
108 boost::make_shared<ParentFiniteElementAdjacencyFunction<2>>(
111 boost::make_shared<ParentFiniteElementAdjacencyFunction<1>>(
113
118 if (addBoundaryFE || !boundaryFields.empty())
121 if (addSkeletonFE || !skeletonFields.empty())
124
126}
127
130 switch (getDim()) {
131 case 1:
132 THROW_MESSAGE("Not implemented");
133 case 2:
134 return setParentAdjacency<2>();
135 case 3:
136 return setParentAdjacency<3>();
137 default:
138 THROW_MESSAGE("Not implemented");
139 }
141}
142
143MoFEMErrorCode Simple::setSkeletonAdjacency(int dim, std::string fe_name) {
145 if (dim == -1)
146 dim = getDim();
147
148 if (fe_name.empty())
149 fe_name = skeletonFE;
150
151 switch (dim) {
152 case 2:
153 return setSkeletonAdjacency<2>(fe_name);
154 case 3:
155 return setSkeletonAdjacency<3>(fe_name);
156 default:
157 SETERRQ(PETSC_COMM_WORLD, MOFEM_NOT_IMPLEMENTED, "Not implemented");
158 }
160}
161
163 : cOre(const_cast<Core &>(core)), bitLevel(BitRefLevel().set(0)),
164 bitLevelMask(BitRefLevel().set()), meshSet(0), boundaryMeshset(0),
165 skeletonMeshset(0), nameOfProblem("SimpleProblem"), domainFE("dFE"),
166 boundaryFE("bFE"), skeletonFE("sFE"), dIm(-1), addSkeletonFE(false),
167 addBoundaryFE(false), addParentAdjacencies(false),
168 bitAdjParent(BitRefLevel().set()), bitAdjParentMask(BitRefLevel().set()),
169 bitAdjEnt(BitRefLevel().set()), bitAdjEntMask(BitRefLevel().set()) {
170 PetscLogEventRegister("SimpSetUp", 0, &MOFEM_EVENT_SimpleSetUP);
171 PetscLogEventRegister("SimpLoadMesh", 0, &MOFEM_EVENT_SimpleLoadMesh);
172 PetscLogEventRegister("SimpBuildFields", 0, &MOFEM_EVENT_SimpleBuildFields);
173 PetscLogEventRegister("SimpBuildFEs", 0,
175 PetscLogEventRegister("SimpSetUp", 0, &MOFEM_EVENT_SimpleBuildProblem);
176 PetscLogEventRegister("SimpKSPSolve", 0, &MOFEM_EVENT_SimpleKSPSolve);
177 strcpy(meshFileName, "mesh.h5m");
178}
179
181 PetscBool flg = PETSC_TRUE;
183 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Simple interface options",
184 "none");
185 ierr = PetscOptionsString("-file_name", "file name", "", "mesh.h5m",
186 meshFileName, 255, &flg);
187 PetscOptionsEnd();
189}
190
191MoFEMErrorCode Simple::loadFile(const std::string options,
192 const std::string mesh_file_name,
193 LoadFileFunc loadFunc) {
194 Interface &m_field = cOre;
196 PetscLogEventBegin(MOFEM_EVENT_SimpleLoadMesh, 0, 0, 0, 0);
197
198 if (!mesh_file_name.empty())
199 strcpy(meshFileName, mesh_file_name.c_str());
200
201 // Invoke the boost function, passing in required arguments
202 CHKERR loadFunc(m_field, meshFileName, options.c_str());
203 CHKERR m_field.rebuild_database();
204
205 // determine problem dimension
206 if (dIm == -1) {
207 int nb_ents_3d;
208 CHKERR m_field.get_moab().get_number_entities_by_dimension(
209 meshSet, 3, nb_ents_3d, true);
210 if (nb_ents_3d > 0) {
211 dIm = 3;
212 } else {
213 int nb_ents_2d;
214 CHKERR m_field.get_moab().get_number_entities_by_dimension(
215 meshSet, 2, nb_ents_2d, true);
216 if (nb_ents_2d > 0) {
217 dIm = 2;
218 } else {
219 dIm = 1;
220 }
221 }
222 }
223
224 if (!boundaryMeshset)
226 if (!skeletonMeshset)
228 if (addSkeletonFE)
230
231 if (bitLevel.any()) {
232 Range ents;
233 CHKERR m_field.get_moab().get_entities_by_dimension(meshSet, dIm, ents,
234 true);
235 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(ents, bitLevel,
236 false);
237 } else {
238 MOFEM_LOG("WORLD", Sev::warning) << "BitRefLevel is none and not set";
239 CHKERR m_field.getInterface<BitRefManager>()->addToDatabaseBitRefLevelByDim(
240 dIm, BitRefLevel().set(), BitRefLevel().set());
241 }
242
243 PetscLogEventEnd(MOFEM_EVENT_SimpleLoadMesh, 0, 0, 0, 0);
245}
246
249 Interface &m_field = cOre;
250 if (m_field.get_comm_size() == 1)
252 else
253 CHKERR loadFile("PARALLEL=READ_PART;"
254 "PARALLEL_RESOLVE_SHARED_ENTS;"
255 "PARTITION=PARALLEL_PARTITION;",
258}
259
261Simple::addDomainField(const std::string name, const FieldSpace space,
262 const FieldApproximationBase base,
263 const FieldCoefficientsNumber nb_of_coefficients,
264 const TagType tag_type, const enum MoFEMTypes bh,
265 int verb) {
266
267 Interface &m_field = cOre;
269 CHKERR m_field.add_field(name, space, base, nb_of_coefficients, tag_type, bh,
270 verb);
271 domainFields.push_back(name);
272 if (space == NOFIELD)
273 SETERRQ(
274 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
275 "NOFIELD space for boundary field not implemented in Simple interface");
276
277
279}
280
282Simple::addDomainBrokenField(const std::string name, const FieldSpace space,
283 const FieldApproximationBase base,
284 const FieldCoefficientsNumber nb_of_coefficients,
285 const TagType tag_type, const enum MoFEMTypes bh,
286 int verb) {
287
288 Interface &m_field = cOre;
290
291 auto get_side_map_hcurl = [&]() {
292 return std::vector<
293
294 std::pair<EntityType,
296
297 >>{
298
299 {MBTRI,
300 [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
302 dofs_side_map);
303 }}
304
305 };
306 };
307
308 auto get_side_map_hdiv = [&]() {
309 return std::vector<
310
311 std::pair<EntityType,
313
314 >>{
315
316 {MBTET,
317 [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
319 dofs_side_map);
320 }}
321
322 };
323 };
324
325 auto get_side_map = [&](auto space) {
326 switch (space) {
327 case HCURL:
328 return get_side_map_hcurl();
329 case HDIV:
330 return get_side_map_hdiv();
331 default:
332 MOFEM_LOG_CHANNEL("WORLD");
333 MOFEM_LOG("WORLD", Sev::warning)
334 << "Side dof map for broken space <" << space << "> not implemented";
335 }
336 return std::vector<
337
338 std::pair<EntityType,
340
341 >>{};
342 };
343
344 CHKERR m_field.add_broken_field(name, space, base, nb_of_coefficients,
345 get_side_map(space), tag_type, bh, verb);
346 domainFields.push_back(name);
347 if (space == NOFIELD)
348 SETERRQ(
349 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
350 "NOFIELD space for boundary field not implemented in Simple interface");
352}
353
355Simple::addBoundaryField(const std::string name, const FieldSpace space,
356 const FieldApproximationBase base,
357 const FieldCoefficientsNumber nb_of_coefficients,
358 const TagType tag_type, const enum MoFEMTypes bh,
359 int verb) {
360 Interface &m_field = cOre;
362 CHKERR m_field.add_field(name, space, base, nb_of_coefficients, tag_type, bh,
363 verb);
364 boundaryFields.push_back(name);
365 if (space == NOFIELD)
366 SETERRQ(
367 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
368 "NOFIELD space for boundary field not implemented in Simple interface");
370}
371
373Simple::addSkeletonField(const std::string name, const FieldSpace space,
374 const FieldApproximationBase base,
375 const FieldCoefficientsNumber nb_of_coefficients,
376 const TagType tag_type, const enum MoFEMTypes bh,
377 int verb) {
378
379 Interface &m_field = cOre;
381 CHKERR m_field.add_field(name, space, base, nb_of_coefficients, tag_type, bh,
382 verb);
383 skeletonFields.push_back(name);
384 if (space == NOFIELD)
385 SETERRQ(
386 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
387 "NOFIELD space for boundary field not implemented in Simple interface");
388
390}
391
393Simple::addDataField(const std::string name, const FieldSpace space,
394 const FieldApproximationBase base,
395 const FieldCoefficientsNumber nb_of_coefficients,
396 const TagType tag_type, const enum MoFEMTypes bh,
397 int verb) {
398
399 Interface &m_field = cOre;
401 CHKERR m_field.add_field(name, space, base, nb_of_coefficients, tag_type, bh,
402 verb);
403 if (space == NOFIELD)
404 noFieldDataFields.push_back(name);
405 else
406 dataFields.push_back(name);
408}
409
411Simple::addMeshsetField(const std::string name, const FieldSpace space,
412 const FieldApproximationBase base,
413 const FieldCoefficientsNumber nb_of_coefficients,
414 const TagType tag_type, const enum MoFEMTypes bh,
415 int verb) {
416
417 Interface &m_field = cOre;
419 CHKERR m_field.add_field(name, space, base, nb_of_coefficients, tag_type, bh,
420 verb);
421 if (space != NOFIELD) {
422 meshsetFields.push_back(name);
423 } else {
424 noFieldFields.push_back(name);
425 }
427}
428
431
432 auto remove_field_from_list = [&](auto &vec) {
433 auto it = std::find(vec.begin(), vec.end(), name);
434 if (it != vec.end())
435 vec.erase(it);
436 };
437
438 remove_field_from_list(domainFields);
439
441}
442
445
446 auto remove_field_from_list = [&](auto &vec) {
447 auto it = std::find(vec.begin(), vec.end(), name);
448 if (it != vec.end())
449 vec.erase(it);
450 };
451
452 remove_field_from_list(boundaryFields);
453
455}
456
459
460 auto remove_field_from_list = [&](auto &vec) {
461 auto it = std::find(vec.begin(), vec.end(), name);
462 if (it != vec.end())
463 vec.erase(it);
464 };
465
466 remove_field_from_list(skeletonFields);
467
469}
470
472 Interface &m_field = cOre;
474
475 auto clear_rows_and_cols = [&](auto &fe_name) {
477 auto fe_ptr = m_field.get_finite_elements();
478 auto &fe_by_name = const_cast<FiniteElement_multiIndex *>(fe_ptr)
479 ->get<FiniteElement_name_mi_tag>();
480 auto it_fe = fe_by_name.find(fe_name);
481 if (it_fe != fe_by_name.end()) {
482
483 if (!fe_by_name.modify(it_fe, FiniteElement_row_change_bit_reset()))
484 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
485 "modification unsuccessful");
486
487 if (!fe_by_name.modify(it_fe, FiniteElement_col_change_bit_reset()))
488 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
489 "modification unsuccessful");
490 }
492 };
493 CHKERR clear_rows_and_cols(domainFE);
494 CHKERR clear_rows_and_cols(boundaryFE);
495 CHKERR clear_rows_and_cols(skeletonFE);
496
497 // Define finite elements
499
500 auto add_fields = [&](auto &fe_name, auto &fields) {
502 for (auto &field : fields) {
503 CHKERR m_field.modify_finite_element_add_field_row(fe_name, field);
504 CHKERR m_field.modify_finite_element_add_field_col(fe_name, field);
505 CHKERR m_field.modify_finite_element_add_field_data(fe_name, field);
506 }
508 };
509
510 auto add_data_fields = [&](auto &fe_name, auto &fields) {
512 for (auto &field : fields)
513 CHKERR m_field.modify_finite_element_add_field_data(fe_name, field);
515 };
516
517 CHKERR add_fields(domainFE, domainFields);
518 CHKERR add_data_fields(domainFE, dataFields);
519 CHKERR add_fields(domainFE, noFieldFields);
520 CHKERR add_data_fields(domainFE, noFieldDataFields);
521
522 if (addBoundaryFE || !boundaryFields.empty()) {
524 CHKERR add_fields(boundaryFE, domainFields);
525 if (!boundaryFields.empty())
526 CHKERR add_fields(boundaryFE, boundaryFields);
527 CHKERR add_data_fields(boundaryFE, dataFields);
528 CHKERR add_data_fields(boundaryFE, noFieldDataFields);
529 CHKERR add_fields(boundaryFE, noFieldFields);
530 }
531 if (addSkeletonFE || !skeletonFields.empty()) {
533 CHKERR add_fields(skeletonFE, domainFields);
534 if (!skeletonFields.empty())
535 CHKERR add_fields(skeletonFE, skeletonFields);
536 CHKERR add_data_fields(skeletonFE, dataFields);
537 CHKERR add_data_fields(skeletonFE, noFieldDataFields);
538 CHKERR add_fields(skeletonFE, noFieldFields);
539 }
540
541 if (noFieldFields.size() || meshsetFields.size()) {
543 CHKERR add_fields(meshsetFE, meshsetFields);
544 CHKERR add_fields(meshsetFE, noFieldFields);
545 CHKERR add_data_fields(meshsetFE, noFieldDataFields);
546 }
547
549}
550
551MoFEMErrorCode Simple::defineProblem(const PetscBool is_partitioned) {
552 Interface &m_field = cOre;
554 // Create dm instance
555 dM = createDM(m_field.get_comm(), "DMMOFEM");
556 // set dm data structure which created mofem data structures
559 CHKERR DMSetFromOptions(dM);
561 if (addBoundaryFE || !boundaryFields.empty()) {
563 }
564 if (addSkeletonFE || !skeletonFields.empty()) {
566 }
567 if (noFieldFields.size() || meshsetFields.size()) {
569 }
571 CHKERR DMMoFEMSetIsPartitioned(dM, is_partitioned);
573}
574
576 const int order, const Range *ents) {
578 fieldsOrder.emplace_back(field_name, order,
579 ents == NULL ? Range() : Range(*ents),
580 ents == NULL ? false : true);
582}
583
585 Interface &m_field = cOre;
587 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildFields, 0, 0, 0, 0);
588
589 auto comm_interface_ptr = m_field.getInterface<CommInterface>();
590 auto bit_ref_ptr = m_field.getInterface<BitRefManager>();
591
592 // Add entities to the fields
593 auto add_ents_to_field = [&](auto meshset, auto dim, auto &fields) {
595 for (auto &field : fields) {
596 CHKERR m_field.add_ents_to_field_by_dim(meshset, dim, field);
597 CHKERR comm_interface_ptr->synchroniseFieldEntities(field, 0);
598 }
600 };
601
602 CHKERR add_ents_to_field(meshSet, dIm, domainFields);
603 CHKERR add_ents_to_field(meshSet, dIm, dataFields);
604 CHKERR add_ents_to_field(boundaryMeshset, dIm - 1, boundaryFields);
605 CHKERR add_ents_to_field(skeletonMeshset, dIm - 1, skeletonFields);
606
607 std::set<std::string> nofield_fields;
608 for (auto &f : noFieldFields)
609 nofield_fields.insert(f);
610 for (auto &f : noFieldDataFields)
611 nofield_fields.insert(f);
612
613 // Set order
614
615 auto get_field_ptr = [&](auto &f) { return m_field.get_field_structure(f); };
616 for (auto &t : fieldsOrder) {
617 const auto f = std::get<0>(t);
618 const auto order = std::get<1>(t);
619
620 MOFEM_TAG_AND_LOG("WORLD", Sev::inform, "Simple")
621 << "Set order to field " << f << " order " << order;
622 MOFEM_LOG_CHANNEL("SYNC");
623 if (std::get<3>(t)) {
624 MOFEM_TAG_AND_LOG("SYNC", Sev::verbose, "Simple")
625 << "To ents: " << std::endl
626 << std::get<2>(t) << std::endl;
627 }
628 MOFEM_LOG_SEVERITY_SYNC(m_field.get_comm(), Sev::verbose);
629
630 if (std::get<3>(t)) {
631
632 CHKERR m_field.set_field_order(std::get<2>(t), f, order);
633
634 } else {
635 auto f_ptr = get_field_ptr(f);
636
637 if (f_ptr->getSpace() == H1) {
638 if (f_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
639 CHKERR m_field.set_field_order(meshSet, MBVERTEX, f, order);
640 } else {
641 CHKERR m_field.set_field_order(meshSet, MBVERTEX, f, 1);
642 }
643 }
644
645 for (auto d = 1; d <= dIm; ++d) {
646 for (EntityType t = CN::TypeDimensionMap[d].first;
647 t <= CN::TypeDimensionMap[d].second; ++t) {
648 CHKERR m_field.set_field_order(meshSet, t, f, order);
649 }
650 }
651 }
652 }
653 MOFEM_LOG_CHANNEL("WORLD");
654 // Build fields
655 CHKERR m_field.build_fields();
656 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildFields, 0, 0, 0, 0);
658}
659
661 Interface &m_field = cOre;
663 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildFiniteElements, 0, 0, 0, 0);
664 // Add finite elements
666 true);
668 if (addBoundaryFE || boundaryFields.size()) {
670 boundaryFE, true);
672 }
673 if (addSkeletonFE || skeletonFields.size()) {
675 skeletonFE, true);
677 }
678
679 if (noFieldFields.size() || meshsetFields.size()) {
680
681 auto create_meshset = [&]() {
682 EntityHandle meshset;
683 auto create = [&]() {
685 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
686 {
687 ParallelComm *pcomm =
688 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
689 Tag part_tag = pcomm->part_tag();
690 int rank = pcomm->rank();
691 CHKERR m_field.get_moab().tag_set_data(part_tag, &meshset, 1,
692 &rank);
693 CHKERR m_field.getInterface<BitRefManager>()->setBitLevelToMeshset(
694 meshset, BitRefLevel().set());
695 }
697 };
698 CHK_THROW_MESSAGE(create(), "Failed to create meshset");
699 return meshset;
700 };
701
702 auto fe_meshset = create_meshset();
704 meshsetFE, false);
705 // Add entities to the meshset, that will make DOFs on those entities
706 // accessible on meshset finite element
707 for (auto ents : meshsetFiniteElementEntities) {
708 CHKERR m_field.get_moab().add_entities(fe_meshset, ents);
709 }
711
712 }
713
714 for (std::vector<std::string>::iterator fit = otherFEs.begin();
715 fit != otherFEs.end(); ++fit) {
716 CHKERR m_field.build_finite_elements(*fit);
717 }
718 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildFiniteElements, 0, 0, 0, 0);
720}
721
723 Interface &m_field = cOre;
725 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
727 // Set problem by the DOFs on the fields rather that by adding DOFs on the
728 // elements
729 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
731 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
732 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
734}
735
736MoFEMErrorCode Simple::setUp(const PetscBool is_partitioned) {
738
739 PetscLogEventBegin(MOFEM_EVENT_SimpleSetUP, 0, 0, 0, 0);
740
742
743 if (addSkeletonFE || !skeletonFields.empty()) {
745 if (addBoundaryFE || !boundaryFields.empty()) {
747 }
748 }
749
752
753 CHKERR defineProblem(is_partitioned);
757
758 PetscLogEventEnd(MOFEM_EVENT_SimpleSetUP, 0, 0, 0, 0);
760}
761
763 Interface &m_field = cOre;
765 PetscLogEventBegin(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
766
767 if (!only_dm) {
770 if (addSkeletonFE || !skeletonFields.empty()) {
772 if (addBoundaryFE || !boundaryFields.empty()) {
774 }
775 }
779 }
780
782
783 const Problem *problem_ptr;
784 CHKERR DMMoFEMGetProblemPtr(dM, &problem_ptr);
785 const auto problem_name = problem_ptr->getName();
786 CHKERR m_field.modify_problem_ref_level_set_bit(problem_name, bitLevel);
789
790 // Set problem by the DOFs on the fields rather that by adding DOFs on the
791 // elements
792 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
794 m_field.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
795
796 PetscLogEventEnd(MOFEM_EVENT_SimpleBuildProblem, 0, 0, 0, 0);
798}
799
802 CHKERR PetscObjectReference(getPetscObject(dM.get()));
803 *dm = dM.get();
805}
806
807/**
808 * @brief Delete dm
809 *
810 * @return MoFEMErrorCode
811 */
817
818/**
819 * @brief Delete finite elements
820 *
821 * @return MoFEMErrorCode
822 */
824 Interface &m_field = cOre;
826 for (auto fe : {domainFE, boundaryFE, skeletonFE}) {
827 if (m_field.check_finite_element(fe)) {
828 CHKERR m_field.delete_finite_element(fe);
829 }
830 }
832}
833
835Simple::addFieldToEmptyFieldBlocks(const std::string row_field,
836 const std::string col_field) const {
837 Interface &m_field = cOre;
840 getProblemName(), row_field, col_field);
842}
843
845 Interface &m_field = cOre;
847 ParallelComm *pcomm =
848 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
849
850 auto get_skin = [&](auto meshset) {
851 // filter not owned entities, those are not on boundary
852
853 Range domain_ents;
854 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, dIm,
855 domain_ents, true);
856 CHKERR pcomm->filter_pstatus(domain_ents,
857 PSTATUS_SHARED | PSTATUS_MULTISHARED,
858 PSTATUS_NOT, -1, nullptr);
859
860 Skinner skin(&m_field.get_moab());
861 Range domain_skin;
862 CHKERR skin.find_skin(0, domain_ents, false, domain_skin);
863 CHKERR pcomm->filter_pstatus(domain_skin,
864 PSTATUS_SHARED | PSTATUS_MULTISHARED,
865 PSTATUS_NOT, -1, nullptr);
866 return domain_skin;
867 };
868
869 auto create_boundary_meshset = [&](auto &&domain_skin) {
871 // create boundary meshset
872 if (boundaryMeshset != 0) {
874 }
875 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, boundaryMeshset);
876 CHKERR m_field.get_moab().add_entities(boundaryMeshset, domain_skin);
877 for (int dd = 0; dd != dIm - 1; dd++) {
878 Range adj;
879 CHKERR m_field.get_moab().get_adjacencies(domain_skin, dd, false, adj,
880 moab::Interface::UNION);
881 CHKERR m_field.get_moab().add_entities(boundaryMeshset, adj);
882 }
884 };
885
886 CHKERR create_boundary_meshset(get_skin(meshSet));
887
889}
890
892 Interface &m_field = cOre;
894
895 ParallelComm *pcomm =
896 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
897
898 auto create_skeleton_meshset = [&](auto meshset) {
900 // create boundary meshset
901 if (skeletonMeshset != 0) {
903 }
904 Range boundary_ents, skeleton_ents;
905 CHKERR m_field.get_moab().get_entities_by_dimension(boundaryMeshset,
906 dIm - 1, boundary_ents);
907 Range domain_ents;
908 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, dIm,
909 domain_ents, true);
910 CHKERR m_field.get_moab().get_adjacencies(
911 domain_ents, dIm - 1, true, skeleton_ents, moab::Interface::UNION);
912 skeleton_ents = subtract(skeleton_ents, boundary_ents);
913 CHKERR pcomm->filter_pstatus(skeleton_ents, PSTATUS_NOT_OWNED, PSTATUS_NOT,
914 -1, nullptr);
915 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, skeletonMeshset);
916 CHKERR m_field.get_moab().add_entities(skeletonMeshset, skeleton_ents);
918 };
919
920 CHKERR create_skeleton_meshset(meshSet);
921
923}
924
926 Interface &m_field = cOre;
928 MOFEM_LOG("WORLD", Sev::verbose) << "Exchange ghost cells";
929
930 ParallelComm *pcomm =
931 ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
932 if (pcomm == NULL)
933 pcomm = new ParallelComm(&m_field.get_moab(), m_field.get_comm());
934
935 CHKERR pcomm->exchange_ghost_cells(getDim(), getDim() - 1, 1,
936 3 /**get all adjacent ghosted entities */,
937 true, true, meshSet ? &meshSet : nullptr);
938
939 Range shared;
940 CHKERR m_field.get_moab().get_entities_by_dimension(0, dIm, shared);
941 for (auto d = dIm - 1; d >= 0; --d) {
942 CHKERR m_field.get_moab().get_adjacencies(shared, d, false, shared,
943 moab::Interface::UNION);
944 }
945 CHKERR pcomm->filter_pstatus(shared, PSTATUS_SHARED | PSTATUS_MULTISHARED,
946 PSTATUS_OR, -1, &shared);
947 Tag part_tag = pcomm->part_tag();
948 CHKERR pcomm->exchange_tags(part_tag, shared);
949 CHKERR m_field.getInterface<CommInterface>()->resolveParentEntities(shared,
950 VERBOSE);
951
953}
954
955} // namespace MoFEM
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
@ VERBOSE
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
@ MF_ZERO
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
Definition definitions.h:82
@ 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
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int order
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
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 DMMoFEM.cpp:114
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:422
PetscErrorCode DMSetUp_MoFEM(DM dm)
Definition DMMoFEM.cpp:1297
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 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_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
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 add_ents_to_finite_element_by_MESHSET(const EntityHandle meshset, const std::string &name, const bool recursive=false)=0
add MESHSET element to finite element database given by name
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 modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
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]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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 Common.hpp:10
PetscObject getPetscObject(T obj)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode add_broken_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const std::vector< std::pair< EntityType, std::function< MoFEMErrorCode(BaseFunction::DofsSideMap &)> > > list_dof_side_map, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
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 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 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 rebuild_database(int verb=DEFAULT_VERBOSITY)=0
Clear database and initialize it once again.
Core (interface) class.
Definition Core.hpp:82
Deprecated interface functions.
keeps basic data about problem
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
boost::shared_ptr< ParentFiniteElementAdjacencyFunction< 2 > > parentAdjFunctionDim2
Definition Simple.hpp:624
EntityHandle boundaryMeshset
meshset with boundary
Definition Simple.hpp:572
MoFEMErrorCode buildProblem()
Build problem.
Definition Simple.cpp:722
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:261
std::vector< std::string > dataFields
Data fields.
Definition Simple.hpp:587
MoFEMErrorCode setParentAdjacency()
Definition Simple.cpp:67
std::list< std::tuple< std::string, int, Range, bool > > fieldsOrder
Definition Simple.hpp:593
std::vector< std::string > domainFields
domain fields
Definition Simple.hpp:584
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition Simple.cpp:471
MoFEM::Core & cOre
Definition Simple.hpp:559
BitRefLevel bitAdjEnt
bit ref level for parent
Definition Simple.hpp:581
char meshFileName[255]
mesh file name
Definition Simple.hpp:606
EntityHandle meshSet
domain meshset
Definition Simple.hpp:571
int getDim() const
Get the problem dimension.
Definition Simple.hpp:333
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
std::vector< std::string > otherFEs
Other finite elements.
Definition Simple.hpp:602
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition Simple.cpp:660
MoFEMErrorCode addMeshsetField(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 meshset field.
Definition Simple.cpp:411
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
Definition Simple.hpp:629
MoFEMErrorCode exchangeGhostCells()
Definition Simple.cpp:925
PetscLogEvent MOFEM_EVENT_SimpleBuildFields
Definition Simple.hpp:566
PetscLogEvent MOFEM_EVENT_SimpleLoadMesh
Definition Simple.hpp:565
MoFEMErrorCode reSetUp(bool only_dm=false)
Rebuild internal MoFEM data structures.
Definition Simple.cpp:762
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:355
std::vector< std::string > noFieldFields
NOFIELD field name.
Definition Simple.hpp:589
std::string meshsetFE
meshset finite element
Definition Simple.hpp:600
bool addSkeletonFE
Add skeleton FE.
Definition Simple.hpp:575
BitRefLevel bitAdjParent
bit ref level for parent
Definition Simple.hpp:579
bool addBoundaryFE
Add boundary FE.
Definition Simple.hpp:576
MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string row_field, const std::string col_field) const
add empty block to problem
Definition Simple.cpp:835
std::string skeletonFE
skeleton finite element
Definition Simple.hpp:599
std::vector< std::string > noFieldDataFields
NOFIELD field name.
Definition Simple.hpp:590
std::vector< std::string > meshsetFields
meshset fields
Definition Simple.hpp:588
std::string boundaryFE
boundary finite element
Definition Simple.hpp:598
std::string nameOfProblem
problem name
Definition Simple.hpp:596
boost::shared_ptr< ParentFiniteElementAdjacencyFunction< 3 > > parentAdjFunctionDim3
Definition Simple.hpp:622
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
std::string domainFE
domain finite element
Definition Simple.hpp:597
std::vector< std::string > boundaryFields
boundary fields
Definition Simple.hpp:585
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 1 > > parentAdjSkeletonFunctionDim1
Definition Simple.hpp:631
MoFEMErrorCode removeSkeletonField(const std::string name)
Remove field form skeleton.
Definition Simple.cpp:457
BitRefLevel bitAdjEntMask
bit ref level for parent parent
Definition Simple.hpp:582
SmartPetscObj< DM > dM
Discrete manager (interface to PETSc/MoFEM functions)
Definition Simple.hpp:610
MoFEMErrorCode removeDomainField(const std::string name)
Remove field form domain.
Definition Simple.cpp:429
MoFEMErrorCode buildFields()
Build fields.
Definition Simple.cpp:584
SmartPetscObj< DM > getDM()
Return smart DM object.
Definition Simple.hpp:323
std::vector< Range > meshsetFiniteElementEntities
Meshset element entities.
Definition Simple.hpp:604
BitRefLevel bitLevel
BitRefLevel of the problem.
Definition Simple.hpp:561
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition Simple.cpp:551
std::vector< std::string > skeletonFields
fields on the skeleton
Definition Simple.hpp:586
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:393
boost::function< MoFEMErrorCode(Interface &, const char *, const char *)> LoadFileFunc
Definition Simple.hpp:43
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition Simple.cpp:8
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:736
MoFEMErrorCode removeBoundaryField(const std::string name)
Remove field form boundary.
Definition Simple.cpp:443
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:410
int dIm
dimension of problem
Definition Simple.hpp:607
MoFEMErrorCode setSkeletonAdjacency(int dim=-1, std::string fe_name="")
Set the skeleton adjacency object.
Definition Simple.cpp:143
MoFEMErrorCode deleteDM()
Delete dm.
Definition Simple.cpp:812
PetscLogEvent MOFEM_EVENT_SimpleSetUP
Definition Simple.hpp:564
PetscLogEvent MOFEM_EVENT_SimpleKSPSolve
Definition Simple.hpp:569
MoFEMErrorCode createBoundaryMeshset()
Definition Simple.cpp:844
PetscLogEvent MOFEM_EVENT_SimpleBuildFiniteElements
Definition Simple.hpp:567
BitRefLevel bitLevelMask
BitRefLevel of the problem.
Definition Simple.hpp:562
BitRefLevel bitAdjParentMask
bit ref level for parent parent
Definition Simple.hpp:580
PetscLogEvent MOFEM_EVENT_SimpleBuildProblem
Definition Simple.hpp:568
boost::shared_ptr< ParentFiniteElementAdjacencyFunction< 1 > > parentAdjFunctionDim1
Definition Simple.hpp:626
MoFEMErrorCode deleteFiniteElements()
Delete finite elements.
Definition Simple.cpp:823
Simple(const MoFEM::Core &core)
Definition Simple.cpp:162
EntityHandle skeletonMeshset
skeleton meshset with boundary
Definition Simple.hpp:573
MoFEMErrorCode createSkeletonMeshset()
Definition Simple.cpp:891
bool addParentAdjacencies
If set to true parent adjacencies are build.
Definition Simple.hpp:577
MoFEMErrorCode addDomainBrokenField(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 broken field on domain.
Definition Simple.cpp:282
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:373
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.