v0.10.0
ForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file ForcesAndSourcesCore.cpp
2 
3 \brief Implementation of Elements on Entities for Forces and Sources
4 */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 #ifdef __cplusplus
21 extern "C" {
22 #endif
23 #include <cblas.h>
24 #include <lapack_wrap.h>
25 // #include <gm_rule.h>
26 #ifdef __cplusplus
27 }
28 #endif
29 
30 namespace MoFEM {
31 
32 static auto cmp_uid_lo(const boost::weak_ptr<FieldEntity> &a, const UId &b) {
33  if (auto a_ptr = a.lock()) {
34  if (a_ptr->getLocalUniqueId() < b)
35  return true;
36  else
37  return false;
38  } else {
39  return false;
40  }
41 }
42 
43 static auto cmp_uid_hi(const UId &b, const boost::weak_ptr<FieldEntity> &a) {
44  if (auto a_ptr = a.lock()) {
45  if (b < a_ptr->getLocalUniqueId())
46  return true;
47  else
48  return false;
49  } else {
50  return true;
51  }
52 }
53 
55  :
56 
57  mField(m_field), getRuleHook(0), setRuleHook(0),
58  dataOnElement{
59 
60  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // NOSPACE,
61  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // NOFIELD
62  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // H1
63  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // HCURL
64  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET), // HDIV
65  boost::make_shared<DataForcesAndSourcesCore>(MBENTITYSET) // L2
66 
67  },
68  derivedDataOnElement{
69 
70  nullptr,
71  boost::make_shared<DerivedDataForcesAndSourcesCore>(
72  dataOnElement[NOFIELD]), // NOFIELD
73  boost::make_shared<DerivedDataForcesAndSourcesCore>(
74  dataOnElement[H1]), // H1
75  boost::make_shared<DerivedDataForcesAndSourcesCore>(
76  dataOnElement[HCURL]), // HCURL
77  boost::make_shared<DerivedDataForcesAndSourcesCore>(
78  dataOnElement[HDIV]), // HDIV
79  boost::make_shared<DerivedDataForcesAndSourcesCore>(
80  dataOnElement[L2]) // L2
81 
82  },
83  dataNoField(*dataOnElement[NOFIELD].get()),
84  dataH1(*dataOnElement[H1].get()), dataHcurl(*dataOnElement[HCURL].get()),
85  dataHdiv(*dataOnElement[HDIV].get()), dataL2(*dataOnElement[L2].get()),
86  lastEvaluatedElementEntityType(MBMAXTYPE), sidePtrFE(nullptr) {}
87 
88 // ** Sense **
89 
91  const EntityType type,
92  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const {
94 
95  auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
96  auto sit = side_table.lower_bound(get_id_for_min_type(type));
97  if (sit != side_table.end()) {
98  auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
99  for (; sit != hi_sit; ++sit) {
100  const int side_number = (*sit)->side_number;
101  if (side_number >= 0) {
102  const int brother_side_number = (*sit)->brother_side_number;
103  const int sense = (*sit)->sense;
104 
105  data[side_number].getSense() = sense;
106  if (brother_side_number != -1)
107  data[brother_side_number].getSense() = sense;
108  }
109  }
110  }
112 }
113 
114 // ** Order **
115 
116 template <typename ENTMULTIINDEX>
117 static inline int getMaxOrder(const ENTMULTIINDEX &multi_index) {
118  int max_order = 0;
119  for (auto ent_field_weak_ptr : multi_index)
120  if (auto e = ent_field_weak_ptr.lock()) {
121  const int order = e->getMaxOrder();
122  max_order = (max_order < order) ? order : max_order;
123  }
124  return max_order;
125 }
126 
128  int max_order = 0;
129  for (auto e : getDataFieldEnts()) {
130  if (auto ptr = e.lock()) {
131  const int order = ptr->getMaxOrder();
132  max_order = (max_order < order) ? order : max_order;
133  }
134  }
135  return max_order;
136 }
137 
139  return getMaxOrder(getRowFieldEnts());
140 }
141 
143  return getMaxOrder(getColFieldEnts());
144 }
145 
147  const EntityType type, const FieldSpace space,
148  boost::ptr_vector<DataForcesAndSourcesCore::EntData> &data) const {
150 
151  auto set_order = [&]() {
153  auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable();
154 
155  for (unsigned int s = 0; s != data.size(); ++s)
156  data[s].getDataOrder() = 0;
157 
158  const FieldEntity_vector_view &data_field_ent = getDataFieldEnts();
159 
160  for (auto r = fieldsPtr->get<BitFieldId_space_mi_tag>().equal_range(space);
161  r.first != r.second; ++r.first) {
162 
163  const auto field_bit_number = (*r.first)->getBitNumber();
164  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
165  field_bit_number, get_id_for_min_type(type));
166  auto lo = std::lower_bound(data_field_ent.begin(), data_field_ent.end(),
167  lo_uid, cmp_uid_lo);
168  if (lo != data_field_ent.end()) {
169  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
170  field_bit_number, get_id_for_max_type(type));
171  auto hi =
172  std::upper_bound(lo, data_field_ent.end(), hi_uid, cmp_uid_hi);
173  for (; lo != hi; ++lo) {
174 
175  if (auto ptr = lo->lock()) {
176 
177  auto &e = *ptr;
178  auto sit = side_table.find(e.getEnt());
179  if (sit != side_table.end()) {
180  auto &side = *sit;
181  const int side_number = side->side_number;
182  if (side_number >= 0) {
183  ApproximationOrder ent_order = e.getMaxOrder();
184  auto &dat = data[side_number];
185  dat.getDataOrder() = dat.getDataOrder() > ent_order
186  ? dat.getDataOrder()
187  : ent_order;
188  }
189  } else
190  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
191  "Entity on side of the element not found");
192  }
193  }
194  }
195  }
196 
198  };
199 
200  auto set_order_on_brother = [&]() {
202  auto &side_table =
203  numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
204  auto sit = side_table.lower_bound(get_id_for_min_type(type));
205  if (sit != side_table.end()) {
206  auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
207  for (; sit != hi_sit; ++sit) {
208  const int brother_side_number = (*sit)->brother_side_number;
209  if (brother_side_number != -1) {
210  const int side_number = (*sit)->side_number;
211  data[brother_side_number].getDataOrder() =
212  data[side_number].getDataOrder();
213  }
214  }
215  }
217  };
218 
219  CHKERR set_order();
220  CHKERR set_order_on_brother();
221 
223 }
224 
225 // ** Indices **
226 
227 template <typename EXTRACTOR>
229  const std::string &field_name, FieldEntity_vector_view &ents_field,
230  VectorInt &nodes_indices, VectorInt &local_nodes_indices,
231  EXTRACTOR &&extractor) const {
233 
234  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
235  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
236 
237  auto bit_number = (*field_it)->getBitNumber();
238  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
239  bit_number, get_id_for_min_type<MBVERTEX>());
240  auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
241  cmp_uid_lo);
242  if (lo != ents_field.end()) {
243  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
244  bit_number, get_id_for_max_type<MBVERTEX>());
245  auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
246 
247  int num_nodes;
248  CHKERR getNumberOfNodes(num_nodes);
249  const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
250  const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
251 
252  int nb_dofs = 0;
253  for (auto it = lo; it != hi; ++it) {
254  if (auto e = it->lock()) {
255  if (auto cache = extractor(e).lock()) {
256  if (cache->loHi[0] != cache->loHi[1]) {
257  nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
258  break;
259  }
260  }
261  }
262  }
263 
264  if (nb_dofs) {
265  nodes_indices.resize(max_nb_dofs, false);
266  local_nodes_indices.resize(max_nb_dofs, false);
267  } else {
268  nodes_indices.resize(0, false);
269  local_nodes_indices.resize(0, false);
270  }
271 
272  if (nb_dofs != max_nb_dofs) {
273  std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
274  std::fill(local_nodes_indices.begin(), local_nodes_indices.end(), -1);
275  }
276 
277  for (auto it = lo; it != hi; ++it) {
278  if (auto e = it->lock()) {
279  auto side_ptr = e->getSideNumberPtr();
280  const auto side_number = side_ptr->side_number;
281  if (side_number >= 0) {
282  const auto brother_side_number = side_ptr->brother_side_number;
283  if (auto cache = extractor(e).lock()) {
284  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
285  auto &dof = **dit;
286  const int idx = dof.getPetscGlobalDofIdx();
287  const int local_idx = dof.getPetscLocalDofIdx();
288  const int pos =
289  side_number * nb_dofs_on_vert + dof.getDofCoeffIdx();
290  nodes_indices[pos] = idx;
291  local_nodes_indices[pos] = local_idx;
292  if (brother_side_number != -1) {
293  const int elem_idx = brother_side_number * nb_dofs_on_vert +
294  (*dit)->getDofCoeffIdx();
295  nodes_indices[elem_idx] = idx;
296  local_nodes_indices[elem_idx] = local_idx;
297  }
298  }
299  }
300  }
301  }
302  }
303  } else {
304  nodes_indices.resize(0, false);
305  local_nodes_indices.resize(0, false);
306  }
307 
308  } else {
309  nodes_indices.resize(0, false);
310  local_nodes_indices.resize(0, false);
311  }
312 
314 }
315 
318  const std::string &field_name) const {
319 
320  struct Extractor {
321  boost::weak_ptr<EntityCacheNumeredDofs>
322  operator()(boost::shared_ptr<FieldEntity> &e) {
323  return e->entityCacheRowDofs;
324  }
325  };
326 
327  return getNodesIndices(field_name, getRowFieldEnts(),
328  data.dataOnEntities[MBVERTEX][0].getIndices(),
329  data.dataOnEntities[MBVERTEX][0].getLocalIndices(),
330  Extractor());
331 }
332 
335  const std::string &field_name) const {
336 
337  struct Extractor {
338  boost::weak_ptr<EntityCacheNumeredDofs>
339  operator()(boost::shared_ptr<FieldEntity> &e) {
340  return e->entityCacheColDofs;
341  }
342  };
343 
344  return getNodesIndices(field_name, getColFieldEnts(),
345  data.dataOnEntities[MBVERTEX][0].getIndices(),
346  data.dataOnEntities[MBVERTEX][0].getLocalIndices(),
347  Extractor());
348 }
349 
350 template <typename EXTRACTOR>
352  DataForcesAndSourcesCore &data, const std::string &field_name,
353  FieldEntity_vector_view &ents_field, const EntityType type_lo,
354  const EntityType type_hi, EXTRACTOR &&extractor) const {
356 
357  for (EntityType t = type_lo; t != type_hi; ++t) {
358  for (auto &dat : data.dataOnEntities[t]) {
359  dat.getIndices().resize(0, false);
360  dat.getLocalIndices().resize(0, false);
361  }
362  }
363 
364  auto bit_number = mField.get_field_bit_number(field_name);
365  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
366  bit_number, get_id_for_min_type(type_lo));
367  auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
368  cmp_uid_lo);
369  if (lo != ents_field.end()) {
370  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
371  bit_number, get_id_for_max_type(type_hi));
372  auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
373  if (lo != hi) {
374 
375  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
376 
377  for (auto it = lo; it != hi; ++it)
378  if (auto e = it->lock()) {
379 
380  const EntityType type = e->getEntType();
381  auto side_ptr = e->getSideNumberPtr();
382  const int side = side_ptr->side_number;
383  if (side >= 0) {
384  const int nb_dofs_on_ent = e->getNbDofsOnEnt();
385  const int brother_side = side_ptr->brother_side_number;
386  auto &dat = data.dataOnEntities[type][side];
387  auto &ent_field_indices = dat.getIndices();
388  auto &ent_field_local_indices = dat.getLocalIndices();
389 
390  ent_field_indices.resize(nb_dofs_on_ent, false);
391  ent_field_local_indices.resize(nb_dofs_on_ent, false);
392  std::fill(ent_field_indices.data().begin(),
393  ent_field_indices.data().end(), -1);
394  std::fill(ent_field_local_indices.data().begin(),
395  ent_field_local_indices.data().end(), -1);
396 
397  if (auto cache = extractor(e).lock()) {
398  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
399  const int idx = (*dit)->getEntDofIdx();
400  ent_field_indices[idx] = (*dit)->getPetscGlobalDofIdx();
401  ent_field_local_indices[idx] = (*dit)->getPetscLocalDofIdx();
402  }
403  }
404 
405  if (brother_side != -1) {
406  auto &dat_brother = data.dataOnEntities[type][brother_side];
407  dat_brother.getIndices().resize(nb_dofs_on_ent, false);
408  dat_brother.getLocalIndices().resize(nb_dofs_on_ent, false);
409  noalias(dat_brother.getIndices()) = dat.getIndices();
410  noalias(dat_brother.getLocalIndices()) = dat.getLocalIndices();
411  }
412  }
413  }
414  }
415  }
416 
418 }
419 
421  DataForcesAndSourcesCore &data, const std::string &field_name,
422  const EntityType type_lo, const EntityType type_hi) const {
423 
424  struct Extractor {
425  boost::weak_ptr<EntityCacheNumeredDofs>
426  operator()(boost::shared_ptr<FieldEntity> &e) {
427  return e->entityCacheRowDofs;
428  }
429  };
430 
431  return getEntityIndices(data, field_name, getRowFieldEnts(), type_lo, type_hi,
432  Extractor());
433 }
434 
436  DataForcesAndSourcesCore &data, const std::string &field_name,
437  const EntityType type_lo, const EntityType type_hi) const {
438 
439  struct Extractor {
440  boost::weak_ptr<EntityCacheNumeredDofs>
441  operator()(boost::shared_ptr<FieldEntity> &e) {
442  return e->entityCacheColDofs;
443  }
444  };
445 
446  return getEntityIndices(data, field_name, getColFieldEnts(), type_lo, type_hi,
447  Extractor());
448 }
449 
451  const std::string &field_name,
452  boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
453  VectorInt &indices) const {
455  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
456  auto dit = dofs->get<Unique_mi_tag>().lower_bound(
457  FieldEntity::getLoBitNumberUId((*field_it)->getBitNumber()));
458  auto hi_dit = dofs->get<Unique_mi_tag>().upper_bound(
459  FieldEntity::getHiBitNumberUId((*field_it)->getBitNumber()));
460  indices.resize(std::distance(dit, hi_dit));
461  for (; dit != hi_dit; dit++) {
462  int idx = (*dit)->getPetscGlobalDofIdx();
463  indices[(*dit)->getDofCoeffIdx()] = idx;
464  }
466 }
467 
469  DataForcesAndSourcesCore &data, const std::string &field_name) const {
471  if (data.dataOnEntities[MBENTITYSET].size() == 0) {
472  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
473  }
474  CHKERR getNoFieldIndices(field_name, getRowDofsPtr(),
475  data.dataOnEntities[MBENTITYSET][0].getIndices());
477 }
478 
480  DataForcesAndSourcesCore &data, const std::string &field_name) const {
482  if (data.dataOnEntities[MBENTITYSET].size() == 0) {
483  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
484  }
485  CHKERR getNoFieldIndices(field_name, getColDofsPtr(),
486  data.dataOnEntities[MBENTITYSET][0].getIndices());
488 }
489 
490 // ** Indices from problem **
491 
493  const std::string &field_name, const NumeredDofEntity_multiIndex &dofs,
494  VectorInt &nodes_indices) const {
496 
497  const Field *field_struture = mField.get_field_structure(field_name);
498  if (field_struture->getSpace() == H1) {
499 
500  int num_nodes;
501  CHKERR getNumberOfNodes(num_nodes);
502  nodes_indices.resize(field_struture->getNbOfCoeffs() * num_nodes, false);
503  std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
504 
505  auto &side_table = const_cast<SideNumber_multiIndex &>(
506  numeredEntFiniteElementPtr->getSideNumberTable());
507  auto siit = side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 0));
508  auto hi_siit =
509  side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 10000));
510 
511  int nn = 0;
512  for (; siit != hi_siit; siit++, nn++) {
513 
514  if (siit->get()->side_number == -1)
515  continue;
516 
517  auto bit_number = mField.get_field_bit_number(field_name);
518  const EntityHandle ent = siit->get()->ent;
519  auto dit = dofs.get<Unique_mi_tag>().lower_bound(
520  FieldEntity::getLoLocalEntityBitNumber(bit_number, ent));
521  auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
522  FieldEntity::getHiLocalEntityBitNumber(bit_number, ent));
523  for (; dit != hi_dit; dit++) {
524  nodes_indices[siit->get()->side_number * (*dit)->getNbOfCoeffs() +
525  (*dit)->getDofCoeffIdx()] =
526  (*dit)->getPetscGlobalDofIdx();
527  }
528  }
529  } else {
530  nodes_indices.resize(0, false);
531  }
532 
534 }
535 
537  const std::string &field_name, const NumeredDofEntity_multiIndex &dofs,
538  EntityType type, int side_number, VectorInt &indices) const {
540 
541  indices.resize(0);
542 
543  auto &side_table = const_cast<SideNumber_multiIndex &>(
544  numeredEntFiniteElementPtr->getSideNumberTable());
545  auto siit =
546  side_table.get<1>().lower_bound(boost::make_tuple(type, side_number));
547  auto hi_siit =
548  side_table.get<1>().upper_bound(boost::make_tuple(type, side_number));
549 
550  for (; siit != hi_siit; siit++) {
551 
552  if (siit->get()->side_number == -1)
553  continue;
554 
555  const EntityHandle ent = siit->get()->ent;
556  auto bit_number = mField.get_field_bit_number(field_name);
557  auto dit = dofs.get<Unique_mi_tag>().lower_bound(
558  FieldEntity::getLoLocalEntityBitNumber(bit_number, ent));
559  auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
560  FieldEntity::getHiLocalEntityBitNumber(bit_number, ent));
561  indices.resize(std::distance(dit, hi_dit));
562  for (; dit != hi_dit; dit++) {
563 
564  indices[(*dit)->getEntDofIdx()] = (*dit)->getPetscGlobalDofIdx();
565  }
566  }
567 
569 }
570 
572  const std::string &field_name, VectorInt &nodes_indices) const {
573  return getProblemNodesIndices(field_name, *(problemPtr->numeredRowDofsPtr),
574  nodes_indices);
575 }
576 
579  EntityType type, int side_number,
580  VectorInt &indices) const {
581  return getProblemTypeIndices(field_name, *(problemPtr->numeredRowDofsPtr),
582  type, side_number, indices);
583 }
584 
586  const std::string &field_name, VectorInt &nodes_indices) const {
587  return getProblemNodesIndices(field_name, *(problemPtr->numeredColDofsPtr),
588  nodes_indices);
589 }
590 
593  EntityType type, int side_number,
594  VectorInt &indices) const {
595  return getProblemTypeIndices(field_name, *(problemPtr->numeredColDofsPtr),
596  type, side_number, indices);
597 }
598 
599 // ** Data **
600 
603  const std::string &field_name) const {
604 
605  auto get_nodes_field_data = [&](VectorDouble &nodes_data,
606  VectorFieldEntities &field_entities,
607  VectorDofs &nodes_dofs, FieldSpace &space,
609  VectorInt &bb_node_order) {
611 
612  nodes_data.resize(0, false);
613  nodes_dofs.resize(0, false);
614  field_entities.resize(0, false);
615 
616  auto field_it = fieldsPtr->get<FieldName_mi_tag>().find(field_name);
617  if (field_it != fieldsPtr->get<FieldName_mi_tag>().end()) {
618 
619  auto bit_number = (*field_it)->getBitNumber();
620  const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
621  space = (*field_it)->getSpace();
622  base = (*field_it)->getApproxBase();
623 
624  auto &field_ents = getDataFieldEnts();
625  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
626  bit_number, get_id_for_min_type<MBVERTEX>());
627  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
628  cmp_uid_lo);
629  if (lo != field_ents.end()) {
630  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
631  bit_number, get_id_for_max_type<MBVERTEX>());
632  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
633  if (lo != hi) {
634 
635  int nb_dofs = 0;
636  for (auto it = lo; it != hi; ++it) {
637  if (auto e = it->lock()) {
638  nb_dofs += e->getNbDofsOnEnt();
639  }
640  }
641 
642  if (nb_dofs) {
643 
644  int num_nodes;
645  CHKERR getNumberOfNodes(num_nodes);
646  bb_node_order.resize(num_nodes, false);
647  bb_node_order.clear();
648  const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
649  nodes_data.resize(max_nb_dofs, false);
650  nodes_dofs.resize(max_nb_dofs, false);
651  field_entities.resize(num_nodes, false);
652  std::fill(nodes_data.begin(), nodes_data.end(), 0);
653  std::fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
654  std::fill(field_entities.begin(), field_entities.end(), nullptr);
655 
656  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
657 
658  for (auto it = lo; it != hi; ++it) {
659  if (auto e = it->lock()) {
660  const auto &sn = e->getSideNumberPtr();
661  const int side_number = sn->side_number;
662  // Some field entities on skeleton can have negative side
663  // numbeer
664  if (side_number >= 0) {
665  const int brother_side_number = sn->brother_side_number;
666 
667  field_entities[side_number] = e.get();
668  if (brother_side_number != -1) {
669  brother_ents_vec.emplace_back(e);
670  field_entities[side_number] = field_entities[side_number];
671  }
672 
673  bb_node_order[side_number] = e->getMaxOrder();
674  int pos = side_number * nb_dofs_on_vert;
675  auto ent_filed_data_vec = e->getEntFieldData();
676  if (auto cache = e->entityCacheDataDofs.lock()) {
677  for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
678  ++dit) {
679  const auto dof_idx = (*dit)->getEntDofIdx();
680  nodes_data[pos + dof_idx] = ent_filed_data_vec[dof_idx];
681  nodes_dofs[pos + dof_idx] =
682  reinterpret_cast<FEDofEntity *>((*dit).get());
683  }
684  }
685  }
686  }
687  }
688 
689  for (auto &it : brother_ents_vec) {
690  if (const auto e = it.lock()) {
691  const auto &sn = e->getSideNumberPtr();
692  const int side_number = sn->side_number;
693  const int brother_side_number = sn->brother_side_number;
694  bb_node_order[brother_side_number] = bb_node_order[side_number];
695  int pos = side_number * nb_dofs_on_vert;
696  int brother_pos = brother_side_number * nb_dofs_on_vert;
697  for (int ii = 0; ii != nb_dofs_on_vert; ++ii) {
698  nodes_data[brother_pos] = nodes_data[pos];
699  nodes_dofs[brother_pos] = nodes_dofs[pos];
700  ++pos;
701  ++brother_pos;
702  }
703  }
704  }
705  }
706  }
707  }
708  }
709 
711  };
712 
713  return get_nodes_field_data(
714  data.dataOnEntities[MBVERTEX][0].getFieldData(),
715  data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
716  data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
717  data.dataOnEntities[MBVERTEX][0].getSpace(),
718  data.dataOnEntities[MBVERTEX][0].getBase(),
719  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder());
720 }
721 
723  DataForcesAndSourcesCore &data, const std::string &field_name,
724  const EntityType type_lo, const EntityType type_hi) const {
726  for (EntityType t = type_lo; t != type_hi; ++t) {
727  for (auto &dat : data.dataOnEntities[t]) {
728  dat.getDataOrder() = 0;
729  dat.getBase() = NOBASE;
730  dat.getSpace() = NOSPACE;
731  dat.getFieldData().resize(0, false);
732  dat.getFieldDofs().resize(0, false);
733  dat.getFieldEntities().resize(0, false);
734  }
735  }
736 
737  auto &field_ents = getDataFieldEnts();
738  auto bit_number = mField.get_field_bit_number(field_name);
739  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
740  bit_number, get_id_for_min_type(type_lo));
741  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
742  cmp_uid_lo);
743  if (lo != field_ents.end()) {
744  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
745  bit_number, get_id_for_max_type(type_hi));
746  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
747  if (lo != hi) {
748 
749  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
750 
751  for (auto it = lo; it != hi; ++it)
752  if (auto e = it->lock()) {
753 
754  const EntityType type = e->getEntType();
755  auto side_ptr = e->getSideNumberPtr();
756  const int side = side_ptr->side_number;
757  if (side >= 0) {
758 
759  auto &dat = data.dataOnEntities[type][side];
760  auto &ent_field = dat.getFieldEntities();
761  auto &ent_field_dofs = dat.getFieldDofs();
762  auto &ent_field_data = dat.getFieldData();
763 
764  const int brother_side = side_ptr->brother_side_number;
765  if (brother_side != -1)
766  brother_ents_vec.emplace_back(e);
767 
768  dat.getBase() = e->getApproxBase();
769  dat.getSpace() = e->getSpace();
770  const int ent_order = e->getMaxOrder();
771  dat.getDataOrder() =
772  dat.getDataOrder() > ent_order ? dat.getDataOrder() : ent_order;
773 
774  auto ent_data = e->getEntFieldData();
775  ent_field_data.resize(ent_data.size(), false);
776  noalias(ent_field_data) = ent_data;
777  ent_field_dofs.resize(ent_data.size(), false);
778  std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
779  ent_field.resize(1, false);
780  ent_field[0] = e.get();
781  if (auto cache = e->entityCacheDataDofs.lock()) {
782  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
783  ent_field_dofs[(*dit)->getEntDofIdx()] =
784  reinterpret_cast<FEDofEntity *>((*dit).get());
785  }
786  }
787  }
788  }
789 
790  for (auto &it : brother_ents_vec) {
791  if (const auto e = it.lock()) {
792  const EntityType type = e->getEntType();
793  const int side = e->getSideNumberPtr()->side_number;
794  const int brother_side = e->getSideNumberPtr()->brother_side_number;
795  auto &dat = data.dataOnEntities[type][side];
796  auto &dat_brother = data.dataOnEntities[type][brother_side];
797  dat_brother.getBase() = dat.getBase();
798  dat_brother.getSpace() = dat.getSpace();
799  dat_brother.getDataOrder() = dat.getDataOrder();
800  dat_brother.getFieldData() = dat.getFieldData();
801  dat_brother.getFieldDofs() = dat.getFieldDofs();
802  dat_brother.getFieldEntities() = dat.getFieldEntities();
803  }
804  }
805  }
806  }
807 
809 }
810 
812  const std::string field_name, VectorDouble &ent_field_data,
813  VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const {
814 
816 
817  ent_field_data.resize(0, false);
818  ent_field_dofs.resize(0, false);
819  ent_field.resize(0, false);
820 
821  auto &field_ents = getDataFieldEnts();
822  auto bit_number = mField.get_field_bit_number(field_name);
823  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
824  bit_number, get_id_for_min_type<MBVERTEX>());
825  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
826  cmp_uid_lo);
827  if (lo != field_ents.end()) {
828 
829  ent_field.resize(field_ents.size(), false);
830  std::fill(ent_field.begin(), ent_field.end(), nullptr);
831 
832  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
833  bit_number, get_id_for_max_type<MBVERTEX>());
834  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
835  if (lo != hi) {
836 
837  int side = 0;
838  for (auto it = lo; it != hi; ++it, ++side)
839  if (auto e = it->lock()) {
840 
841  const auto size = e->getNbDofsOnEnt();
842  ent_field_data.resize(size, false);
843  ent_field_dofs.resize(size, false);
844  ent_field[side] = e.get();
845  noalias(ent_field_data) = e->getEntFieldData();
846 
847  if (auto cache = e->entityCacheDataDofs.lock()) {
848  for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
849  ent_field_dofs[(*dit)->getEntDofIdx()] =
850  reinterpret_cast<FEDofEntity *>((*dit).get());
851  }
852  }
853  }
854  }
855  }
856 
858 }
859 
862  const std::string field_name) const {
864  if (data.dataOnEntities[MBENTITYSET].size() == 0)
865  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
866  "No space to insert data");
867 
869  field_name, data.dataOnEntities[MBENTITYSET][0].getFieldData(),
870  data.dataOnEntities[MBENTITYSET][0].getFieldDofs(),
871  data.dataOnEntities[MBENTITYSET][0].getFieldEntities());
873 }
874 
875 // ** Face **
876 
880  // PetscAttachDebugger();
881  data.facesNodes.resize(4, 3, false);
882  auto &side_table = const_cast<SideNumber_multiIndex &>(
883  numeredEntFiniteElementPtr->getSideNumberTable());
884  auto siit = side_table.get<1>().lower_bound(boost::make_tuple(MBTRI, 0));
885  auto hi_siit = side_table.get<1>().upper_bound(boost::make_tuple(MBTRI, 4));
886  if (std::distance(siit, hi_siit) != 4) {
887  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
888  "Should be 4 triangles on tet, side_table not initialized");
889  }
890  const int canonical_face_sense_p1[4][3] = {
891  {0, 1, 3},
892  {1, 2, 3},
893  {0, 3, 2} /**/,
894  {0, 2, 1} /**/}; // second index is offset (positive sense)
895  const int canonical_face_sense_m1[4][3] = {
896  {0, 3, 1},
897  {1, 3, 2},
898  {0, 2, 3},
899  {0, 1, 2}}; // second index is offset (negative sense
900  for (; siit != hi_siit; siit++) {
901  const boost::shared_ptr<SideNumber> side = *siit;
902  int face_conn[3] = {-1, -1, -1};
903  if (side->offset == 0) {
904  face_conn[0] = side->sense == 1
905  ? canonical_face_sense_p1[(int)side->side_number][0]
906  : canonical_face_sense_m1[(int)side->side_number][0];
907  face_conn[1] = side->sense == 1
908  ? canonical_face_sense_p1[(int)side->side_number][1]
909  : canonical_face_sense_m1[(int)side->side_number][1];
910  face_conn[2] = side->sense == 1
911  ? canonical_face_sense_p1[(int)side->side_number][2]
912  : canonical_face_sense_m1[(int)side->side_number][2];
913  }
914  if (side->offset == 1) {
915  face_conn[0] =
916  side->sense == 1
917  ? canonical_face_sense_p1[(int)side->side_number][1]
918  : canonical_face_sense_m1[(int)side->side_number][2] /**/;
919  face_conn[1] = side->sense == 1
920  ? canonical_face_sense_p1[(int)side->side_number][2]
921  : canonical_face_sense_m1[(int)side->side_number][0];
922  face_conn[2] = side->sense == 1
923  ? canonical_face_sense_p1[(int)side->side_number][0]
924  : canonical_face_sense_m1[(int)side->side_number][1];
925  }
926  if (side->offset == 2) {
927  face_conn[0] =
928  side->sense == 1
929  ? canonical_face_sense_p1[(int)side->side_number][2]
930  : canonical_face_sense_m1[(int)side->side_number][1] /**/;
931  face_conn[1] = side->sense == 1
932  ? canonical_face_sense_p1[(int)side->side_number][0]
933  : canonical_face_sense_m1[(int)side->side_number][2];
934  face_conn[2] = side->sense == 1
935  ? canonical_face_sense_p1[(int)side->side_number][1]
936  : canonical_face_sense_m1[(int)side->side_number][0];
937  }
938  for (int nn = 0; nn < 3; nn++)
939  data.facesNodes(side->side_number, nn) = face_conn[nn];
940  {
941  const EntityHandle *conn_tet;
942  int num_nodes_tet;
944  CHKERR mField.get_moab().get_connectivity(ent, conn_tet, num_nodes_tet,
945  true);
946  if (num_nodes_tet != 4)
947  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
948  "data inconsistency");
949  int num_nodes_face;
950  const EntityHandle *conn_face;
951  CHKERR mField.get_moab().get_connectivity(side->ent, conn_face,
952  num_nodes_face, true);
953  if (num_nodes_face != 3)
954  SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
955  if (conn_face[0] != conn_tet[data.facesNodes(side->side_number, 0)])
956  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
957  "data inconsistency");
958  if (conn_face[1] != conn_tet[data.facesNodes(side->side_number, 1)])
959  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
960  "data inconsistency");
961  if (conn_face[2] != conn_tet[data.facesNodes(side->side_number, 2)])
962  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
963  "data inconsistency");
964  }
965  }
967 }
968 
969 // ** Space and Base **
970 
972  DataForcesAndSourcesCore &data) const {
974 
975  if (nInTheLoop == 0) {
976  data.sPace.reset();
977  data.bAse.reset();
978  for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
979  data.spacesOnEntities[t].reset();
980  data.basesOnEntities[t].reset();
981  }
982  for (int s = 0; s != LASTSPACE; ++s) {
983  data.basesOnSpaces[s].reset();
984  }
985  }
986 
987  auto fe_type = numeredEntFiniteElementPtr->getEntType();
988 
989  if (getDataFieldEntsPtr())
990  for (auto e : getDataFieldEnts()) {
991  if (auto ptr = e.lock()) {
992  // get data from entity
993  const EntityType type = ptr->getEntType();
995  const FieldSpace space = ptr->getSpace();
996  const FieldApproximationBase approx = ptr->getApproxBase();
997  // set data
998  data.sPace.set(space);
999  data.bAse.set(approx);
1000  data.spacesOnEntities[type].set(space);
1001  data.basesOnEntities[type].set(approx);
1002  data.basesOnSpaces[space].set(approx);
1003  }
1004  }
1005  }
1006  else
1007  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1008  "data fields ents not allocated on element");
1009 
1011 }
1012 
1014  const FieldApproximationBase b) {
1016  if (dataOnElement[H1]->bAse.test(b)) {
1017  switch (static_cast<FieldApproximationBase>(b)) {
1018  case NOBASE:
1019  break;
1021  // BERNSTEIN_BEZIER_BASE is not hierarchical base
1022  break;
1025  case DEMKOWICZ_JACOBI_BASE:
1026  if (!getElementPolynomialBase())
1027  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1028  "Functions genrating approximation base not defined");
1029 
1030  for (int space = H1; space != LASTSPACE; ++space) {
1031  if (dataOnElement[H1]->sPace.test(space) &&
1032  dataOnElement[H1]->bAse.test(b) &&
1033  dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1034  CHKERR getElementPolynomialBase()->getValue(
1035  gaussPts,
1036  boost::make_shared<EntPolynomialBaseCtx>(
1037  *dataOnElement[space], static_cast<FieldSpace>(space),
1038  static_cast<FieldApproximationBase>(b), NOBASE));
1039  }
1040  }
1041  break;
1042  case USER_BASE:
1043  if (!getUserPolynomialBase())
1044  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1045  "Functions genrating user approximation base not defined");
1046 
1047  for (int space = H1; space != LASTSPACE; ++space)
1048  if (dataOnElement[H1]->sPace.test(space) &&
1049  dataOnElement[H1]->bAse.test(b) &&
1050  dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1051  CHKERR getUserPolynomialBase()->getValue(
1052  gaussPts,
1053  boost::make_shared<EntPolynomialBaseCtx>(
1054  *dataOnElement[space], static_cast<FieldSpace>(space),
1055  static_cast<FieldApproximationBase>(b), NOBASE));
1056  }
1057  break;
1058  default:
1060  "Base <%s> not yet implemented",
1062  }
1063  }
1065 }
1066 
1069  /// Use the some node base. Node base is usually used for construction other
1070  /// bases.
1071  for (int space = HCURL; space != LASTSPACE; ++space) {
1072  dataOnElement[space]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
1073  dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
1074  dataOnElement[space]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1075  NOBASE) =
1076  dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1077  NOBASE);
1078  }
1079  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
1081  static_cast<FieldApproximationBase>(b));
1082  }
1084 }
1085 
1089 
1090  auto get_nodal_base_data = [&](DataForcesAndSourcesCore &data,
1091  auto field_ptr) {
1093  auto &space = data.dataOnEntities[MBVERTEX][0].getSpace();
1094  auto &base = data.dataOnEntities[MBVERTEX][0].getBase();
1095  auto &bb_node_order = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder();
1096 
1097  auto &field_ents = getDataFieldEnts();
1098  auto bit_number = field_ptr->getBitNumber();
1099  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1100  bit_number, get_id_for_min_type<MBVERTEX>());
1101  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1102  cmp_uid_lo);
1103  if (lo != field_ents.end()) {
1104  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1105  bit_number, get_id_for_max_type<MBVERTEX>());
1106  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1107  if (lo != hi) {
1108 
1109  for (auto it = lo; it != hi; ++it)
1110  if (auto first_e = it->lock()) {
1111  space = first_e->getSpace();
1112  base = first_e->getApproxBase();
1113  int num_nodes;
1114  CHKERR getNumberOfNodes(num_nodes);
1115  bb_node_order.resize(num_nodes, false);
1116  bb_node_order.clear();
1117  const int nb_dof_idx = first_e->getNbOfCoeffs();
1118 
1119  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1120 
1121  for (; it != hi; ++it) {
1122  if (auto e = it->lock()) {
1123  const auto &sn = e->getSideNumberPtr();
1124  const int side_number = sn->side_number;
1125  const int brother_side_number = sn->brother_side_number;
1126  if (brother_side_number != -1)
1127  brother_ents_vec.emplace_back(e);
1128  bb_node_order[side_number] = e->getMaxOrder();
1129  }
1130  }
1131 
1132  for (auto &it : brother_ents_vec) {
1133  if (const auto e = it.lock()) {
1134  const auto &sn = e->getSideNumberPtr();
1135  const int side_number = sn->side_number;
1136  const int brother_side_number = sn->brother_side_number;
1137  bb_node_order[brother_side_number] = bb_node_order[side_number];
1138  }
1139  }
1140 
1141  break;
1142  }
1143  }
1144  }
1145 
1147  };
1148 
1149  auto get_entity_base_data = [&](DataForcesAndSourcesCore &data,
1150  auto field_ptr, const EntityType type_lo,
1151  const EntityType type_hi) {
1153  for (EntityType t = type_lo; t != type_hi; ++t) {
1154  for (auto &dat : data.dataOnEntities[t]) {
1155  dat.getDataOrder() = 0;
1156  dat.getBase() = NOBASE;
1157  dat.getSpace() = NOSPACE;
1158  dat.getFieldData().resize(0, false);
1159  dat.getFieldDofs().resize(0, false);
1160  }
1161  }
1162 
1163  auto &field_ents = getDataFieldEnts();
1164  auto bit_number = field_ptr->getBitNumber();
1165  const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1166  bit_number, get_id_for_min_type(type_lo));
1167  auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1168  cmp_uid_lo);
1169  if (lo != field_ents.end()) {
1170  const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1171  bit_number, get_id_for_max_type(type_hi));
1172  auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1173 
1174  std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1175  for (; lo != hi; ++lo) {
1176  if (auto e = lo->lock()) {
1177  if (auto cache = e->entityCacheDataDofs.lock()) {
1178  if (cache->loHi[0] != cache->loHi[1]) {
1179  const EntityType type = e->getEntType();
1180  const int side = e->getSideNumberPtr()->side_number;
1181  auto &dat = data.dataOnEntities[type][side];
1182  const int brother_side =
1183  e->getSideNumberPtr()->brother_side_number;
1184  if (brother_side != -1)
1185  brother_ents_vec.emplace_back(e);
1186  dat.getBase() = e->getApproxBase();
1187  dat.getSpace() = e->getSpace();
1188  const auto ent_order = e->getMaxOrder();
1189  dat.getDataOrder() = dat.getDataOrder() > ent_order
1190  ? dat.getDataOrder()
1191  : ent_order;
1192  }
1193  }
1194  }
1195  }
1196 
1197  for (auto &ent_ptr : brother_ents_vec) {
1198  if (auto e = ent_ptr.lock()) {
1199  const EntityType type = e->getEntType();
1200  const int side = e->getSideNumberPtr()->side_number;
1201  const int brother_side = e->getSideNumberPtr()->brother_side_number;
1202  auto &dat = data.dataOnEntities[type][side];
1203  auto &dat_brother = data.dataOnEntities[type][brother_side];
1204  dat_brother.getBase() = dat.getBase();
1205  dat_brother.getSpace() = dat.getSpace();
1206  dat_brother.getDataOrder() = dat.getDataOrder();
1207  }
1208  }
1209  }
1211  };
1212 
1213  for (auto &e : getDataFieldEnts()) {
1214  if (auto ent_data_ptr = e.lock()) {
1215  if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1216  auto space = ent_data_ptr->getSpace();
1217  for (EntityType t = MBVERTEX; t != MBPOLYHEDRON; ++t) {
1218  for (auto &dat : (*dataOnElement[space]).dataOnEntities[t]) {
1219  for (auto &ptr : dat.getBBAlphaIndicesByOrderArray())
1220  ptr.reset();
1221  for (auto &ptr : dat.getBBNByOrderArray())
1222  ptr.reset();
1223  for (auto &ptr : dat.getBBDiffNByOrderArray())
1224  ptr.reset();
1225  }
1226  }
1227  }
1228  }
1229  }
1230 
1231  std::set<string> fields_list;
1232  for (auto &e : getDataFieldEnts()) {
1233  if (auto ent_data_ptr = e.lock()) {
1234  if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1235  auto field_name = ent_data_ptr->getName();
1236  if (fields_list.find(field_name) == fields_list.end()) {
1237  auto field_ptr = ent_data_ptr->getFieldRawPtr();
1238  auto space = ent_data_ptr->getSpace();
1239  CHKERR get_nodal_base_data(*dataOnElement[space], field_ptr);
1240  CHKERR get_entity_base_data(*dataOnElement[space], field_ptr, MBEDGE,
1241  MBPOLYHEDRON);
1242  CHKERR getElementPolynomialBase()->getValue(
1243  gaussPts, boost::make_shared<EntPolynomialBaseCtx>(
1244  *dataOnElement[space], field_name,
1245  static_cast<FieldSpace>(space),
1247  fields_list.insert(field_name);
1248  }
1249  }
1250  }
1251  }
1252 
1254 };
1255 
1258 
1259  const EntityType type = numeredEntFiniteElementPtr->getEntType();
1262 
1263  // Data on elements for proper spaces
1264  for (int space = H1; space != LASTSPACE; ++space) {
1265  dataOnElement[space]->setElementType(type);
1266  derivedDataOnElement[space]->setElementType(type);
1267  }
1268 
1270 
1272 }
1273 
1274 #define FUNCTION_NAME_WITH_OP_NAME(OP) \
1275  std::ostringstream ss; \
1276  ss << "(Calling user data operator " \
1277  << boost::typeindex::type_id_runtime(OP).pretty_name() << " rowField " \
1278  << (OP).rowFieldName << " colField " << (OP).colFieldName << ") "
1279 
1280 #define CATCH_OP_ERRORS(OP) \
1281  catch (MoFEMExceptionInitial const &ex) { \
1282  FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1283  return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1284  ex.errorCode, PETSC_ERROR_INITIAL, ex.what()); \
1285  } \
1286  catch (MoFEMExceptionRepeat const &ex) { \
1287  FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1288  return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1289  ex.errorCode, PETSC_ERROR_REPEAT, " "); \
1290  } \
1291  catch (MoFEMException const &ex) { \
1292  FUNCTION_NAME_WITH_OP_NAME(OP) << ex.errorMessage; \
1293  SETERRQ(PETSC_COMM_SELF, ex.errorCode, ss.str().c_str()); \
1294  } \
1295  catch (std::exception const &ex) { \
1296  std::string message("Error: " + std::string(ex.what()) + " at " + \
1297  boost::lexical_cast<std::string>(__LINE__) + " : " + \
1298  std::string(__FILE__) + " in " + \
1299  std::string(PETSC_FUNCTION_NAME)); \
1300  FUNCTION_NAME_WITH_OP_NAME(OP) << message; \
1301  SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str()); \
1302  }
1303 
1306 
1307  const EntityType type = numeredEntFiniteElementPtr->getEntType();
1310  std::vector<std::string> last_eval_field_name(2);
1311 
1312  boost::ptr_vector<UserDataOperator>::iterator oit, hi_oit;
1313  oit = opPtrVector.begin();
1314  hi_oit = opPtrVector.end();
1315 
1316  for (; oit != hi_oit; oit++) {
1317 
1318  try {
1319 
1320  CHKERR oit->setPtrFE(this);
1321 
1322  if (oit->opType == UserDataOperator::OPLAST) {
1323 
1324  // Set field
1325  switch (oit->sPace) {
1326  case NOSPACE:
1327  CHKERR oit->doWork(
1328  0, MBENTITYSET,
1329  dataOnElement[oit->sPace]->dataOnEntities[MBENTITYSET][0]);
1330  break;
1331  case NOFIELD:
1332  case H1:
1333  case HCURL:
1334  case HDIV:
1335  case L2:
1336  break;
1337  default:
1338  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1339  "Not implemented for this space", oit->sPace);
1340  }
1341 
1342  // Reseat all data which all field dependent
1343  dataOnElement[oit->sPace]->resetFieldDependentData();
1344  last_eval_field_name[0] = "";
1345 
1346  // Run operator
1347  try {
1348  CHKERR oit->opRhs(*dataOnElement[oit->sPace], false);
1349  }
1350  CATCH_OP_ERRORS(*oit);
1351 
1352  } else {
1353 
1354  boost::shared_ptr<DataForcesAndSourcesCore> op_data[2];
1355  std::array<bool, 2> base_swap;
1356  std::array<std::pair<std::string, FieldApproximationBase>, 2>
1357  base_swap_data;
1358  auto swap_bases = [&]() {
1360  for (size_t ss = 0; ss != 2; ++ss)
1361  if (base_swap[ss])
1362  CHKERR op_data[ss]->baseSwap(base_swap_data[ss].first,
1363  base_swap_data[ss].second);
1365  };
1366 
1367  for (size_t ss = 0; ss != 2; ss++) {
1368 
1369  const std::string field_name =
1370  !ss ? oit->rowFieldName : oit->colFieldName;
1371  if (field_name.empty()) {
1372  SETERRQ2(
1373  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1374  "No field name in operator %d (0-row, 1-column) in operator %s",
1375  ss,
1376  (boost::typeindex::type_id_runtime(*oit).pretty_name())
1377  .c_str());
1378  }
1379  const Field *field_struture = mField.get_field_structure(field_name);
1380  const BitFieldId data_id = field_struture->getId();
1381  const FieldSpace space = field_struture->getSpace();
1382  const FieldApproximationBase base = field_struture->getApproxBase();
1383  op_data[ss] =
1384  !ss ? dataOnElement[space] : derivedDataOnElement[space];
1385 
1386  switch (base) {
1388  base_swap_data[ss] = std::pair<std::string, FieldApproximationBase>(
1389  field_name, AINSWORTH_BERNSTEIN_BEZIER_BASE);
1390  base_swap[ss] = true;
1391  break;
1392  default:
1393  base_swap[ss] = false;
1394  };
1395 
1396  if ((oit->getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
1397  data_id)
1398  .none()) {
1399  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1400  "no data field < %s > on finite element < %s >",
1401  field_name.c_str(), feName.c_str());
1402  }
1403 
1404  if (oit->getOpType() & types[ss] ||
1405  oit->getOpType() & UserDataOperator::OPROWCOL) {
1406 
1407  switch (space) {
1408  case NOSPACE:
1409  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1410  "unknown space");
1411  break;
1412  case NOFIELD:
1413  case H1:
1414  case HCURL:
1415  case HDIV:
1416  case L2:
1417  break;
1418  default:
1419  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1420  "Not implemented for this space", space);
1421  }
1422 
1423  if (last_eval_field_name[ss] != field_name) {
1424 
1425  CHKERR getEntityFieldData(*op_data[ss], field_name, MBEDGE);
1426  if (!ss)
1427  CHKERR getEntityRowIndices(*op_data[ss], field_name, MBEDGE);
1428  else
1429  CHKERR getEntityColIndices(*op_data[ss], field_name, MBEDGE);
1430 
1431  switch (space) {
1432  case NOSPACE:
1433  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1434  "unknown space");
1435  break;
1436  case H1:
1437  if (!ss)
1438  CHKERR getRowNodesIndices(*op_data[ss], field_name);
1439  else
1440  CHKERR getColNodesIndices(*op_data[ss], field_name);
1441  CHKERR getNodesFieldData(*op_data[ss], field_name);
1442  break;
1443  case HCURL:
1444  case HDIV:
1445  break;
1446  case L2:
1447  switch (type) {
1448  case MBVERTEX:
1449  CHKERR getNodesFieldData(*op_data[ss], field_name);
1450  break;
1451  default:
1452  break;
1453  }
1454  break;
1455  case NOFIELD:
1456  if (!getNinTheLoop()) {
1457  // NOFIELD data are the same for each element, can be
1458  // retrieved only once
1459  if (!ss) {
1460  CHKERR getNoFieldRowIndices(*op_data[ss], field_name);
1461  } else {
1462  CHKERR getNoFieldColIndices(*op_data[ss], field_name);
1463  }
1464  CHKERR getNoFieldFieldData(*op_data[ss], field_name);
1465  }
1466  break;
1467  default:
1468  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1469  "Not implemented for this space", space);
1470  }
1471  last_eval_field_name[ss] = field_name;
1472  }
1473  }
1474  }
1475 
1476  CHKERR swap_bases();
1477 
1478  if (oit->getOpType() & UserDataOperator::OPROW) {
1479  try {
1480  CHKERR oit->opRhs(*op_data[0], false);
1481  }
1482  CATCH_OP_ERRORS(*oit);
1483  }
1484 
1485  if (oit->getOpType() & UserDataOperator::OPCOL) {
1486  try {
1487  CHKERR oit->opRhs(*op_data[1], false);
1488  }
1489  CATCH_OP_ERRORS(*oit);
1490  }
1491 
1492  if (oit->getOpType() & UserDataOperator::OPROWCOL) {
1493  try {
1494  CHKERR oit->opLhs(*op_data[0], *op_data[1]);
1495  }
1496  CATCH_OP_ERRORS(*oit);
1497  }
1498 
1499  CHKERR swap_bases();
1500  }
1501  }
1502  CATCH_OP_ERRORS(*oit);
1503  }
1505 }
1506 
1508  const std::string field_name, const EntityType type, const int side,
1509  VectorInt &indices) const {
1511  if (ptrFE == NULL)
1512  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1513 
1514  switch (type) {
1515  case MBVERTEX:
1516  CHKERR ptrFE->getProblemNodesRowIndices(field_name, indices);
1517  break;
1518  default:
1519  CHKERR ptrFE->getProblemTypeRowIndices(field_name, type, side, indices);
1520  }
1522 }
1523 
1525  const std::string field_name, const EntityType type, const int side,
1526  VectorInt &indices) const {
1528  if (ptrFE == NULL)
1529  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1530 
1531  switch (type) {
1532  case MBVERTEX:
1533  CHKERR ptrFE->getProblemNodesColIndices(field_name, indices);
1534  break;
1535  default:
1536  CHKERR ptrFE->getProblemTypeColIndices(field_name, type, side, indices);
1537  }
1539 }
1540 
1544  sidePtrFE = const_cast<ForcesAndSourcesCore *>(side_fe_ptr);
1546 }
1547 
1549  const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t side_dim,
1550  const EntityHandle ent_for_side) {
1552  const EntityHandle ent = ent_for_side ? ent_for_side : getFEEntityHandle();
1553 
1554  const Problem *problem_ptr = getFEMethod()->problemPtr;
1555  Range adjacent_ents;
1556  CHKERR ptrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
1557  ent, side_dim, adjacent_ents);
1558  typedef NumeredEntFiniteElement_multiIndex::index<
1559  Composite_Name_And_Ent_mi_tag>::type FEByComposite;
1560  FEByComposite &numered_fe = problem_ptr->numeredFiniteElementsPtr
1562 
1563  side_fe->feName = fe_name;
1564 
1565  CHKERR side_fe->setSideFEPtr(ptrFE);
1566  CHKERR side_fe->copyBasicMethod(*getFEMethod());
1567  CHKERR side_fe->copyKsp(*getFEMethod());
1568  CHKERR side_fe->copySnes(*getFEMethod());
1569  CHKERR side_fe->copyTs(*getFEMethod());
1570 
1571  CHKERR side_fe->preProcess();
1572 
1573  int nn = 0;
1574  side_fe->loopSize = adjacent_ents.size();
1575  for (Range::iterator vit = adjacent_ents.begin(); vit != adjacent_ents.end();
1576  vit++) {
1577  FEByComposite::iterator miit =
1578  numered_fe.find(boost::make_tuple(fe_name, *vit));
1579  if (miit != numered_fe.end()) {
1580  side_fe->nInTheLoop = nn++;
1581  side_fe->numeredEntFiniteElementPtr = *miit;
1582  CHKERR (*side_fe)();
1583  }
1584  }
1585 
1586  CHKERR side_fe->postProcess();
1587 
1589 }
1590 
1591 int ForcesAndSourcesCore::getRule(int order_row, int order_col,
1592  int order_data) {
1593  return getRuleHook ? getRuleHook(order_row, order_col, order_data)
1594  : getRule(order_data);
1595 }
1596 
1598  int order_data) {
1599  return setRuleHook ? setRuleHook(order_row, order_col, order_data)
1600  : setGaussPts(order_data);
1601 }
1602 
1604 
1605 /** \deprecated setGaussPts(int row_order, int col_order, int data order);
1606  */
1609  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Sorry, not implemented");
1611 }
1612 
1614  const char type,
1615  const bool symm)
1616  : DataOperator(symm), opType(type), sPace(space), ptrFE(nullptr) {}
1617 
1619  const std::string &field_name, const char type, const bool symm)
1620  : DataOperator(symm), opType(type), rowFieldName(field_name),
1621  colFieldName(field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1622 
1624  const std::string &row_field_name, const std::string &col_field_name,
1625  const char type, const bool symm)
1626  : DataOperator(symm), opType(type), rowFieldName(row_field_name),
1627  colFieldName(col_field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1628 
1631  if (preProcessHook) {
1632  ierr = preProcessHook();
1633  CHKERRG(ierr);
1634  }
1636 }
1639  if (operatorHook) {
1640  ierr = operatorHook();
1641  CHKERRG(ierr);
1642  }
1644 }
1647  if (postProcessHook) {
1648  ierr = postProcessHook();
1649  CHKERRG(ierr);
1650  }
1652 }
1653 
1657  if (!(ptrFE = dynamic_cast<ForcesAndSourcesCore *>(ptr)))
1658  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1659  "User operator and finite element do not work together");
1661 }
1662 
1663 } // namespace MoFEM
NOSPACE
@ NOSPACE
Definition: definitions.h:175
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MoFEM::BasicMethod::getNinTheLoop
int getNinTheLoop() const
get number of evaluated element in the loop
Definition: LoopMethods.hpp:251
MoFEM::KspMethod::copyKsp
MoFEMErrorCode copyKsp(const KspMethod &ksp)
copy data form another method
Definition: LoopMethods.cpp:62
MoFEM::FEMethod::getDataFieldEntsPtr
boost::shared_ptr< FieldEntity_vector_view > & getDataFieldEntsPtr() const
Definition: LoopMethods.hpp:405
MoFEM::DataForcesAndSourcesCore::basesOnEntities
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
Definition: DataStructures.hpp:1036
MoFEM::FieldEntity::getLoBitNumberUId
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:198
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
EntityHandle
MoFEM::ForcesAndSourcesCore::getRowNodesIndices
MoFEMErrorCode getRowNodesIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get row node indices from FENumeredDofEntity_multiIndex
Definition: ForcesAndSourcesCore.cpp:317
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:98
H1
@ H1
continuous field
Definition: definitions.h:177
MoFEM::FEDofEntity
keeps information about indexed dofs for the finite element
Definition: DofsMultiIndices.hpp:281
LASTBASE
@ LASTBASE
Definition: definitions.h:161
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:451
MoFEM::BasicMethod::nInTheLoop
int nInTheLoop
number currently of processed method
Definition: LoopMethods.hpp:242
MoFEM::ForcesAndSourcesCore::getNoFieldRowIndices
MoFEMErrorCode getNoFieldRowIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get col NoField indices
Definition: ForcesAndSourcesCore.cpp:468
MoFEM::BasicMethod::loopSize
int loopSize
local number oe methods to process
Definition: LoopMethods.hpp:247
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:36
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into DataForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:34
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::ForcesAndSourcesCore::getNodesIndices
MoFEMErrorCode getNodesIndices(const std::string &field_name, FieldEntity_vector_view &ents_field, VectorInt &nodes_indices, VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const
get node indices
Definition: ForcesAndSourcesCore.cpp:228
MoFEM::VectorFieldEntities
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
Definition: DataStructures.hpp:43
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition: definitions.h:164
MoFEM::cmp_uid_hi
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
Definition: ForcesAndSourcesCore.cpp:43
MoFEM::ForcesAndSourcesCore::UserDataOperator::getProblemColIndices
MoFEMErrorCode getProblemColIndices(const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
Get col indices.
Definition: ForcesAndSourcesCore.cpp:1524
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPCOL
@ OPCOL
Definition: ForcesAndSourcesCore.hpp:100
NOBASE
@ NOBASE
Definition: definitions.h:151
MoFEM::ForcesAndSourcesCore::getProblemNodesIndices
MoFEMErrorCode getProblemNodesIndices(const std::string &field_name, const NumeredDofEntity_multiIndex &dofs, VectorInt &nodes_indices) const
get indices of nodal indices which are declared for problem but not this particular element
Definition: ForcesAndSourcesCore.cpp:492
L2
@ L2
field with C-1 continuity
Definition: definitions.h:180
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1597
MoFEM::ForcesAndSourcesCore::calBernsteinBezierBaseFunctionsOnElement
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
Definition: ForcesAndSourcesCore.cpp:1087
MoFEM::ForcesAndSourcesCore::getNoFieldIndices
MoFEMErrorCode getNoFieldIndices(const std::string &field_name, boost::shared_ptr< FENumeredDofEntity_multiIndex > dofs, VectorInt &nodes_indices) const
get NoField indices
Definition: ForcesAndSourcesCore.cpp:450
MoFEM::Composite_Name_And_Ent_mi_tag
Definition: TagMultiIndices.hpp:95
MoFEM::getMaxOrder
static int getMaxOrder(const ENTMULTIINDEX &multi_index)
Definition: ForcesAndSourcesCore.cpp:117
MoFEM::ForcesAndSourcesCore::getEntityColIndices
MoFEMErrorCode getEntityColIndices(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
Definition: ForcesAndSourcesCore.cpp:435
MoFEM::BasicMethod::fieldsPtr
const Field_multiIndex * fieldsPtr
raw pointer to fields container
Definition: LoopMethods.hpp:269
MoFEM::Types::ApproximationOrder
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:37
MoFEM::ForcesAndSourcesCore::UserDataOperator::OpType
OpType
Controls loop over entities on element.
Definition: ForcesAndSourcesCore.hpp:98
MoFEM::TSMethod::copyTs
MoFEMErrorCode copyTs(const TSMethod &ts)
Copy TS solver data.
Definition: LoopMethods.cpp:112
MoFEM::cmp_uid_lo
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
Definition: ForcesAndSourcesCore.cpp:32
MoFEM::BasicMethod::postProcessHook
boost::function< MoFEMErrorCode()> postProcessHook
Hook function for post-processing.
Definition: LoopMethods.hpp:316
MoFEM::Unique_mi_tag
Definition: TagMultiIndices.hpp:31
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::ForcesAndSourcesCore::postProcess
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: ForcesAndSourcesCore.cpp:1645
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:160
MoFEM::ForcesAndSourcesCore::getMaxRowOrder
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Definition: ForcesAndSourcesCore.cpp:138
MoFEM::FEMethod::getColFieldEnts
auto & getColFieldEnts() const
Definition: LoopMethods.hpp:419
MoFEM::ForcesAndSourcesCore::lastEvaluatedElementEntityType
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
Definition: ForcesAndSourcesCore.hpp:825
MoFEM::ForcesAndSourcesCore::getUserPolynomialBase
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:439
MoFEM::DefaultElementAdjacency::getDefTypeMap
static bool getDefTypeMap(const EntityType fe_type, const EntityType ent_type)
Definition: FEMultiIndices.hpp:314
NumeredDofEntity_multiIndex
multi_index_container< boost::shared_ptr< NumeredDofEntity >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, UId, &NumeredDofEntity::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredDofEntity, unsigned int, &NumeredDofEntity::pArt > >, ordered_non_unique< tag< Idx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::dofIdx > >, ordered_non_unique< tag< PetscGlobalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscGloablDofIdx > >, ordered_non_unique< tag< PetscLocalIdx_mi_tag >, member< NumeredDofEntity, DofIdx, &NumeredDofEntity::petscLocalDofIdx > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredDofEntity::interface_type_DofEntity, EntityHandle, &NumeredDofEntity::getEnt > > > > NumeredDofEntity_multiIndex
MultiIndex container keeps NumeredDofEntity.
Definition: DofsMultiIndices.hpp:483
MoFEM::ForcesAndSourcesCore::getColNodesIndices
MoFEMErrorCode getColNodesIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get col node indices from FENumeredDofEntity_multiIndex
Definition: ForcesAndSourcesCore.cpp:334
MoFEM::DataForcesAndSourcesCore::sPace
std::bitset< LASTSPACE > sPace
spaces on element
Definition: DataStructures.hpp:1029
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:174
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPLAST
@ OPLAST
Definition: ForcesAndSourcesCore.hpp:102
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1613
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEM::Field::getApproxBase
FieldApproximationBase getApproxBase() const
Get approximation base.
Definition: FieldMultiIndices.hpp:274
MoFEM::ForcesAndSourcesCore::operator()
virtual MoFEMErrorCode operator()()
function is run for every finite element
Definition: ForcesAndSourcesCore.cpp:1637
MoFEM::ForcesAndSourcesCore::getRule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
Definition: ForcesAndSourcesCore.cpp:1591
MoFEM::ForcesAndSourcesCore::getFaceTriNodes
MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const
Get nodes on triangles.
Definition: ForcesAndSourcesCore.cpp:878
MoFEM::FEMethod::getRowDofsPtr
auto getRowDofsPtr() const
Definition: LoopMethods.hpp:427
convert.type
type
Definition: convert.py:66
MoFEM::ForcesAndSourcesCore::UserDataOperator::setPtrFE
virtual MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Definition: ForcesAndSourcesCore.cpp:1655
MoFEM::ForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ForcesAndSourcesCore.cpp:1304
MoFEM::BasicMethod::copyBasicMethod
MoFEMErrorCode copyBasicMethod(const BasicMethod &basic)
Copy data from other base method to this base method.
Definition: LoopMethods.cpp:137
MoFEM::ForcesAndSourcesCore::derivedDataOnElement
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
Definition: ForcesAndSourcesCore.hpp:804
MoFEM::FEMethod::getNumberOfNodes
MoFEMErrorCode getNumberOfNodes(int &num_nodes) const
Get number of DOFs on element.
Definition: LoopMethods.cpp:200
MoFEM::ForcesAndSourcesCore::getEntityFieldData
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
Definition: ForcesAndSourcesCore.cpp:722
MoFEM::FieldEntity::getLocalUniqueIdCalculate
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
Definition: FieldEntsMultiIndices.hpp:144
MoFEM::ForcesAndSourcesCore::getEntitySense
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< DataForcesAndSourcesCore::EntData > &data) const
get sense (orientation) of entity
Definition: ForcesAndSourcesCore.cpp:90
MoFEM::ForcesAndSourcesCore::getNoFieldColIndices
MoFEMErrorCode getNoFieldColIndices(DataForcesAndSourcesCore &data, const std::string &field_name) const
get col NoField indices
Definition: ForcesAndSourcesCore.cpp:479
MoFEM::DataForcesAndSourcesCore::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: DataStructures.hpp:1031
MoFEM::FEMethod::getDataFieldEnts
const FieldEntity_vector_view & getDataFieldEnts() const
Definition: LoopMethods.hpp:400
ElectroPhysiology::b
const double b
Definition: ElecPhysOperators.hpp:30
MoFEM::ForcesAndSourcesCore::getProblemNodesRowIndices
MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name, VectorInt &nodes_indices) const
Definition: ForcesAndSourcesCore.cpp:571
MoFEM::Problem::numeredColDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
Definition: ProblemsMultiIndices.hpp:100
lapack_wrap.h
StdRDOperators::order
const int order
Definition: rd_stdOperators.hpp:28
MoFEM::BasicMethod::operatorHook
boost::function< MoFEMErrorCode()> operatorHook
Hook function for operator.
Definition: LoopMethods.hpp:321
MoFEM::get_id_for_max_type
EntityHandle get_id_for_max_type()
Definition: RefEntsMultiIndices.hpp:26
MoFEM::ForcesAndSourcesCore::getNodesFieldData
MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data, const std::string &field_name) const
Get data on nodes.
Definition: ForcesAndSourcesCore.cpp:602
LASTSPACE
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:181
CATCH_OP_ERRORS
#define CATCH_OP_ERRORS(OP)
Definition: ForcesAndSourcesCore.cpp:1280
MoFEM::Field::getSpace
FieldSpace getSpace() const
Get field approximation space.
Definition: FieldMultiIndices.hpp:268
MoFEM::FieldName_mi_tag
MultiIndex Tag for field name.
Definition: TagMultiIndices.hpp:80
SideNumber_multiIndex
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
Definition: RefEntsMultiIndices.hpp:92
MoFEM::BitFieldId_space_mi_tag
Definition: TagMultiIndices.hpp:81
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:154
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEM::FieldEntity::getLoLocalEntityBitNumber
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition: FieldEntsMultiIndices.hpp:230
MoFEM::FEMethod::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Definition: LoopMethods.hpp:517
MoFEM::ForcesAndSourcesCore::getElementPolynomialBase
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:432
MoFEM::interface_RefEntity::getSideNumberPtr
boost::shared_ptr< SideNumber > getSideNumberPtr() const
Get the Side number.
Definition: RefEntsMultiIndices.hpp:592
MoFEM::ForcesAndSourcesCore::calHierarchicalBaseFunctionsOnElement
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
Definition: ForcesAndSourcesCore.cpp:1067
MoFEM::FEMethod::getColDofsPtr
auto getColDofsPtr() const
Definition: LoopMethods.hpp:431
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
Definition: ForcesAndSourcesCore.hpp:101
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:389
cblas.h
MoFEM::ForcesAndSourcesCore::getRuleHook
RuleHookFun getRuleHook
Hook to get rule.
Definition: ForcesAndSourcesCore.hpp:50
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:156
MoFEM::DataForcesAndSourcesCore
data structure for finite element entity
Definition: DataStructures.hpp:52
MoFEM::Field::getId
const BitFieldId & getId() const
Get unique field id.
Definition: FieldMultiIndices.hpp:246
MoFEM::ForcesAndSourcesCore::getNoFieldFieldData
MoFEMErrorCode getNoFieldFieldData(const std::string field_name, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
Definition: ForcesAndSourcesCore.cpp:811
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< DataForcesAndSourcesCore >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:797
MoFEM::BasicMethod::problemPtr
const Problem * problemPtr
raw pointer to problem
Definition: LoopMethods.hpp:267
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:64
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
Definition: ForcesAndSourcesCore.hpp:99
MoFEM::DataForcesAndSourcesCore::basesOnSpaces
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
Definition: DataStructures.hpp:1038
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
MoFEM::ForcesAndSourcesCore::getProblemTypeRowIndices
MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
Definition: ForcesAndSourcesCore.cpp:578
MoFEM::DataForcesAndSourcesCore::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: DataStructures.hpp:1040
MoFEM::ForcesAndSourcesCore::getMaxDataOrder
int getMaxDataOrder() const
Get max order of approximation for data fields.
Definition: ForcesAndSourcesCore.cpp:127
MoFEM::get_id_for_min_type
EntityHandle get_id_for_min_type()
Definition: RefEntsMultiIndices.hpp:31
MoFEM::FieldEntity::getHiBitNumberUId
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
Definition: FieldEntsMultiIndices.hpp:207
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:34
MoFEM::ForcesAndSourcesCore::getEntityDataOrder
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< DataForcesAndSourcesCore::EntData > &data) const
Get the entity data order.
Definition: ForcesAndSourcesCore.cpp:146
MoFEM::ForcesAndSourcesCore::getProblemNodesColIndices
MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name, VectorInt &nodes_indices) const
Definition: ForcesAndSourcesCore.cpp:585
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:971
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:178
MoFEM::ForcesAndSourcesCore::UserDataOperator::getProblemRowIndices
MoFEMErrorCode getProblemRowIndices(const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
Get row indices.
Definition: ForcesAndSourcesCore.cpp:1507
MoFEM::ForcesAndSourcesCore::sidePtrFE
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
Definition: ForcesAndSourcesCore.hpp:843
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
MoFEM::Types::VectorInt
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
MoFEM::SnesMethod::copySnes
MoFEMErrorCode copySnes(const SnesMethod &snes)
Copy snes data.
Definition: LoopMethods.cpp:86
MoFEM::ForcesAndSourcesCore::getMaxColOrder
int getMaxColOrder() const
Get max order of approximation for field in columns.
Definition: ForcesAndSourcesCore.cpp:142
MoFEM::DataForcesAndSourcesCore::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: DataStructures.hpp:1034
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: DataStructures.hpp:35
MoFEM::DataForcesAndSourcesCore::bAse
std::bitset< LASTBASE > bAse
bases on element
Definition: DataStructures.hpp:1030
MoFEM::Types::BitFieldId
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:150
MoFEM::Field::getNbOfCoeffs
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
Definition: FieldMultiIndices.hpp:290
MoFEM::ForcesAndSourcesCore::createDataOnElement
MoFEMErrorCode createDataOnElement()
Create a entity data on element object.
Definition: ForcesAndSourcesCore.cpp:1256
MoFEM::Problem::numeredFiniteElementsPtr
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
Definition: ProblemsMultiIndices.hpp:102
MoFEM::ForcesAndSourcesCore::getProblemTypeColIndices
MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
Definition: ForcesAndSourcesCore.cpp:592
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:79
MoFEM::ForcesAndSourcesCore::opPtrVector
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
Definition: ForcesAndSourcesCore.hpp:816
MoFEM::DataOperator
base operator to do operations at Gauss Pt. level
Definition: DataOperators.hpp:36
MoFEM::ForcesAndSourcesCore::preProcess
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: ForcesAndSourcesCore.cpp:1629
MoFEM::FieldEntity::getHiLocalEntityBitNumber
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition: FieldEntsMultiIndices.hpp:241
MoFEM::FieldEntity_vector_view
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
Definition: FieldEntsMultiIndices.hpp:449
MoFEM::ForcesAndSourcesCore::setSideFEPtr
MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr)
Set the pointer to face element on the side.
Definition: ForcesAndSourcesCore.cpp:1542
MoFEM::BasicMethod::preProcessHook
boost::function< MoFEMErrorCode()> preProcessHook
Hook function for pre-processing.
Definition: LoopMethods.hpp:311
MoFEM::ForcesAndSourcesCore::ForcesAndSourcesCore
ForcesAndSourcesCore(Interface &m_field)
Definition: ForcesAndSourcesCore.cpp:54
MoFEM::FEMethod::getRowFieldEnts
auto & getRowFieldEnts() const
Definition: LoopMethods.hpp:411
convert.int
int
Definition: convert.py:66
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
MoFEM::ForcesAndSourcesCore::getEntityIndices
MoFEMErrorCode getEntityIndices(DataForcesAndSourcesCore &data, const std::string &field_name, FieldEntity_vector_view &ents_field, const EntityType type_lo, const EntityType type_hi, EXTRACTOR &&extractor) const
Definition: ForcesAndSourcesCore.cpp:351
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:179
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:124
MoFEM::ForcesAndSourcesCore::UserDataOperator::loopSide
MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t dim, const EntityHandle ent_for_side=0)
User call this function to loop over elements on the side of face. This function calls finite element...
Definition: ForcesAndSourcesCore.cpp:1548
MoFEM::Types::VectorDouble
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
MoFEM::ForcesAndSourcesCore::getProblemTypeIndices
MoFEMErrorCode getProblemTypeIndices(const std::string &field_name, const NumeredDofEntity_multiIndex &dofs, EntityType type, int side_number, VectorInt &indices) const
get indices by type (generic function) which are declared for problem but not this particular element
Definition: ForcesAndSourcesCore.cpp:536
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
MoFEM::ForcesAndSourcesCore::getEntityRowIndices
MoFEMErrorCode getEntityRowIndices(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
Definition: ForcesAndSourcesCore.cpp:420
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:552
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
MoFEM::Types::UId
uint128_t UId
Unique Id.
Definition: Types.hpp:42
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
MoFEM::FEMethod::feName
std::string feName
Name of finite element.
Definition: LoopMethods.hpp:386
NOFIELD
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:176
MoFEM::ForcesAndSourcesCore::setRuleHook
RuleHookFun setRuleHook
Set function to calculate integration rule.
Definition: ForcesAndSourcesCore.hpp:56
MoFEM::ForcesAndSourcesCore::UserDataOperator::ptrFE
ForcesAndSourcesCore * ptrFE
Definition: ForcesAndSourcesCore.hpp:389
cache
BlockParamData * cache
Definition: multifield_plasticity.cpp:81