v0.12.1
ProblemsCore.cpp
Go to the documentation of this file.
1 /** \file ProblemsCore.cpp
2  * \brief Managing complexities for problem
3  */
4 
5 /* MoFEM is free software: you can redistribute it and/or modify it under
6  * the terms of the GNU Lesser General Public License as published by the
7  * Free Software Foundation, either version 3 of the License, or (at your
8  * option) any later version.
9  *
10  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13  * License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
17  */
18 
19 #include <MoFEM.hpp>
20 
21 #define ProblemCoreFunctionBegin \
22  MoFEMFunctionBegin; \
23  MOFEM_LOG_CHANNEL("WORLD"); \
24  MOFEM_LOG_CHANNEL("SYNC"); \
25  MOFEM_LOG_FUNCTION(); \
26  MOFEM_LOG_TAG("SYNC", "ProblemCore"); \
27  MOFEM_LOG_TAG("WORLD", "ProblemCore")
28 
29 namespace MoFEM {
30 
31 bool Core::check_problem(const string name) {
32  Problem_multiIndex::index<Problem_mi_tag>::type::iterator pit;
33  pit = pRoblems.get<Problem_mi_tag>().find(name);
34  if (pit == pRoblems.get<Problem_mi_tag>().end()) {
35  return false;
36  }
37  return true;
38 }
39 
40 MoFEMErrorCode Core::addProblem(const BitProblemId id, const std::string &name,
41  int verb) {
43 
44  if (verb == -1)
45  verb = verbose;
46  EntityHandle meshset;
47  CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
48  CHKERR get_moab().tag_set_data(th_ProblemId, &meshset, 1, &id);
49 
50  // Add problem meshset to partion meshset. In case of no elements
51  // on processor part, when mesh file is read, finite element meshset is
52  // prevented from deletion by moab reader.
53  auto add_meshset_to_partition = [&](auto meshset) {
55  const void *tag_vals[] = {&rAnk};
56  ParallelComm *pcomm = ParallelComm::get_pcomm(
57  &get_moab(), get_basic_entity_data_ptr()->pcommID);
58  Tag part_tag = pcomm->part_tag();
59  Range tagged_sets;
60  CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
61  tag_vals, 1, tagged_sets,
62  moab::Interface::UNION);
63  for (auto s : tagged_sets)
64  CHKERR get_moab().add_entities(s, &meshset, 1);
66  };
67  CHKERR add_meshset_to_partition(meshset);
68 
69  void const *tag_data[] = {name.c_str()};
70  int tag_sizes[1];
71  tag_sizes[0] = name.size();
72  CHKERR get_moab().tag_set_by_ptr(th_ProblemName, &meshset, 1, tag_data,
73  tag_sizes);
74  // create entry
75  auto p = pRoblems.insert(Problem(moab, meshset));
76  if (!p.second) {
79  MOFEM_LOG("SELF", Sev::error) << "Following problem can not be added:";
80  MOFEM_LOG("SELF", Sev::error)
81  << "Problem " << name << " id(" << id.to_ulong() << ") " << id
82  << " added meshset " << meshset;
83  MOFEM_LOG("SELF", Sev::error) << "List of problems already in databse:";
84  for (auto &p : pRoblems) {
85  MOFEM_LOG("SELF", Sev::error)
86  << p.getName() << " id(" << p.getId().to_ulong() << ") " << id
87  << " on meshset " << p.meshset;
88  }
89  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problem not added");
90  }
91 
92  MOFEM_LOG("WORLD", Sev::inform) << "Add problem " << name;
93 
95 }
96 
97 MoFEMErrorCode Core::add_problem(const std::string &name, enum MoFEMTypes bh,
98  int verb) {
100  if (verb == -1)
101  verb = verbose;
102  auto miit = pRoblems.get<Problem_mi_tag>().find(name);
103  if (miit == pRoblems.get<Problem_mi_tag>().end()) {
104 
105  int p_shift = 0;
106  for (; pRoblems.get<BitProblemId_mi_tag>().find(BitProblemId().set(
107  p_shift)) != pRoblems.get<BitProblemId_mi_tag>().end();
108  ++p_shift) {
109  }
110 
111  auto id = BitProblemId().set(p_shift);
112  CHKERR addProblem(id, name, verb);
113 
114  } else if (bh == MF_EXCL) {
115  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem is in database %s",
116  name.c_str());
117  }
119 }
120 
121 MoFEMErrorCode Core::delete_problem(const std::string name) {
123  auto p_miit = pRoblems.get<Problem_mi_tag>().find(name);
124  if (p_miit == pRoblems.get<Problem_mi_tag>().end()) {
125  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "no such problem like < %s >",
126  name.c_str());
127  }
128  const EntityHandle meshset = p_miit->meshset;
129  pRoblems.get<Problem_mi_tag>().erase(p_miit);
130  CHKERR get_moab().delete_entities(&meshset, 1);
132 }
133 
134 BitProblemId Core::getBitProblemId(const std::string &name) const {
135  auto p_miit = pRoblems.get<Problem_mi_tag>().find(name);
136  if (p_miit == pRoblems.get<Problem_mi_tag>().end()) {
137  THROW_MESSAGE("no such problem like " + name + " >");
138  }
139  return p_miit->getId();
140 }
141 
145  const ProblemById &set_id = pRoblems.get<BitProblemId_mi_tag>();
146  ProblemById::iterator miit = set_id.begin();
147  for (; miit != set_id.end(); miit++) {
148  std::ostringstream ss;
149  ss << *miit << std::endl;
150  PetscPrintf(mofemComm, ss.str().c_str());
151  }
153 }
154 
156 Core::modify_problem_add_finite_element(const std::string &name_problem,
157  const std::string &fe_name) {
160  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
161  ProblemsByName::iterator miit = set.find(name_problem);
162  if (miit == set.end()) {
163  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is not there",
164  name_problem.c_str());
165  }
166  BitFEId f_id = getBitFEId(fe_name);
167  bool success = set.modify(miit, ProblemFiniteElementChangeBitAdd(f_id));
168  if (!success)
169  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
170  "modification unsuccessful");
172 }
173 
175 Core::modify_problem_unset_finite_element(const std::string &name_problem,
176  const std::string &fe_name) {
179  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
180  ProblemsByName::iterator miit = set.find(name_problem);
181  if (miit == set.end()) {
182  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is not there",
183  name_problem.c_str());
184  }
185  BitFEId f_id = getBitFEId(fe_name);
186  bool success = set.modify(miit, ProblemFiniteElementChangeBitUnSet(f_id));
187  if (!success)
188  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
189  "modification unsuccessful");
191 }
192 
194 Core::modify_problem_ref_level_add_bit(const std::string &name_problem,
195  const BitRefLevel &bit) {
198  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
199  ProblemsByName::iterator miit = set.find(name_problem);
200  std::ostringstream ss;
201  ss << name_problem;
202  if (miit == set.end())
203  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
204  ss.str().c_str());
205  bool success = set.modify(miit, ProblemChangeRefLevelBitAdd(bit));
206  if (!success)
207  SETERRQ(PETSC_COMM_SELF, 1, "modification unsuccessful");
209 }
210 
212 Core::modify_problem_ref_level_set_bit(const std::string &name_problem,
213  const BitRefLevel &bit) {
216  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
217  ProblemsByName::iterator miit = set.find(name_problem);
218  std::ostringstream ss;
219  ss << name_problem;
220  if (miit == set.end())
221  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
222  ss.str().c_str());
223  bool success = set.modify(miit, ProblemChangeRefLevelBitSet(bit));
224  if (!success)
225  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
226  "modification unsuccessful");
228 }
229 
231 Core::modify_problem_mask_ref_level_add_bit(const std::string &name_problem,
232  const BitRefLevel &bit) {
235  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
236  ProblemsByName::iterator miit = set.find(name_problem);
237  if (miit == set.end()) {
238  std::ostringstream ss;
239  ss << name_problem;
240  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
241  ss.str().c_str());
242  }
243  bool success = set.modify(miit, ProblemChangeRefLevelBitDofMaskAdd(bit));
244  if (!success)
245  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
246  "modification unsuccessful");
248 }
249 
251 Core::modify_problem_mask_ref_level_set_bit(const std::string &name_problem,
252  const BitRefLevel &bit) {
254  ProblemsByName;
255  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
256  ProblemsByName::iterator miit = set.find(name_problem);
257  if (miit == set.end()) {
258  std::ostringstream ss;
259  ss << name_problem;
260  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
261  ss.str().c_str());
262  }
263  bool success = set.modify(miit, ProblemChangeRefLevelBitDofMaskSet(bit));
264  if (!success)
265  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
266  "modification unsuccessful");
268 }
269 
272  if (verb == -1)
273  verb = verbose;
274  Problem_multiIndex::iterator p_miit = pRoblems.begin();
275  for (; p_miit != pRoblems.end(); p_miit++) {
276  CHKERR getInterface<ProblemsManager>()->buildProblemOnDistributedMesh(
277  const_cast<Problem *>(&*p_miit), verb);
278  }
280 }
281 
282 MoFEMErrorCode Core::clear_problem(const std::string problem_name, int verb) {
284  if (verb == -1)
285  verb = verbose;
287  ProblemsByName &prob_by_name = pRoblems.get<Problem_mi_tag>();
288  ProblemsByName::iterator p_miit = prob_by_name.find(problem_name);
289  if (p_miit == prob_by_name.end()) {
290  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
291  "problem < %s > not found, (top tip: check spelling)",
292  problem_name.c_str());
293  }
294  // zero rows
295  bool success = prob_by_name.modify(p_miit, ProblemZeroNbRowsChange());
296  if (!success)
297  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
298  "modification unsuccessful");
299  // zero cols
300  success = prob_by_name.modify(p_miit, ProblemZeroNbColsChange());
301  if (!success)
302  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
303  "modification unsuccessful");
304  // clear finite elements
305  success =
306  prob_by_name.modify(p_miit, ProblemClearNumeredFiniteElementsChange());
307  if (!success)
308  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
309  "modification unsuccessful");
310  // clear data structures
311  success = prob_by_name.modify(p_miit, ProblemClearSubProblemData());
312  if (!success)
313  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
314  "modification unsuccessful");
315  success = prob_by_name.modify(p_miit, ProblemClearComposedProblemData());
316  if (!success)
317  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
318  "modification unsuccessful");
319  if (p_miit->getRowDofsSequence())
320  p_miit->getRowDofsSequence()->clear();
321  if (p_miit->getColDofsSequence())
322  p_miit->getColDofsSequence()->clear();
323  if (p_miit->getSubData())
324  p_miit->getSubData().reset();
325  if (p_miit->getComposedProblemsData())
326  p_miit->getComposedProblemsData().reset();
328 }
329 
332  if (verb == -1)
333  verb = verbose;
334  if (!((*buildMoFEM) & BUILD_FIELD))
335  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
336  if (!((*buildMoFEM) & BUILD_FE))
337  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
338  if (!((*buildMoFEM) & BUILD_ADJ))
339  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
340  // iterate problems
341  Problem_multiIndex::iterator p_miit = pRoblems.begin();
342  for (; p_miit != pRoblems.end(); p_miit++) {
343  Problem *problem_ptr = const_cast<Problem *>(&*p_miit);
344  CHKERR getInterface<ProblemsManager>()->buildProblem(problem_ptr, false,
345  verb);
346  }
347  *buildMoFEM |= BUILD_PROBLEM;
349 }
350 
353  if (verb == -1)
354  verb = verbose;
355  // iterate problems
356  for (auto p_miit = pRoblems.begin(); p_miit != pRoblems.end(); p_miit++)
357  CHKERR clear_problem(p_miit->getName(), verb);
359 }
360 
361 #define SET_BASIC_METHOD(METHOD, PROBLEM_PTR) \
362  { \
363  METHOD.rAnk = rAnk; \
364  METHOD.sIze = sIze; \
365  METHOD.problemPtr = PROBLEM_PTR; \
366  METHOD.fieldsPtr = &fIelds; \
367  METHOD.refinedEntitiesPtr = &refinedEntities; \
368  METHOD.entitiesPtr = &entsFields; \
369  METHOD.dofsPtr = &dofsField; \
370  METHOD.refinedFiniteElementsPtr = &refinedFiniteElements; \
371  METHOD.finiteElementsPtr = &finiteElements; \
372  METHOD.finiteElementsEntitiesPtr = &entsFiniteElements; \
373  METHOD.adjacenciesPtr = &entFEAdjacencies; \
374  }
375 
377  BasicMethod &method,
378  int verb) {
380  if (verb == -1)
381  verb = verbose;
382  // finite element
383  SET_BASIC_METHOD(method, problem_ptr)
384  PetscLogEventBegin(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
385  CHKERR method.preProcess();
386  PetscLogEventEnd(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
388 }
389 
391 Core::problem_basic_method_preProcess(const std::string &problem_name,
392  BasicMethod &method, int verb) {
394  if (verb == -1)
395  verb = verbose;
397  // find p_miit
398  ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
399  ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
400  if (p_miit == pRoblems_set.end())
401  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem is not in database %s",
402  problem_name.c_str());
403  CHKERR problem_basic_method_preProcess(&*p_miit, method, verb);
405 }
406 
409  BasicMethod &method, int verb) {
411  SET_BASIC_METHOD(method, problem_ptr)
412 
413  PetscLogEventBegin(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
414  CHKERR method.postProcess();
415  PetscLogEventEnd(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
416 
418 }
419 
421 Core::problem_basic_method_postProcess(const std::string &problem_name,
422  BasicMethod &method, int verb) {
424  if (verb == -1)
425  verb = verbose;
427 
428  // find p_miit
429  ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
430  ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
431  if (p_miit == pRoblems_set.end())
432  SETERRQ1(mofemComm, 1, "problem is not in database %s",
433  problem_name.c_str());
434 
435  CHKERR problem_basic_method_postProcess(&*p_miit, method, verb);
436 
438 }
439 
441  const Problem *problem_ptr, const std::string &fe_name, FEMethod &method,
442  int lower_rank, int upper_rank,
443  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
444  CacheTupleWeakPtr cache_ptr, int verb) {
446  if (verb == DEFAULT_VERBOSITY)
447  verb = verbose;
448 
449  CacheTupleSharedPtr tmp_cache_ptr;
450  if (!cache_ptr.use_count()) {
451  tmp_cache_ptr = boost::make_shared<CacheTuple>();
452  CHKERR cache_problem_entities(problem_ptr->getName(), tmp_cache_ptr);
453  }
454 
455  method.feName = fe_name;
456  SET_BASIC_METHOD(method, &*problem_ptr)
457  PetscLogEventBegin(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
458  CHKERR method.preProcess();
459  PetscLogEventEnd(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
460 
461  if (!fe_ptr)
462  fe_ptr = problem_ptr->numeredFiniteElementsPtr;
463 
464  auto miit = fe_ptr->get<Composite_Name_And_Part_mi_tag>().lower_bound(
465  boost::make_tuple(fe_name, lower_rank));
466  auto hi_miit = fe_ptr->get<Composite_Name_And_Part_mi_tag>().upper_bound(
467  boost::make_tuple(fe_name, upper_rank));
468 
469  if (miit == hi_miit && (bh & MF_EXIST)) {
470  if (!check_finite_element(fe_name)) {
471  SETERRQ1(mofemComm, MOFEM_NOT_FOUND, "finite element < %s > not found",
472  fe_name.c_str());
473  }
474  }
475 
476  method.loopSize = std::distance(miit, hi_miit);
477  for (int nn = 0; miit != hi_miit; miit++, nn++) {
478 
479  method.nInTheLoop = nn;
480  method.numeredEntFiniteElementPtr = *miit;
481 
482  PetscLogEventBegin(MOFEM_EVENT_operator, 0, 0, 0, 0);
483  CHKERR method();
484  PetscLogEventEnd(MOFEM_EVENT_operator, 0, 0, 0, 0);
485  }
486 
487  PetscLogEventBegin(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
488  CHKERR method.postProcess();
489  PetscLogEventEnd(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
490 
492 }
493 
495  const std::string &problem_name, const std::string &fe_name,
496  FEMethod &method,
497  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
498  CacheTupleWeakPtr cache_ptr, int verb) {
500  if (verb == DEFAULT_VERBOSITY)
501  verb = verbose;
502 
503  CHKERR loop_finite_elements(problem_name, fe_name, method, rAnk, rAnk, fe_ptr,
504  bh, cache_ptr, verb);
505 
507 }
508 
510  const std::string &problem_name, const std::string &fe_name,
511  FEMethod &method, int lower_rank, int upper_rank,
512  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
513  CacheTupleWeakPtr cache_ptr, int verb) {
515  if (verb == DEFAULT_VERBOSITY)
516  verb = verbose;
517 
518  auto &prb_by_name = pRoblems.get<Problem_mi_tag>();
519  auto p_miit = prb_by_name.find(problem_name);
520  if (p_miit == prb_by_name.end())
521  SETERRQ1(mofemComm, MOFEM_INVALID_DATA, "Problem <%s> is not in database",
522  problem_name.c_str());
523 
524  CHKERR loop_finite_elements(&*p_miit, fe_name, method, lower_rank, upper_rank,
525  fe_ptr, bh, cache_ptr, verb);
526 
528 }
529 
531  const std::string &field_name, RowColData rc,
532  DofMethod &method, int lower_rank,
533  int upper_rank, int verb) {
535  SET_BASIC_METHOD(method, &*problem_ptr);
537  NumeredDofsByUId;
538  NumeredDofsByUId *dofs;
539  switch (rc) {
540  case ROW:
541  dofs = &problem_ptr->numeredRowDofsPtr->get<Unique_mi_tag>();
542  break;
543  case COL:
544  dofs = &problem_ptr->numeredColDofsPtr->get<Unique_mi_tag>();
545  break;
546  default:
547  SETERRQ(mofemComm, MOFEM_DATA_INCONSISTENCY, "Not implemented");
548  }
549 
550  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
551  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
552  method.fieldPtr = *field_it;
553  } else {
554  SETERRQ(mofemComm, MOFEM_NOT_FOUND,
555  ("Field not found " + field_name).c_str());
556  }
557 
558  auto miit = dofs->lower_bound(
559  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
560  auto hi_miit = dofs->upper_bound(
561  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
562 
563  CHKERR method.preProcess();
564  for (; miit != hi_miit; miit++) {
565  if ((*miit)->getPart() >= lower_rank && (*miit)->getPart() <= upper_rank) {
566  method.dofPtr = miit->get()->getDofEntityPtr();
567  method.dofNumeredPtr = *miit;
568  CHKERR method();
569  }
570  }
571  CHKERR method.postProcess();
573 }
574 
576  const std::string &problem_name, const std::string &field_name,
577  RowColData rc, // ROW or COL
578  DofMethod &method, // Finite element instance processed on each DOF
579  int lower_rank, // Only DOFs on processor higher or equal to this are
580  // processed
581  int upper_rank, // Only DOFs lowest or higher to this are processed
582  int verb // verbosity level
583 ) {
585  if (verb == DEFAULT_VERBOSITY)
586  verb = verbose;
588  // find p_miit
589  ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
590  ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
591  if (p_miit == pRoblems_set.end())
592  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
593  problem_name.c_str());
594  CHKERR loop_dofs(&*p_miit, field_name, rc, method, lower_rank, upper_rank,
595  verb);
597 }
598 
599 MoFEMErrorCode Core::loop_dofs(const std::string &problem_name,
600  const std::string &field_name, RowColData rc,
601  DofMethod &method, int verb) {
603  if (verb == DEFAULT_VERBOSITY)
604  verb = verbose;
605  CHKERR loop_dofs(problem_name, field_name, rc, method, 0, sIze, verb);
607 }
608 
609 MoFEMErrorCode Core::loop_dofs(const std::string &field_name, DofMethod &method,
610  int verb) {
612  if (verb == DEFAULT_VERBOSITY)
613  verb = verbose;
614  SET_BASIC_METHOD(method, nullptr);
615 
616  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
617  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
618  method.fieldPtr = *field_it;
619  } else {
620  SETERRQ(mofemComm, MOFEM_NOT_FOUND,
621  ("Field not found " + field_name).c_str());
622  }
623 
624  auto miit = dofsField.get<Unique_mi_tag>().lower_bound(
625  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
626  auto hi_miit = dofsField.get<Unique_mi_tag>().upper_bound(
627  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
628 
629  method.loopSize = std::distance(miit, hi_miit);
630  CHKERR method.preProcess();
631  for (int nn = 0; miit != hi_miit; miit++, nn++) {
632  method.nInTheLoop = nn;
633  method.dofPtr = *miit;
634  CHKERR method();
635  }
636  CHKERR method.postProcess();
638 }
639 
641  const std::string field_name, RowColData rc,
642  EntityMethod &method, int lower_rank,
643  int upper_rank, int verb) {
645  if (verb == DEFAULT_VERBOSITY)
646  verb = verbose;
647  decltype(problem_ptr->numeredRowDofsPtr) dofs;
648  switch (rc) {
649  case ROW:
650  dofs = problem_ptr->numeredRowDofsPtr;
651  break;
652  case COL:
653  dofs = problem_ptr->numeredColDofsPtr;
654  break;
655  default:
656  SETERRQ(mofemComm, MOFEM_DATA_INCONSISTENCY,
657  "It works only with rows or columns");
658  }
659 
660  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
661  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
662  method.fieldPtr = *field_it;
663  } else {
664  SETERRQ(mofemComm, MOFEM_NOT_FOUND,
665  ("Field not found " + field_name).c_str());
666  }
667 
668  auto miit = dofs->lower_bound(
669  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
670  auto hi_miit = dofs->upper_bound(
671  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
672 
673  using FieldEntity_view_multiIndex = multi_index_container<
674 
675  boost::shared_ptr<FieldEntity>,
676  indexed_by<
677 
678  ordered_unique<
679 
680  tag<Ent_mi_tag>,
681  const_mem_fun<FieldEntity::interface_type_RefEntity, EntityHandle,
683 
684  >>;
685 
686  FieldEntity_view_multiIndex ents_view;
687  auto hint = ents_view.begin();
688  for (; miit != hi_miit; ++miit)
689  if ((*miit)->getPart() >= lower_rank && (*miit)->getPart() <= upper_rank)
690  ents_view.emplace_hint(hint, (*miit)->getFieldEntityPtr());
691 
692  method.loopSize = ents_view.size();
693  SET_BASIC_METHOD(method, problem_ptr);
694  CHKERR method.preProcess();
695  method.nInTheLoop = 0;
696  for (auto &field_ent : ents_view) {
697  method.entPtr = field_ent;
698  CHKERR method();
699  ++method.nInTheLoop;
700  }
701  CHKERR method.postProcess();
703 }
704 
705 MoFEMErrorCode Core::loop_entities(const std::string problem_name,
706  const std::string field_name, RowColData rc,
707  EntityMethod &method, int lower_rank,
708  int upper_rank, int verb) {
710  if (verb == DEFAULT_VERBOSITY)
711  verb = verbose;
712  // find p_miit
713  auto &prb = pRoblems.get<Problem_mi_tag>();
714  auto p_miit = prb.find(problem_name);
715  if (p_miit == prb.end())
716  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
717  problem_name.c_str());
718  CHKERR loop_entities(&*p_miit, field_name, rc, method, lower_rank, upper_rank,
719  verb);
721 }
722 
723 MoFEMErrorCode Core::loop_entities(const std::string problem_name,
724  const std::string field_name, RowColData rc,
725  EntityMethod &method, int verb) {
726  return loop_entities(problem_name, field_name, rc, method, rAnk, rAnk, verb);
727 }
728 
729 MoFEMErrorCode Core::loop_entities(const std::string field_name,
730  EntityMethod &method,
731  Range const *const ents, int verb) {
733  if (verb == DEFAULT_VERBOSITY)
734  verb = verbose;
735  SET_BASIC_METHOD(method, nullptr);
736 
737  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
738  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
739  method.fieldPtr = *field_it;
740  } else {
741  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found %s",
742  field_name.c_str());
743  }
744 
745  auto lo_eit = entsFields.get<Unique_mi_tag>().lower_bound(
746  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
747  auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
748  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
749 
750  typedef multi_index_container<
751  boost::shared_ptr<FieldEntity>,
752  indexed_by<ordered_unique<
753  tag<Ent_mi_tag>, const_mem_fun<FieldEntity::interface_RefEntity,
754  EntityHandle, &FieldEntity::getEnt>>>>
755  FieldEntity_view_multiIndex;
756 
757  FieldEntity_view_multiIndex ents_view;
758  ents_view.insert(lo_eit, hi_eit);
759 
760  method.loopSize = ents_view.size();
761  CHKERR method.preProcess();
762  method.nInTheLoop = 0;
763 
764  if (ents)
765  for (auto p = ents->const_pair_begin(); p != ents->const_pair_end(); ++p)
766  for (auto feit = ents_view.lower_bound(p->first);
767  feit != ents_view.upper_bound(p->second); ++feit) {
768  method.entPtr = *feit;
769  CHKERR method();
770  ++method.nInTheLoop;
771  }
772  else
773  for (auto &field_ent : ents_view) {
774  method.entPtr = field_ent;
775  CHKERR method();
776  ++method.nInTheLoop;
777  }
778 
779  CHKERR method.postProcess();
781 }
782 
783 MoFEMErrorCode Core::cache_problem_entities(const std::string prb_name,
784  CacheTupleWeakPtr cache_weak_ptr) {
786 
787  if (auto cache_ptr = cache_weak_ptr.lock()) {
788  auto p_miit = pRoblems.get<Problem_mi_tag>().find(prb_name);
789  if (p_miit == pRoblems.get<Problem_mi_tag>().end())
790  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
791  prb_name.c_str());
792 
793  const BitRefLevel prb_bit = p_miit->getBitRefLevel();
794  const BitRefLevel prb_mask = p_miit->getMaskBitRefLevel();
795  const BitFEId prb_fe_id = p_miit->getBitFEId();
796 
797  auto &row_dofs = p_miit->numeredRowDofsPtr;
798  auto &col_dofs = p_miit->numeredColDofsPtr;
799 
800  auto &cache_data = std::get<0>(*cache_ptr);
801  auto &cache_row = std::get<1>(*cache_ptr);
802  auto &cache_col = std::get<2>(*cache_ptr);
803 
804  cache_row.resize(entsFields.size());
805  if (row_dofs != col_dofs)
806  cache_col.resize(entsFields.size());
807  cache_data.resize(entsFields.size());
808 
809  size_t idx = 0;
810  for (auto it = entsFields.begin(); it != entsFields.end(); ++it, ++idx) {
811 
812  const auto uid = (*it)->getLocalUniqueId();
813  auto r = entFEAdjacencies.get<Unique_mi_tag>().equal_range(uid);
814  for (auto lo = r.first; lo != r.second; ++lo) {
815 
816  if ((lo->getBitFEId() & prb_fe_id).any()) {
817 
818  const BitRefLevel fe_bit = lo->entFePtr->getBitRefLevel();
819 
820  // if entity is not problem refinement level
821  if (((fe_bit & prb_mask) != fe_bit) ||
822  ((fe_bit & prb_bit) != prb_bit))
823  continue;
824 
825  auto cache_numered_dofs = [&](auto &numered_dofs, auto &cache_vec,
826  auto &ent_cache) {
827  auto dit = numered_dofs->lower_bound(uid);
828  decltype(dit) hi_dit;
829  if (dit != numered_dofs->end())
830  hi_dit = numered_dofs->upper_bound(
831  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
832  else
833  hi_dit = dit;
834 
835  ent_cache = boost::shared_ptr<EntityCacheNumeredDofs>(
836  cache_ptr, &(cache_vec[idx]));
837  cache_vec[idx].loHi = {dit, hi_dit};
838  };
839 
840  auto cache_dofs = [&](auto &dofs, auto &cache_vec, auto &ent_cache) {
841  auto dit = dofs.lower_bound(uid);
842  decltype(dit) hi_dit;
843  if (dit != dofs.end())
844  hi_dit = 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<EntityCacheDofs>(cache_ptr,
850  &(cache_vec[idx]));
851  cache_vec[idx].loHi = {dit, hi_dit};
852  };
853 
854  cache_numered_dofs(row_dofs, cache_row, (*it)->entityCacheRowDofs);
855  if (row_dofs != col_dofs) {
856  if (cache_col.size() != entsFields.size())
857  cache_col.resize(entsFields.size());
858  cache_numered_dofs(col_dofs, cache_col, (*it)->entityCacheColDofs);
859  } else {
860  (*it)->entityCacheColDofs = (*it)->entityCacheRowDofs;
861  }
862 
863  cache_dofs(dofsField, cache_data, (*it)->entityCacheDataDofs);
864 
865  break;
866  }
867  }
868  }
869  } else {
870  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cache not allocated");
871  }
872 
874 }
875 
876 } // namespace MoFEM
static Index< 'p', 3 > p
#define ProblemCoreFunctionBegin
#define SET_BASIC_METHOD(METHOD, PROBLEM_PTR)
@ DEFAULT_VERBOSITY
Definition: definitions.h:282
RowColData
RowColData.
Definition: definitions.h:198
@ COL
Definition: definitions.h:198
@ ROW
Definition: definitions.h:198
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:190
@ MF_EXIST
Definition: definitions.h:193
@ MF_EXCL
Definition: definitions.h:192
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:311
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:522
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:421
@ MOFEM_NOT_FOUND
Definition: definitions.h:126
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:127
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:124
@ MOFEM_INVALID_DATA
Definition: definitions.h:129
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:491
#define CHKERR
Inline error check.
Definition: definitions.h:610
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:515
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:636
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:300
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:329
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:54
std::bitset< BITPROBLEMID_SIZE > BitProblemId
Problem Id.
Definition: Types.hpp:55
uint128_t UId
Unique Id.
Definition: Types.hpp:42
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
boost::shared_ptr< CacheTuple > CacheTupleSharedPtr
boost::weak_ptr< CacheTuple > CacheTupleWeakPtr
const double r
rate factor
Data structure to exchange data between mofem and User Loop Methods.
int loopSize
local number oe methods to process
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
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)
build problem data structures, assuming that mesh is distributed (collective)
bool check_problem(const std::string name)
check if problem exist
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 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
BitProblemId getBitProblemId(const std::string &name) const
DEPRECATED MoFEMErrorCode build_problems(int verb=DEFAULT_VERBOSITY)
build problem data structures
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.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
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
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
std::string getName() const
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< REFENT > &sptr)
EntityHandle getEnt() const
Get the entity handle.