v0.10.0
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)
77  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Problem not added");
78 
79  MOFEM_LOG("WORLD", Sev::inform) << "Add problem " << name;
80 
82 }
83 
84 MoFEMErrorCode 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  BitProblemId id = getProblemShift();
92  CHKERR addProblem(id, name, verb);
93  } else if (bh == MF_EXCL) {
94  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem is in database %s",
95  name.c_str());
96  }
98 }
99 
100 MoFEMErrorCode Core::delete_problem(const std::string name) {
102  auto p_miit = pRoblems.get<Problem_mi_tag>().find(name);
103  if (p_miit == pRoblems.get<Problem_mi_tag>().end()) {
104  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "no such problem like < %s >",
105  name.c_str());
106  }
107  const EntityHandle meshset = p_miit->meshset;
108  pRoblems.get<Problem_mi_tag>().erase(p_miit);
109  CHKERR get_moab().delete_entities(&meshset, 1);
111 }
112 
113 BitProblemId Core::getBitProblemId(const std::string &name) const {
114  auto p_miit = pRoblems.get<Problem_mi_tag>().find(name);
115  if (p_miit == pRoblems.get<Problem_mi_tag>().end()) {
116  THROW_MESSAGE("no such problem like " + name + " >");
117  }
118  return p_miit->getId();
119 }
120 
124  const ProblemById &set_id = pRoblems.get<BitProblemId_mi_tag>();
125  ProblemById::iterator miit = set_id.begin();
126  for (; miit != set_id.end(); miit++) {
127  std::ostringstream ss;
128  ss << *miit << std::endl;
129  PetscPrintf(cOmm, ss.str().c_str());
130  }
132 }
133 
135 Core::modify_problem_add_finite_element(const std::string &name_problem,
136  const std::string &fe_name) {
139  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
140  ProblemsByName::iterator miit = set.find(name_problem);
141  if (miit == set.end()) {
142  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is not there",
143  name_problem.c_str());
144  }
145  BitFEId f_id = getBitFEId(fe_name);
146  bool success = set.modify(miit, ProblemFiniteElementChangeBitAdd(f_id));
147  if (!success)
148  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
149  "modification unsuccessful");
151 }
152 
154 Core::modify_problem_unset_finite_element(const std::string &name_problem,
155  const std::string &fe_name) {
158  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
159  ProblemsByName::iterator miit = set.find(name_problem);
160  if (miit == set.end()) {
161  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is not there",
162  name_problem.c_str());
163  }
164  BitFEId f_id = getBitFEId(fe_name);
165  bool success = set.modify(miit, ProblemFiniteElementChangeBitUnSet(f_id));
166  if (!success)
167  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
168  "modification unsuccessful");
170 }
171 
173 Core::modify_problem_ref_level_add_bit(const std::string &name_problem,
174  const BitRefLevel &bit) {
177  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
178  ProblemsByName::iterator miit = set.find(name_problem);
179  std::ostringstream ss;
180  ss << name_problem;
181  if (miit == set.end())
182  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
183  ss.str().c_str());
184  bool success = set.modify(miit, ProblemChangeRefLevelBitAdd(bit));
185  if (!success)
186  SETERRQ(PETSC_COMM_SELF, 1, "modification unsuccessful");
188 }
189 
191 Core::modify_problem_ref_level_set_bit(const std::string &name_problem,
192  const BitRefLevel &bit) {
195  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
196  ProblemsByName::iterator miit = set.find(name_problem);
197  std::ostringstream ss;
198  ss << name_problem;
199  if (miit == set.end())
200  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
201  ss.str().c_str());
202  bool success = set.modify(miit, ProblemChangeRefLevelBitSet(bit));
203  if (!success)
204  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
205  "modification unsuccessful");
207 }
208 
210 Core::modify_problem_mask_ref_level_add_bit(const std::string &name_problem,
211  const BitRefLevel &bit) {
214  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
215  ProblemsByName::iterator miit = set.find(name_problem);
216  if (miit == set.end()) {
217  std::ostringstream ss;
218  ss << name_problem;
219  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
220  ss.str().c_str());
221  }
222  bool success = set.modify(miit, ProblemChangeRefLevelBitDofMaskAdd(bit));
223  if (!success)
224  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
225  "modification unsuccessful");
227 }
228 
230 Core::modify_problem_mask_ref_level_set_bit(const std::string &name_problem,
231  const BitRefLevel &bit) {
233  ProblemsByName;
234  ProblemsByName &set = pRoblems.get<Problem_mi_tag>();
235  ProblemsByName::iterator miit = set.find(name_problem);
236  if (miit == set.end()) {
237  std::ostringstream ss;
238  ss << name_problem;
239  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "this problem <%s> is there",
240  ss.str().c_str());
241  }
242  bool success = set.modify(miit, ProblemChangeRefLevelBitDofMaskSet(bit));
243  if (!success)
244  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
245  "modification unsuccessful");
247 }
248 
251  if (verb == -1)
252  verb = verbose;
253  Problem_multiIndex::iterator p_miit = pRoblems.begin();
254  for (; p_miit != pRoblems.end(); p_miit++) {
255  CHKERR getInterface<ProblemsManager>()->buildProblemOnDistributedMesh(
256  const_cast<Problem *>(&*p_miit), verb);
257  }
259 }
260 
261 MoFEMErrorCode Core::clear_problem(const std::string problem_name, int verb) {
263  if (verb == -1)
264  verb = verbose;
266  ProblemsByName &prob_by_name = pRoblems.get<Problem_mi_tag>();
267  ProblemsByName::iterator p_miit = prob_by_name.find(problem_name);
268  if (p_miit == prob_by_name.end()) {
269  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
270  "problem < %s > not found, (top tip: check spelling)",
271  problem_name.c_str());
272  }
273  // zero rows
274  bool success = prob_by_name.modify(p_miit, ProblemZeroNbRowsChange());
275  if (!success)
276  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
277  "modification unsuccessful");
278  // zero cols
279  success = prob_by_name.modify(p_miit, ProblemZeroNbColsChange());
280  if (!success)
281  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
282  "modification unsuccessful");
283  // clear finite elements
284  success =
285  prob_by_name.modify(p_miit, ProblemClearNumeredFiniteElementsChange());
286  if (!success)
287  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
288  "modification unsuccessful");
289  // clear data structures
290  success = prob_by_name.modify(p_miit, ProblemClearSubProblemData());
291  if (!success)
292  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
293  "modification unsuccessful");
294  success = prob_by_name.modify(p_miit, ProblemClearComposedProblemData());
295  if (!success)
296  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
297  "modification unsuccessful");
298  if (p_miit->getRowDofsSequence())
299  p_miit->getRowDofsSequence()->clear();
300  if (p_miit->getColDofsSequence())
301  p_miit->getColDofsSequence()->clear();
302  if (p_miit->getSubData())
303  p_miit->getSubData().reset();
304  if (p_miit->getComposedProblemsData())
305  p_miit->getComposedProblemsData().reset();
307 }
308 
311  if (verb == -1)
312  verb = verbose;
313  if (!((*buildMoFEM) & BUILD_FIELD))
314  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "fields not build");
315  if (!((*buildMoFEM) & BUILD_FE))
316  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "FEs not build");
317  if (!((*buildMoFEM) & BUILD_ADJ))
318  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "adjacencies not build");
319  // iterate problems
320  Problem_multiIndex::iterator p_miit = pRoblems.begin();
321  for (; p_miit != pRoblems.end(); p_miit++) {
322  Problem *problem_ptr = const_cast<Problem *>(&*p_miit);
323  CHKERR getInterface<ProblemsManager>()->buildProblem(problem_ptr, false,
324  verb);
325  }
326  *buildMoFEM |= BUILD_PROBLEM;
328 }
329 
332  if (verb == -1)
333  verb = verbose;
334  // iterate problems
335  for (auto p_miit = pRoblems.begin(); p_miit != pRoblems.end(); p_miit++)
336  CHKERR clear_problem(p_miit->getName(), verb);
338 }
339 
340 #define SET_BASIC_METHOD(METHOD, PROBLEM_PTR) \
341  { \
342  METHOD.rAnk = rAnk; \
343  METHOD.sIze = sIze; \
344  METHOD.problemPtr = PROBLEM_PTR; \
345  METHOD.fieldsPtr = &fIelds; \
346  METHOD.refinedEntitiesPtr = &refinedEntities; \
347  METHOD.entitiesPtr = &entsFields; \
348  METHOD.dofsPtr = &dofsField; \
349  METHOD.refinedFiniteElementsPtr = &refinedFiniteElements; \
350  METHOD.finiteElementsPtr = &finiteElements; \
351  METHOD.finiteElementsEntitiesPtr = &entsFiniteElements; \
352  METHOD.adjacenciesPtr = &entFEAdjacencies; \
353  }
354 
356  BasicMethod &method,
357  int verb) {
359  if (verb == -1)
360  verb = verbose;
361  // finite element
362  SET_BASIC_METHOD(method, problem_ptr)
363  PetscLogEventBegin(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
364  CHKERR method.preProcess();
365  PetscLogEventEnd(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
367 }
368 
370 Core::problem_basic_method_preProcess(const std::string &problem_name,
371  BasicMethod &method, int verb) {
373  if (verb == -1)
374  verb = verbose;
376  // find p_miit
377  ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
378  ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
379  if (p_miit == pRoblems_set.end())
380  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem is not in database %s",
381  problem_name.c_str());
382  CHKERR problem_basic_method_preProcess(&*p_miit, method, verb);
384 }
385 
388  BasicMethod &method, int verb) {
390  SET_BASIC_METHOD(method, problem_ptr)
391 
392  PetscLogEventBegin(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
393  CHKERR method.postProcess();
394  PetscLogEventEnd(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
395 
397 }
398 
400 Core::problem_basic_method_postProcess(const std::string &problem_name,
401  BasicMethod &method, int verb) {
403  if (verb == -1)
404  verb = verbose;
406 
407  // find p_miit
408  ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
409  ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
410  if (p_miit == pRoblems_set.end())
411  SETERRQ1(cOmm, 1, "problem is not in database %s", problem_name.c_str());
412 
413  CHKERR problem_basic_method_postProcess(&*p_miit, method, verb);
414 
416 }
417 
419  const Problem *problem_ptr, const std::string &fe_name, FEMethod &method,
420  int lower_rank, int upper_rank,
421  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
422  CacheTupleWeakPtr cache_ptr, int verb) {
424  if (verb == DEFAULT_VERBOSITY)
425  verb = verbose;
426 
427  CacheTupleSharedPtr tmp_cache_ptr;
428  if(!cache_ptr.use_count()) {
429  tmp_cache_ptr = boost::make_shared<CacheTuple>();
430  CHKERR cache_problem_entities(problem_ptr->getName(), tmp_cache_ptr);
431  }
432 
433  method.feName = fe_name;
434  SET_BASIC_METHOD(method, &*problem_ptr)
435  PetscLogEventBegin(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
436  CHKERR method.preProcess();
437  PetscLogEventEnd(MOFEM_EVENT_preProcess, 0, 0, 0, 0);
438 
439  if (!fe_ptr)
440  fe_ptr = problem_ptr->numeredFiniteElementsPtr;
441 
442  auto miit = fe_ptr->get<Composite_Name_And_Part_mi_tag>().lower_bound(
443  boost::make_tuple(fe_name, lower_rank));
444  auto hi_miit = fe_ptr->get<Composite_Name_And_Part_mi_tag>().upper_bound(
445  boost::make_tuple(fe_name, upper_rank));
446 
447  if (miit == hi_miit && (bh & MF_EXIST)) {
448  if (!check_finite_element(fe_name)) {
449  SETERRQ1(cOmm, MOFEM_NOT_FOUND, "finite element < %s > not found",
450  fe_name.c_str());
451  }
452  }
453 
454  method.loopSize = std::distance(miit, hi_miit);
455  for (int nn = 0; miit != hi_miit; miit++, nn++) {
456 
457  method.nInTheLoop = nn;
458  method.numeredEntFiniteElementPtr = *miit;
459 
460  PetscLogEventBegin(MOFEM_EVENT_operator, 0, 0, 0, 0);
461  CHKERR method();
462  PetscLogEventEnd(MOFEM_EVENT_operator, 0, 0, 0, 0);
463  }
464 
465  PetscLogEventBegin(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
466  CHKERR method.postProcess();
467  PetscLogEventEnd(MOFEM_EVENT_postProcess, 0, 0, 0, 0);
468 
470 }
471 
473  const std::string &problem_name, const std::string &fe_name,
474  FEMethod &method,
475  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
476  CacheTupleWeakPtr cache_ptr, int verb) {
478  if (verb == DEFAULT_VERBOSITY)
479  verb = verbose;
480 
481  CHKERR loop_finite_elements(problem_name, fe_name, method, rAnk, rAnk, fe_ptr,
482  bh, cache_ptr, verb);
483 
485 }
486 
488  const std::string &problem_name, const std::string &fe_name,
489  FEMethod &method, int lower_rank, int upper_rank,
490  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> fe_ptr, MoFEMTypes bh,
491  CacheTupleWeakPtr cache_ptr,
492  int verb) {
494  if (verb == DEFAULT_VERBOSITY)
495  verb = verbose;
496 
497  auto &prb_by_name = pRoblems.get<Problem_mi_tag>();
498  auto p_miit = prb_by_name.find(problem_name);
499  if (p_miit == prb_by_name.end())
500  SETERRQ1(cOmm, MOFEM_INVALID_DATA, "Problem <%s> is not in database",
501  problem_name.c_str());
502 
503  CHKERR loop_finite_elements(&*p_miit, fe_name, method, lower_rank, upper_rank,
504  fe_ptr, bh, cache_ptr, verb);
505 
507 }
508 
510  const std::string &field_name, RowColData rc,
511  DofMethod &method, int lower_rank,
512  int upper_rank, int verb) {
514  SET_BASIC_METHOD(method, &*problem_ptr);
516  NumeredDofsByUId;
517  NumeredDofsByUId *dofs;
518  switch (rc) {
519  case ROW:
520  dofs = &problem_ptr->numeredRowDofsPtr->get<Unique_mi_tag>();
521  break;
522  case COL:
523  dofs = &problem_ptr->numeredColDofsPtr->get<Unique_mi_tag>();
524  break;
525  default:
526  SETERRQ(cOmm, MOFEM_DATA_INCONSISTENCY, "Not implemented");
527  }
528 
529  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
530  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
531  method.fieldPtr = *field_it;
532  } else {
533  SETERRQ(cOmm, MOFEM_NOT_FOUND, ("Field not found " + field_name).c_str());
534  }
535 
536  auto miit = dofs->lower_bound(
537  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
538  auto hi_miit = dofs->upper_bound(
539  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
540 
541  CHKERR method.preProcess();
542  for (; miit != hi_miit; miit++) {
543  if ((*miit)->getPart() >= lower_rank && (*miit)->getPart() <= upper_rank) {
544  method.dofPtr = miit->get()->getDofEntityPtr();
545  method.dofNumeredPtr = *miit;
546  CHKERR method();
547  }
548  }
549  CHKERR method.postProcess();
551 }
552 
554  const std::string &problem_name, const std::string &field_name,
555  RowColData rc, // ROW or COL
556  DofMethod &method, // Finite element instance processed on each DOF
557  int lower_rank, // Only DOFs on processor higher or equal to this are
558  // processed
559  int upper_rank, // Only DOFs lowest or higher to this are processed
560  int verb // verbosity level
561 ) {
563  if (verb == DEFAULT_VERBOSITY)
564  verb = verbose;
566  // find p_miit
567  ProblemsByName &pRoblems_set = pRoblems.get<Problem_mi_tag>();
568  ProblemsByName::iterator p_miit = pRoblems_set.find(problem_name);
569  if (p_miit == pRoblems_set.end())
570  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
571  problem_name.c_str());
572  CHKERR loop_dofs(&*p_miit, field_name, rc, method, lower_rank, upper_rank,
573  verb);
575 }
576 
577 MoFEMErrorCode Core::loop_dofs(const std::string &problem_name,
578  const std::string &field_name, RowColData rc,
579  DofMethod &method, int verb) {
581  if (verb == DEFAULT_VERBOSITY)
582  verb = verbose;
583  CHKERR loop_dofs(problem_name, field_name, rc, method, 0, sIze, verb);
585 }
586 
587 MoFEMErrorCode Core::loop_dofs(const std::string &field_name, DofMethod &method,
588  int verb) {
590  if (verb == DEFAULT_VERBOSITY)
591  verb = verbose;
592  SET_BASIC_METHOD(method, nullptr);
593 
594  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
595  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
596  method.fieldPtr = *field_it;
597  } else {
598  SETERRQ(cOmm, MOFEM_NOT_FOUND, ("Field not found " + field_name).c_str());
599  }
600 
601  auto miit = dofsField.get<Unique_mi_tag>().lower_bound(
602  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
603  auto hi_miit = dofsField.get<Unique_mi_tag>().upper_bound(
604  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
605 
606  method.loopSize = std::distance(miit, hi_miit);
607  CHKERR method.preProcess();
608  for (int nn = 0; miit != hi_miit; miit++, nn++) {
609  method.nInTheLoop = nn;
610  method.dofPtr = *miit;
611  CHKERR method();
612  }
613  CHKERR method.postProcess();
615 }
616 
618  const std::string field_name, RowColData rc,
619  EntityMethod &method, int lower_rank,
620  int upper_rank, int verb) {
622  if (verb == DEFAULT_VERBOSITY)
623  verb = verbose;
624  decltype(problem_ptr->numeredRowDofsPtr) dofs;
625  switch (rc) {
626  case ROW:
627  dofs = problem_ptr->numeredRowDofsPtr;
628  break;
629  case COL:
630  dofs = problem_ptr->numeredColDofsPtr;
631  break;
632  default:
633  SETERRQ(cOmm, MOFEM_DATA_INCONSISTENCY,
634  "It works only with rows or columns");
635  }
636 
637  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
638  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
639  method.fieldPtr = *field_it;
640  } else {
641  SETERRQ(cOmm, MOFEM_NOT_FOUND, ("Field not found " + field_name).c_str());
642  }
643 
644  auto miit = dofs->lower_bound(
645  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
646  auto hi_miit = dofs->upper_bound(
647  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
648 
649  using FieldEntity_view_multiIndex = multi_index_container<
650 
651  boost::shared_ptr<FieldEntity>,
652  indexed_by<
653 
654  ordered_unique<
655 
656  tag<Ent_mi_tag>,
659 
660  >>;
661 
662  FieldEntity_view_multiIndex ents_view;
663  auto hint = ents_view.begin();
664  for (; miit != hi_miit; ++miit)
665  if ((*miit)->getPart() >= lower_rank && (*miit)->getPart() <= upper_rank)
666  ents_view.emplace_hint(hint, (*miit)->getFieldEntityPtr());
667 
668  method.loopSize = ents_view.size();
669  SET_BASIC_METHOD(method, problem_ptr);
670  CHKERR method.preProcess();
671  method.nInTheLoop = 0;
672  for (auto &field_ent : ents_view) {
673  method.entPtr = field_ent;
674  CHKERR method();
675  ++method.nInTheLoop;
676  }
677  CHKERR method.postProcess();
679 }
680 
681 MoFEMErrorCode Core::loop_entities(const std::string problem_name,
682  const std::string field_name, RowColData rc,
683  EntityMethod &method, int lower_rank,
684  int upper_rank, int verb) {
686  if (verb == DEFAULT_VERBOSITY)
687  verb = verbose;
688  // find p_miit
689  auto &prb = pRoblems.get<Problem_mi_tag>();
690  auto p_miit = prb.find(problem_name);
691  if (p_miit == prb.end())
692  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
693  problem_name.c_str());
694  CHKERR loop_entities(&*p_miit, field_name, rc, method, lower_rank, upper_rank,
695  verb);
697 }
698 
699 MoFEMErrorCode Core::loop_entities(const std::string problem_name,
700  const std::string field_name, RowColData rc,
701  EntityMethod &method, int verb) {
702  return loop_entities(problem_name, field_name, rc, method, rAnk, rAnk, verb);
703 }
704 
705 MoFEMErrorCode Core::loop_entities(const std::string field_name,
706  EntityMethod &method,
707  Range const *const ents, int verb) {
709  if (verb == DEFAULT_VERBOSITY)
710  verb = verbose;
711  SET_BASIC_METHOD(method, nullptr);
712 
713  auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
714  if (field_it != fIelds.get<FieldName_mi_tag>().end()) {
715  method.fieldPtr = *field_it;
716  } else {
717  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found %s",
718  field_name.c_str());
719  }
720 
721  auto lo_eit = entsFields.get<Unique_mi_tag>().lower_bound(
722  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
723  auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
724  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
725 
726  typedef multi_index_container<
727  boost::shared_ptr<FieldEntity>,
728  indexed_by<ordered_unique<
729  tag<Ent_mi_tag>, const_mem_fun<FieldEntity::interface_RefEntity,
731  FieldEntity_view_multiIndex;
732 
733  FieldEntity_view_multiIndex ents_view;
734  ents_view.insert(lo_eit, hi_eit);
735 
736  method.loopSize = ents_view.size();
737  CHKERR method.preProcess();
738  method.nInTheLoop = 0;
739 
740  if (ents)
741  for (auto p = ents->const_pair_begin(); p != ents->const_pair_end(); ++p)
742  for (auto feit = ents_view.lower_bound(p->first);
743  feit != ents_view.upper_bound(p->second); ++feit) {
744  method.entPtr = *feit;
745  CHKERR method();
746  ++method.nInTheLoop;
747  }
748  else
749  for (auto &field_ent : ents_view) {
750  method.entPtr = field_ent;
751  CHKERR method();
752  ++method.nInTheLoop;
753  }
754 
755  CHKERR method.postProcess();
757 }
758 
759 MoFEMErrorCode Core::cache_problem_entities(const std::string prb_name,
760  CacheTupleWeakPtr cache_weak_ptr) {
762 
763  if (auto cache_ptr = cache_weak_ptr.lock()) {
764  auto p_miit = pRoblems.get<Problem_mi_tag>().find(prb_name);
765  if (p_miit == pRoblems.get<Problem_mi_tag>().end())
766  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "problem not in database %s",
767  prb_name.c_str());
768 
769  const BitRefLevel prb_bit = p_miit->getBitRefLevel();
770  const BitRefLevel prb_mask = p_miit->getMaskBitRefLevel();
771  const BitFEId prb_fe_id = p_miit->getBitFEId();
772 
773  auto &row_dofs = p_miit->numeredRowDofsPtr;
774  auto &col_dofs = p_miit->numeredColDofsPtr;
775 
776  auto &cache_data = std::get<0>(*cache_ptr);
777  auto &cache_row = std::get<1>(*cache_ptr);
778  auto &cache_col = std::get<2>(*cache_ptr);
779 
780  cache_row.resize(entsFields.size());
781  if (row_dofs != col_dofs)
782  cache_col.resize(entsFields.size());
783  cache_data.resize(entsFields.size());
784 
785  size_t idx = 0;
786  for (auto it = entsFields.begin(); it != entsFields.end(); ++it, ++idx) {
787 
788  const auto uid = (*it)->getLocalUniqueId();
789  auto r = entFEAdjacencies.get<Unique_mi_tag>().equal_range(uid);
790  for (auto lo = r.first; lo != r.second; ++lo) {
791 
792  if ((lo->getBitFEId() & prb_fe_id).any()) {
793 
794  const BitRefLevel fe_bit = lo->entFePtr->getBitRefLevel();
795 
796  // if entity is not problem refinement level
797  if (((fe_bit & prb_mask) != fe_bit) ||
798  ((fe_bit & prb_bit) != prb_bit))
799  continue;
800 
801  auto cache_numered_dofs = [&](auto &numered_dofs, auto &cache_vec,
802  auto &ent_cache) {
803  auto dit = numered_dofs->lower_bound(uid);
804  decltype(dit) hi_dit;
805  if (dit != numered_dofs->end())
806  hi_dit = numered_dofs->upper_bound(
807  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
808  else
809  hi_dit = dit;
810 
811  ent_cache = boost::shared_ptr<EntityCacheNumeredDofs>(
812  cache_ptr, &(cache_vec[idx]));
813  cache_vec[idx].loHi = {dit, hi_dit};
814  };
815 
816  auto cache_dofs = [&](auto &dofs, auto &cache_vec, auto &ent_cache) {
817  auto dit = dofs.lower_bound(uid);
818  decltype(dit) hi_dit;
819  if (dit != dofs.end())
820  hi_dit = dofs.upper_bound(
821  uid | static_cast<UId>(MAX_DOFS_ON_ENTITY - 1));
822  else
823  hi_dit = dit;
824 
825  ent_cache = boost::shared_ptr<EntityCacheDofs>(cache_ptr,
826  &(cache_vec[idx]));
827  cache_vec[idx].loHi = {dit, hi_dit};
828  };
829 
830  cache_numered_dofs(row_dofs, cache_row, (*it)->entityCacheRowDofs);
831  if (row_dofs != col_dofs) {
832  if (cache_col.size() != entsFields.size())
833  cache_col.resize(entsFields.size());
834  cache_numered_dofs(col_dofs, cache_col, (*it)->entityCacheColDofs);
835  } else {
836  (*it)->entityCacheColDofs = (*it)->entityCacheRowDofs;
837  }
838 
839  cache_dofs(dofsField, cache_data, (*it)->entityCacheDataDofs);
840 
841  break;
842  }
843  }
844  }
845  } else {
846  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cache not allocated");
847  }
848 
850 }
851 
852 } // namespace MoFEM
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
boost::shared_ptr< DofEntity > dofPtr
Clear sub-problem data structure.
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)
Cache variables.
MoFEMErrorCode clear_problems(int verb=DEFAULT_VERBOSITY)
clear problems
MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)
Set data for BasicMethodThis function set data about problem, adjacencies and other multi-indices in ...
boost::shared_ptr< Field > fieldPtr
MoFEMErrorCode modify_problem_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)
set ref level for problem
int loopSize
local number oe methods to process
Data structure to exchange data between mofem and User Loop Methods on entities.It allows to exchange...
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
boost::shared_ptr< CacheTuple > CacheTupleSharedPtr
static Index< 'p', 3 > p
boost::weak_ptr< CacheTuple > CacheTupleWeakPtr
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
std::string feName
Name of finite element.
boost::shared_ptr< Field > fieldPtr
#define THROW_MESSAGE(a)
Throw MoFEM exception.
Definition: definitions.h:628
std::bitset< BITPROBLEMID_SIZE > BitProblemId
Problem Id.
Definition: Types.hpp:55
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.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
MoFEMErrorCode modify_problem_mask_ref_level_set_bit(const std::string &name_problem, const BitRefLevel &bit)
set dof mask ref level for problem
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
Data structure to exchange data between mofem and User Loop Methods on entities.It allows to exchange...
keeps basic data about problemThis is low level structure with information about problem,...
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
#define SET_BASIC_METHOD(METHOD, PROBLEM_PTR)
RowColData
RowColData.
Definition: definitions.h:192
MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)
Set data for BasicMethod.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
BitProblemId getBitProblemId(const std::string &name) const
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
Add problem.
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:54
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode delete_problem(const std::string name)
Delete problem.
DEPRECATED MoFEMErrorCode build_problems(int verb=DEFAULT_VERBOSITY)
build problem data structures
remove finite element from problem
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
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.
DEPRECATED MoFEMErrorCode build_problem_on_distributed_mesh(int verb=DEFAULT_VERBOSITY)
build problem data structures, assuming that mesh is distributed (collective)
interface_RefEntity(const boost::shared_ptr< RefEntity > &sptr)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEMErrorCode clear_problem(const std::string name, int verb=DEFAULT_VERBOSITY)
clear problem
#define ProblemCoreFunctionBegin
MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)
add ref level to problem
MoFEMErrorCode addProblem(const BitProblemId id, const std::string &name, int verb=DEFAULT_VERBOSITY)
add problem
interface_RefEntity< RefEntity > interface_type_RefEntity
MoFEMErrorCode list_problem() const
list problems
bool check_problem(const std::string name)
check if problem exist
#define CHKERR
Inline error check.
Definition: definitions.h:604
const double r
rate factor
int nInTheLoop
number currently of processed method
boost::shared_ptr< NumeredDofEntity > dofNumeredPtr
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
boost::shared_ptr< FieldEntity > entPtr
MultiIndex Tag for field name.
EntityHandle getEnt() const
Get the entity handle.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
MoFEMErrorCode modify_problem_mask_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)
set dof mask ref level for problem
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:305
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
std::string getName() const
Clear composed problem data structure.
Data structure to exchange data between mofem and User Loop Methods.It allows to exchange data betwee...
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:189
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
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 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...