v0.13.1
ProblemsCore.cpp
Go to the documentation of this file.
1/** \file ProblemsCore.cpp
2 * \brief Managing complexities for problem
3 */
4
5
6#include <MoFEM.hpp>
7
8#define ProblemCoreFunctionBegin \
9 MoFEMFunctionBegin; \
10 MOFEM_LOG_CHANNEL("WORLD"); \
11 MOFEM_LOG_CHANNEL("SYNC"); \
12 MOFEM_LOG_FUNCTION(); \
13 MOFEM_LOG_TAG("SYNC", "ProblemCore"); \
14 MOFEM_LOG_TAG("WORLD", "ProblemCore")
15
16namespace MoFEM {
17
18bool Core::check_problem(const string name) {
19 Problem_multiIndex::index<Problem_mi_tag>::type::iterator pit;
20 pit = pRoblems.get<Problem_mi_tag>().find(name);
21 if (pit == pRoblems.get<Problem_mi_tag>().end()) {
22 return false;
23 }
24 return true;
25}
26
27MoFEMErrorCode Core::addProblem(const BitProblemId id, const std::string &name,
28 int verb) {
30
31 if (verb == -1)
32 verb = verbose;
33 EntityHandle meshset;
34 CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
35 CHKERR get_moab().tag_set_data(th_ProblemId, &meshset, 1, &id);
36
37 // Add problem meshset to partion meshset. In case of no elements
38 // on processor part, when mesh file is read, finite element meshset is
39 // prevented from deletion by moab reader.
40 auto add_meshset_to_partition = [&](auto meshset) {
42 const void *tag_vals[] = {&rAnk};
43 ParallelComm *pcomm = ParallelComm::get_pcomm(
44 &get_moab(), get_basic_entity_data_ptr()->pcommID);
45 Tag part_tag = pcomm->part_tag();
46 Range tagged_sets;
47 CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
48 tag_vals, 1, tagged_sets,
49 moab::Interface::UNION);
50 for (auto s : tagged_sets)
51 CHKERR get_moab().add_entities(s, &meshset, 1);
53 };
54 CHKERR add_meshset_to_partition(meshset);
55
56 void const *tag_data[] = {name.c_str()};
57 int tag_sizes[1];
58 tag_sizes[0] = name.size();
59 CHKERR get_moab().tag_set_by_ptr(th_ProblemName, &meshset, 1, tag_data,
60 tag_sizes);
61 // create entry
62 auto p = pRoblems.insert(Problem(moab, meshset));
63 if (!p.second) {
66 MOFEM_LOG("SELF", Sev::error) << "Following problem can not be added:";
67 MOFEM_LOG("SELF", Sev::error)
68 << "Problem " << name << " id(" << id.to_ulong() << ") " << id
69 << " added meshset " << meshset;
70 MOFEM_LOG("SELF", Sev::error) << "List of problems already in databse:";
71 for (auto &p : pRoblems) {
72 MOFEM_LOG("SELF", Sev::error)
73 << p.getName() << " id(" << p.getId().to_ulong() << ") " << id
74 << " on meshset " << p.meshset;
75 }
76 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problem not added");
77 }
78
79 MOFEM_LOG("WORLD", Sev::inform) << "Add problem " << name;
80
82}
83
84MoFEMErrorCode Core::add_problem(const std::string &name, enum MoFEMTypes bh,
85 int verb) {
87 if (verb == -1)
88 verb = verbose;
89 auto miit = pRoblems.get<Problem_mi_tag>().find(name);
90 if (miit == pRoblems.get<Problem_mi_tag>().end()) {
91
92 int p_shift = 0;
93 for (; pRoblems.get<BitProblemId_mi_tag>().find(BitProblemId().set(
94 p_shift)) != pRoblems.get<BitProblemId_mi_tag>().end();
95 ++p_shift) {
96 }
97
98 auto id = BitProblemId().set(p_shift);
99 CHKERR addProblem(id, name, verb);
100
101 } else if (bh == MF_EXCL) {
102 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem is in database %s",
103 name.c_str());
104 }
106}
107
108MoFEMErrorCode Core::delete_problem(const std::string name) {
110 auto p_miit = pRoblems.get<Problem_mi_tag>().find(name);
111 if (p_miit == pRoblems.get<Problem_mi_tag>().end()) {
112 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "no such problem like < %s >",
113 name.c_str());
114 }
115 const EntityHandle meshset = p_miit->meshset;
116 pRoblems.get<Problem_mi_tag>().erase(p_miit);
117 CHKERR get_moab().delete_entities(&meshset, 1);
119}
120
121BitProblemId Core::getBitProblemId(const std::string &name) const {
122 auto p_miit = pRoblems.get<Problem_mi_tag>().find(name);
123 if (p_miit == pRoblems.get<Problem_mi_tag>().end()) {
124 THROW_MESSAGE("no such problem like " + name + " >");
125 }
126 return p_miit->getId();
127}
128
132 const ProblemById &set_id = pRoblems.get<BitProblemId_mi_tag>();
133 ProblemById::iterator miit = set_id.begin();
134 for (; miit != set_id.end(); miit++) {
135 std::ostringstream ss;
136 ss << *miit << std::endl;
137 PetscPrintf(mofemComm, ss.str().c_str());
138 }
140}
141
143Core::modify_problem_add_finite_element(const std::string &name_problem,
144 const std::string &fe_name) {
147 ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
148 ProblemsByName::iterator miit = set.find(name_problem);
149 if (miit == set.end()) {
150 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is not there",
151 name_problem.c_str());
152 }
153 BitFEId f_id = getBitFEId(fe_name);
154 bool success = set.modify(miit, ProblemFiniteElementChangeBitAdd(f_id));
155 if (!success)
156 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
157 "modification unsuccessful");
159}
160
162Core::modify_problem_unset_finite_element(const std::string &name_problem,
163 const std::string &fe_name) {
166 ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
167 ProblemsByName::iterator miit = set.find(name_problem);
168 if (miit == set.end()) {
169 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is not there",
170 name_problem.c_str());
171 }
172 BitFEId f_id = getBitFEId(fe_name);
173 bool success = set.modify(miit, ProblemFiniteElementChangeBitUnSet(f_id));
174 if (!success)
175 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
176 "modification unsuccessful");
178}
179
181Core::modify_problem_ref_level_add_bit(const std::string &name_problem,
182 const BitRefLevel &bit) {
185 ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
186 ProblemsByName::iterator miit = set.find(name_problem);
187 std::ostringstream ss;
188 ss << name_problem;
189 if (miit == set.end())
190 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
191 ss.str().c_str());
192 bool success = set.modify(miit, ProblemChangeRefLevelBitAdd(bit));
193 if (!success)
194 SETERRQ(PETSC_COMM_SELF, 1, "modification unsuccessful");
196}
197
199Core::modify_problem_ref_level_set_bit(const std::string &name_problem,
200 const BitRefLevel &bit) {
203 ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
204 ProblemsByName::iterator miit = set.find(name_problem);
205 std::ostringstream ss;
206 ss << name_problem;
207 if (miit == set.end())
208 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
209 ss.str().c_str());
210 bool success = set.modify(miit, ProblemChangeRefLevelBitSet(bit));
211 if (!success)
212 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
213 "modification unsuccessful");
215}
216
218Core::modify_problem_mask_ref_level_add_bit(const std::string &name_problem,
219 const BitRefLevel &bit) {
222 ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
223 ProblemsByName::iterator miit = set.find(name_problem);
224 if (miit == set.end()) {
225 std::ostringstream ss;
226 ss << name_problem;
227 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
228 ss.str().c_str());
229 }
230 bool success = set.modify(miit, ProblemChangeRefLevelBitDofMaskAdd(bit));
231 if (!success)
232 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
233 "modification unsuccessful");
235}
236
238Core::modify_problem_mask_ref_level_set_bit(const std::string &name_problem,
239 const BitRefLevel &bit) {
241 ProblemsByName;
242 ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
243 ProblemsByName::iterator miit = set.find(name_problem);
244 if (miit == set.end()) {
245 std::ostringstream ss;
246 ss << name_problem;
247 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
248 ss.str().c_str());
249 }
250 bool success = set.modify(miit, ProblemChangeRefLevelBitDofMaskSet(bit));
251 if (!success)
252 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
253 "modification unsuccessful");
255}
256
259 if (verb == -1)
260 verb = verbose;
261 Problem_multiIndex::iterator p_miit = pRoblems.begin();
262 for (; p_miit != pRoblems.end(); p_miit++) {
263 CHKERR getInterface<ProblemsManager>()->buildProblemOnDistributedMesh(
264 const_cast<Problem *>(&*p_miit), verb);
265 }
267}
268
269MoFEMErrorCode Core::clear_problem(const std::string problem_name, int verb) {
271 if (verb == -1)
272 verb = verbose;
274 ProblemsByName &prob_by_name = pRoblems.get<Problem_mi_tag>();
275 ProblemsByName::iterator p_miit = prob_by_name.find(problem_name);
276 if (p_miit == prob_by_name.end()) {
277 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
278 "problem < %s > not found, (top tip: check spelling)",
279 problem_name.c_str());
280 }
281 // zero rows
282 bool success = prob_by_name.modify(p_miit, ProblemZeroNbRowsChange());
283 if (!success)
284 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
285 "modification unsuccessful");
286 // zero cols
287 success = prob_by_name.modify(p_miit, ProblemZeroNbColsChange());
288 if (!success)
289 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
290 "modification unsuccessful");
291 // clear finite elements
292 success =
293 prob_by_name.modify(p_miit, ProblemClearNumeredFiniteElementsChange());
294 if (!success)
295 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
296 "modification unsuccessful");
297 // clear data structures
298 success = prob_by_name.modify(p_miit, ProblemClearSubProblemData());
299 if (!success)
300 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
301 "modification unsuccessful");
302 success = prob_by_name.modify(p_miit, ProblemClearComposedProblemData());
303 if (!success)
304 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
305 "modification unsuccessful");
306 if (p_miit->getRowDofsSequence())
307 p_miit->getRowDofsSequence()->clear();
308 if (p_miit->getColDofsSequence())
309 p_miit->getColDofsSequence()->clear();
310 if (p_miit->getSubData())
311 p_miit->getSubData().reset();
312 if (p_miit->getComposedProblemsData())
313 p_miit->getComposedProblemsData().reset();
315}
316
319 if (verb == -1)
320 verb = verbose;
321 if (!((*buildMoFEM) & BUILD_FIELD))
322 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
323 if (!((*buildMoFEM) & BUILD_FE))
324 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
325 if (!((*buildMoFEM) & BUILD_ADJ))
326 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
327 // iterate problems
328 Problem_multiIndex::iterator p_miit = pRoblems.begin();
329 for (; p_miit != pRoblems.end(); p_miit++) {
330 Problem *problem_ptr = const_cast<Problem *>(&*p_miit);
331 CHKERR getInterface<ProblemsManager>()->buildProblem(problem_ptr, false,
332 verb);
333 }
334 *buildMoFEM |= BUILD_PROBLEM;
336}
337
340 if (verb == -1)
341 verb = verbose;
342 // iterate problems
343 for (auto p_miit = pRoblems.begin(); p_miit != pRoblems.end(); p_miit++)
344 CHKERR clear_problem(p_miit->getName(), verb);
346}
347
348#define SET_BASIC_METHOD(METHOD, PROBLEM_PTR) \
349 { \
350 METHOD.rAnk = rAnk; \
351 METHOD.sIze = sIze; \
352 METHOD.problemPtr = PROBLEM_PTR; \
353 METHOD.fieldsPtr = &fIelds; \
354 METHOD.refinedEntitiesPtr = &refinedEntities; \
355 METHOD.entitiesPtr = &entsFields; \
356 METHOD.dofsPtr = &dofsField; \
357 METHOD.refinedFiniteElementsPtr = &refinedFiniteElements; \
358 METHOD.finiteElementsPtr = &finiteElements; \
359 METHOD.finiteElementsEntitiesPtr = &entsFiniteElements; \
360 METHOD.adjacenciesPtr = &entFEAdjacencies; \
361 }
362
364 BasicMethod &method,
365 int verb) {
367 if (verb == -1)
368 verb = verbose;
369 // finite element
370 SET_BASIC_METHOD(method, problem_ptr)
371 PetscLogEventBegin(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
372 CHKERR method.preProcess();
373 PetscLogEventEnd(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
375}
376
378Core::problem_basic_method_preProcess(const std::string &problem_name,
379 BasicMethod &method, int verb) {
381 if (verb == -1)
382 verb = verbose;
384 // find p_miit
385 ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
386 ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
387 if (p_miit == pRoblems_set.end())
388 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem is not in database %s",
389 problem_name.c_str());
390 CHKERR problem_basic_method_preProcess(&*p_miit, method, verb);
392}
393
396 BasicMethod &method, int verb) {
398 SET_BASIC_METHOD(method, problem_ptr)
399
400 PetscLogEventBegin(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
401 CHKERR method.postProcess();
402 PetscLogEventEnd(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
403
405}
406
408Core::problem_basic_method_postProcess(const std::string &problem_name,
409 BasicMethod &method, int verb) {
411 if (verb == -1)
412 verb = verbose;
414
415 // find p_miit
416 ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
417 ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
418 if (p_miit == pRoblems_set.end())
419 SETERRQ1(mofemComm, 1, "problem is not in database %s",
420 problem_name.c_str());
421
422 CHKERR problem_basic_method_postProcess(&*p_miit, method, verb);
423
425}
426
428 const Problem *problem_ptr, const std::string &fe_name, FEMethod &method,
429 int lower_rank, int upper_rank,
430 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
431 CacheTupleWeakPtr cache_ptr, int verb) {
433 if (verb == DEFAULT_VERBOSITY)
434 verb = verbose;
435
436 CacheTupleSharedPtr tmp_cache_ptr;
437 if (!cache_ptr.use_count()) {
438 tmp_cache_ptr = boost::make_shared<CacheTuple>();
439 CHKERR cache_problem_entities(problem_ptr->getName(), tmp_cache_ptr);
440 method.cacheWeakPtr = tmp_cache_ptr;
441 } else {
442 method.cacheWeakPtr = cache_ptr;
443 }
444
445 if (!fe_ptr)
446 fe_ptr = problem_ptr->numeredFiniteElementsPtr;
447
448 auto miit = fe_ptr->get<Composite_Name_And_Part_mi_tag>().lower_bound(
449 boost::make_tuple(fe_name, lower_rank));
450 auto hi_miit = fe_ptr->get<Composite_Name_And_Part_mi_tag>().upper_bound(
451 boost::make_tuple(fe_name, upper_rank));
452
453 if (miit == hi_miit && (bh & MF_EXIST)) {
454 if (!check_finite_element(fe_name)) {
455 SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "finite element < %s > not found",
456 fe_name.c_str());
457 }
458 }
459
460 method.feName = fe_name;
461 method.loopSize =
462 std::distance(miit, hi_miit); // Set numbers of elements in the loop
463 method.loHiFERank = std::make_pair(lower_rank, upper_rank);
464
465 SET_BASIC_METHOD(method, &*problem_ptr)
466
467 PetscLogEventBegin(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
468 CHKERR method.preProcess();
469 PetscLogEventEnd(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
470
471 PetscLogEventBegin(MOFEM_EVENT_operator, 0, 0, 0, 0);
472 for (int nn = 0; miit != hi_miit; miit++, nn++) {
473
474 method.nInTheLoop = nn; // Index of element in the loop
475 method.numeredEntFiniteElementPtr = *miit;
476
477 if (method.exeTestHook) {
478 if (method.exeTestHook(&method)) {
479 CHKERR method();
480 }
481 } else {
482 CHKERR method();
483 }
484
485 }
486 PetscLogEventEnd(MOFEM_EVENT_operator, 0, 0, 0, 0);
487
488 PetscLogEventBegin(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
489 CHKERR method.postProcess();
490 PetscLogEventEnd(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
491
493}
494
496 const std::string &problem_name, const std::string &fe_name,
497 FEMethod &method,
498 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
499 CacheTupleWeakPtr cache_ptr, int verb) {
501 if (verb == DEFAULT_VERBOSITY)
502 verb = verbose;
503
504 CHKERR loop_finite_elements(problem_name, fe_name, method, rAnk, rAnk, fe_ptr,
505 bh, cache_ptr, verb);
506
508}
509
511 const std::string &problem_name, const std::string &fe_name,
512 FEMethod &method, int lower_rank, int upper_rank,
513 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
514 CacheTupleWeakPtr cache_ptr, int verb) {
516 if (verb == DEFAULT_VERBOSITY)
517 verb = verbose;
518
519 auto &prb_by_name = pRoblems.get<Problem_mi_tag>();
520 auto p_miit = prb_by_name.find(problem_name);
521 if (p_miit == prb_by_name.end())
522 SETERRQ1(mofemComm, MOFEM_INVALID_DATA, "Problem <%s> is not in database",
523 problem_name.c_str());
524
525 CHKERR loop_finite_elements(&*p_miit, fe_name, method, lower_rank, upper_rank,
526 fe_ptr, bh, cache_ptr, verb);
527
529}
530
532 const std::string &field_name, RowColData rc,
533 DofMethod &method, int lower_rank,
534 int upper_rank, int verb) {
536 SET_BASIC_METHOD(method, &*problem_ptr);
538 NumeredDofsByUId;
539 NumeredDofsByUId *dofs;
540 switch (rc) {
541 case ROW:
542 dofs = &problem_ptr->numeredRowDofsPtr->get<Unique_mi_tag>();
543 break;
544 case COL:
545 dofs = &problem_ptr->numeredColDofsPtr->get<Unique_mi_tag>();
546 break;
547 default:
548 SETERRQ(mofemComm, MOFEM_DATA_INCONSISTENCY, "Not implemented");
549 }
550
551 auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
552 if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
553 method.fieldPtr = *field_it;
554 } else {
555 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
556 ("Field not found " + field_name).c_str());
557 }
558
559 auto miit = dofs->lower_bound(
560 FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
561 auto hi_miit = dofs->upper_bound(
562 FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
563
564 method.loopSize = std::distance(miit, hi_miit);
565 method.loHiFERank = std::make_pair(lower_rank, upper_rank);
566
567 CHKERR method.preProcess();
568
569 int nn = 0;
570 for (; miit != hi_miit; miit++, nn++) {
571 if ((*miit)->getPart() >= lower_rank && (*miit)->getPart() <= upper_rank) {
572 method.nInTheLoop = nn; // Index of element in the loop
573 method.dofPtr = miit->get()->getDofEntityPtr();
574 method.dofNumeredPtr = *miit;
575 CHKERR method();
576 }
577 }
578
579 CHKERR method.postProcess();
581}
582
584 const std::string &problem_name, const std::string &field_name,
585 RowColData rc, // ROW or COL
586 DofMethod &method, // Finite element instance processed on each DOF
587 int lower_rank, // Only DOFs on processor higher or equal to this are
588 // processed
589 int upper_rank, // Only DOFs lowest or higher to this are processed
590 int verb // verbosity level
591) {
593 if (verb == DEFAULT_VERBOSITY)
594 verb = verbose;
596 // find p_miit
597 ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
598 ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
599 if (p_miit == pRoblems_set.end())
600 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
601 problem_name.c_str());
602 CHKERR loop_dofs(&*p_miit, field_name, rc, method, lower_rank, upper_rank,
603 verb);
605}
606
607MoFEMErrorCode Core::loop_dofs(const std::string &problem_name,
608 const std::string &field_name, RowColData rc,
609 DofMethod &method, int verb) {
611 if (verb == DEFAULT_VERBOSITY)
612 verb = verbose;
613 CHKERR loop_dofs(problem_name, field_name, rc, method, 0, sIze, verb);
615}
616
618 int verb) {
620 if (verb == DEFAULT_VERBOSITY)
621 verb = verbose;
622 SET_BASIC_METHOD(method, nullptr);
623
624 auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
625 if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
626 method.fieldPtr = *field_it;
627 } else {
628 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
629 ("Field not found " + field_name).c_str());
630 }
631
632 auto miit = dofsField.get<Unique_mi_tag>().lower_bound(
633 FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
634 auto hi_miit = dofsField.get<Unique_mi_tag>().upper_bound(
635 FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
636
637 method.loopSize = std::distance(miit, hi_miit);
638 method.loHiFERank = std::make_pair(0, get_comm_size());
639
640 CHKERR method.preProcess();
641 for (int nn = 0; miit != hi_miit; miit++, nn++) {
642 method.nInTheLoop = nn;
643 method.dofPtr = *miit;
644 CHKERR method();
645 }
646 CHKERR method.postProcess();
648}
649
651 const std::string field_name, RowColData rc,
652 EntityMethod &method, int lower_rank,
653 int upper_rank, int verb) {
655 if (verb == DEFAULT_VERBOSITY)
656 verb = verbose;
657 decltype(problem_ptr->numeredRowDofsPtr) dofs;
658 switch (rc) {
659 case ROW:
660 dofs = problem_ptr->numeredRowDofsPtr;
661 break;
662 case COL:
663 dofs = problem_ptr->numeredColDofsPtr;
664 break;
665 default:
666 SETERRQ(mofemComm, MOFEM_DATA_INCONSISTENCY,
667 "It works only with rows or columns");
668 }
669
670 auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
671 if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
672 method.fieldPtr = *field_it;
673 } else {
674 SETERRQ(mofemComm, MOFEM_NOT_FOUND,
675 ("Field not found " + field_name).c_str());
676 }
677
678 auto miit = dofs->lower_bound(
679 FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
680 auto hi_miit = dofs->upper_bound(
681 FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
682
683 using FieldEntity_view_multiIndex = multi_index_container<
684
685 boost::shared_ptr<FieldEntity>,
686 indexed_by<
687
688 ordered_unique<
689
690 tag<Ent_mi_tag>,
691 const_mem_fun<FieldEntity::interface_type_RefEntity, EntityHandle,
693
694 >>;
695
696 FieldEntity_view_multiIndex ents_view;
697 auto hint = ents_view.begin();
698 for (; miit != hi_miit; ++miit)
699 if ((*miit)->getPart() >= lower_rank && (*miit)->getPart() <= upper_rank)
700 ents_view.emplace_hint(hint, (*miit)->getFieldEntityPtr());
701
702 SET_BASIC_METHOD(method, problem_ptr);
703
704 method.loopSize = ents_view.size();
705 method.loHiFERank = std::make_pair(lower_rank, upper_rank);
706
707 CHKERR method.preProcess();
708 method.nInTheLoop = 0;
709 for (auto &field_ent : ents_view) {
710 method.entPtr = field_ent;
711 CHKERR method();
712 ++method.nInTheLoop;
713 }
714 CHKERR method.postProcess();
716}
717
718MoFEMErrorCode Core::loop_entities(const std::string problem_name,
719 const std::string field_name, RowColData rc,
720 EntityMethod &method, int lower_rank,
721 int upper_rank, int verb) {
723 if (verb == DEFAULT_VERBOSITY)
724 verb = verbose;
725 // find p_miit
726 auto &prb = pRoblems.get<Problem_mi_tag>();
727 auto p_miit = prb.find(problem_name);
728 if (p_miit == prb.end())
729 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
730 problem_name.c_str());
731 CHKERR loop_entities(&*p_miit, field_name, rc, method, lower_rank, upper_rank,
732 verb);
734}
735
736MoFEMErrorCode Core::loop_entities(const std::string problem_name,
737 const std::string field_name, RowColData rc,
738 EntityMethod &method, int verb) {
739 return loop_entities(problem_name, field_name, rc, method, rAnk, rAnk, verb);
740}
741
743 EntityMethod &method,
744 Range const *const ents, int verb) {
746 if (verb == DEFAULT_VERBOSITY)
747 verb = verbose;
748 SET_BASIC_METHOD(method, nullptr);
749
750 auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
751 if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
752 method.fieldPtr = *field_it;
753 } else {
754 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found %s",
755 field_name.c_str());
756 }
757
758 auto lo_eit = entsFields.get<Unique_mi_tag>().lower_bound(
759 FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
760 auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
761 FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
762
763 typedef multi_index_container<
764 boost::shared_ptr<FieldEntity>,
765 indexed_by<ordered_unique<
766 tag<Ent_mi_tag>, const_mem_fun<FieldEntity::interface_RefEntity,
767 EntityHandle, &FieldEntity::getEnt>>>>
768 FieldEntity_view_multiIndex;
769
770 FieldEntity_view_multiIndex ents_view;
771 ents_view.insert(lo_eit, hi_eit);
772
773 method.loopSize = ents_view.size();
774 method.loHiFERank = std::make_pair(0, get_comm_size());
775
776 CHKERR method.preProcess();
777 method.nInTheLoop = 0;
778
779 if (ents)
780 for (auto p = ents->const_pair_begin(); p != ents->const_pair_end(); ++p)
781 for (auto feit = ents_view.lower_bound(p->first);
782 feit != ents_view.upper_bound(p->second); ++feit) {
783 method.entPtr = *feit;
784 CHKERR method();
785 ++method.nInTheLoop;
786 }
787 else
788 for (auto &field_ent : ents_view) {
789 method.entPtr = field_ent;
790 CHKERR method();
791 ++method.nInTheLoop;
792 }
793
794 CHKERR method.postProcess();
796}
797
799 CacheTupleWeakPtr cache_weak_ptr) {
801
802 if (auto cache_ptr = cache_weak_ptr.lock()) {
803 auto p_miit = pRoblems.get<Problem_mi_tag>().find(prb_name);
804 if (p_miit == pRoblems.get<Problem_mi_tag>().end())
805 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
806 prb_name.c_str());
807
808 const BitRefLevel &prb_bit = p_miit->getBitRefLevel();
809 const BitRefLevel &prb_mask = p_miit->getBitRefLevelMask();
810 const BitFEId &prb_fe_id = p_miit->getBitFEId();
811
812 auto &row_dofs = p_miit->numeredRowDofsPtr;
813 auto &col_dofs = p_miit->numeredColDofsPtr;
814
815 auto &cache_data = std::get<0>(*cache_ptr);
816 auto &cache_row = std::get<1>(*cache_ptr);
817 auto &cache_col = std::get<2>(*cache_ptr);
818
819 cache_row.resize(entsFields.size());
820 if (row_dofs != col_dofs)
821 cache_col.resize(entsFields.size());
822 cache_data.resize(entsFields.size());
823
824 size_t idx = 0;
825 for (auto it = entsFields.begin(); it != entsFields.end(); ++it, ++idx) {
826
827 const auto uid = (*it)->getLocalUniqueId();
828 auto r = entFEAdjacencies.get<Unique_mi_tag>().equal_range(uid);
829 for (auto lo = r.first; lo != r.second; ++lo) {
830
831 if ((lo->getBitFEId() & prb_fe_id).any()) {
832
833 const BitRefLevel &fe_bit = lo->entFePtr->getBitRefLevel();
834
835 // if entity is not problem refinement level
836 if (((fe_bit & prb_mask) != fe_bit) || ((fe_bit & prb_bit).none()))
837 continue;
838
839 auto cache_numered_dofs = [&](auto &numered_dofs, auto &cache_vec,
840 auto &ent_cache) {
841 auto dit = numered_dofs->lower_bound(uid);
842 decltype(dit) hi_dit;
843 if (dit != numered_dofs->end())
844 hi_dit = numered_dofs->upper_bound(
845 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
846 else
847 hi_dit = dit;
848
849 ent_cache = boost::shared_ptr<EntityCacheNumeredDofs>(
850 cache_ptr, &(cache_vec[idx]));
851 cache_vec[idx].loHi = {dit, hi_dit};
852 };
853
854 auto cache_dofs = [&](auto &dofs, auto &cache_vec, auto &ent_cache) {
855 auto dit = dofs.lower_bound(uid);
856 decltype(dit) hi_dit;
857 if (dit != dofs.end())
858 hi_dit = dofs.upper_bound(
859 uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
860 else
861 hi_dit = dit;
862
863 ent_cache = boost::shared_ptr<EntityCacheDofs>(cache_ptr,
864 &(cache_vec[idx]));
865 cache_vec[idx].loHi = {dit, hi_dit};
866 };
867
868 cache_numered_dofs(row_dofs, cache_row, (*it)->entityCacheRowDofs);
869 if (row_dofs != col_dofs) {
870 if (cache_col.size() != entsFields.size())
871 cache_col.resize(entsFields.size());
872 cache_numered_dofs(col_dofs, cache_col, (*it)->entityCacheColDofs);
873 } else {
874 (*it)->entityCacheColDofs = (*it)->entityCacheRowDofs;
875 }
876
877 cache_dofs(dofsField, cache_data, (*it)->entityCacheDataDofs);
878
879 break;
880 }
881 }
882 }
883 } else {
884 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cache not allocated");
885 }
886
888}
889
890} // namespace MoFEM
static Index< 'p', 3 > p
#define ProblemCoreFunctionBegin
Definition: ProblemsCore.cpp:8
#define SET_BASIC_METHOD(METHOD, PROBLEM_PTR)
@ DEFAULT_VERBOSITY
Definition: definitions.h:207
RowColData
RowColData.
Definition: definitions.h:123
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:97
@ MF_EXIST
Definition: definitions.h:100
@ MF_EXCL
Definition: definitions.h:99
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#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
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:289
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:318
auto bit
set bit
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:43
std::bitset< BITPROBLEMID_SIZE > BitProblemId
Problem Id.
Definition: Types.hpp:44
uint128_t UId
Unique Id.
Definition: Types.hpp:31
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
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
boost::shared_ptr< CacheTuple > CacheTupleSharedPtr
boost::weak_ptr< CacheTuple > CacheTupleWeakPtr
const double r
rate factor
constexpr auto field_name
Data structure to exchange data between mofem and User Loop Methods.
int loopSize
local number oe methods to process
std::pair< int, int > loHiFERank
Llo and hi processor rank of iterated entities.
int nInTheLoop
number currently of processed method
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::weak_ptr< CacheTuple > cacheWeakPtr
MoFEMErrorCode modify_problem_mask_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)
set dof mask ref level for problem
MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)
Set data for BasicMethod.
MoFEMErrorCode loop_entities(const Problem *problem_ptr, const std::string field_name, RowColData rc, EntityMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)
Loop over field entities in the problem.
MoFEMErrorCode modify_problem_unset_finite_element(const std::string &name_problem, const std::string &MoFEMFiniteElement_name)
unset finite element from problem, this remove entities assigned to finite element to a particular pr...
MoFEMErrorCode list_problem() const
list problems
DEPRECATED MoFEMErrorCode build_problem_on_distributed_mesh(int verb=DEFAULT_VERBOSITY)
MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
Add problem.
MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)
Set data for BasicMethod.
MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)
add ref level to problem
MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)
Make a loop over finite elements on partitions from upper to lower rank.
MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)
clear problem
MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &MoFEMFiniteElement_name)
add finite element to problem, this add entities assigned to finite element to a particular problem
MoFEMErrorCode addProblem(const BitProblemId id, const std::string &name, int verb=DEFAULT_VERBOSITY)
add problem
MoFEMErrorCode clear_problems(int verb=DEFAULT_VERBOSITY)
clear problems
MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)
Make a loop over dofs.
MoFEMErrorCode modify_problem_mask_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)
set dof mask ref level for problem
bool check_problem(const std::string name)
check if problem exist
BitProblemId getBitProblemId(const std::string &name) const
DEPRECATED MoFEMErrorCode build_problems(int verb=DEFAULT_VERBOSITY)
MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)
Cache variables.
MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)
set ref level for problem
MoFEMErrorCode delete_problem(const std::string name)
Delete problem.
Data structure to exchange data between mofem and User Loop Methods on entities.
boost::shared_ptr< DofEntity > dofPtr
boost::shared_ptr< NumeredDofEntity > dofNumeredPtr
boost::shared_ptr< Field > fieldPtr
Data structure to exchange data between mofem and User Loop Methods on entities.
boost::shared_ptr< Field > fieldPtr
boost::shared_ptr< FieldEntity > entPtr
structure for User Loop Methods on finite elements
std::string feName
Name of finite element.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
boost::function< bool(FEMethod *fe_method_ptr)> exeTestHook
Tet if element to skip element.
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
interface_RefEntity< RefEntity > interface_type_RefEntity
MultiIndex Tag for field name.
Clear composed problem data structure.
Clear sub-problem data structure.
remove finite element from problem
keeps basic data about problem
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
interface_RefEntity(const boost::shared_ptr< RefEntity > &sptr)
EntityHandle getEnt() const
Get the entity handle.