v0.13.1
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
21extern "C" {
22#endif
23#include <cblas.h>
24#include <lapack_wrap.h>
25// #include <gm_rule.h>
26#ifdef __cplusplus
27}
28#endif
29
30namespace MoFEM {
31
32static 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
43static 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<EntitiesFieldData>(MBENTITYSET), // NOSPACE,
61 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // NOFIELD
62 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // H1
63 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HCURL
64 boost::make_shared<EntitiesFieldData>(MBENTITYSET), // HDIV
65 boost::make_shared<EntitiesFieldData>(MBENTITYSET) // L2
66
67 },
68 derivedDataOnElement{
69
70 nullptr,
71 boost::make_shared<DerivedEntitiesFieldData>(
72 dataOnElement[NOFIELD]), // NOFIELD
73 boost::make_shared<DerivedEntitiesFieldData>(dataOnElement[H1]), // H1
74 boost::make_shared<DerivedEntitiesFieldData>(
75 dataOnElement[HCURL]), // HCURL
76 boost::make_shared<DerivedEntitiesFieldData>(
77 dataOnElement[HDIV]), // HDIV
78 boost::make_shared<DerivedEntitiesFieldData>(dataOnElement[L2]) // L2
79
80 },
81 dataNoField(*dataOnElement[NOFIELD].get()),
82 dataH1(*dataOnElement[H1].get()), dataHcurl(*dataOnElement[HCURL].get()),
83 dataHdiv(*dataOnElement[HDIV].get()), dataL2(*dataOnElement[L2].get()),
84 lastEvaluatedElementEntityType(MBMAXTYPE), sidePtrFE(nullptr),
85 refinePtrFE(nullptr) {
86
87 dataOnElement[NOSPACE]->dataOnEntities[MBENTITYSET].push_back(
89
90 dataOnElement[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
92 derivedDataOnElement[NOFIELD]->dataOnEntities[MBENTITYSET].push_back(
94}
95
96// ** Sense **
97
99 const EntityType type,
100 boost::ptr_vector<EntitiesFieldData::EntData> &data) const {
102
103 auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
104 auto sit = side_table.lower_bound(get_id_for_min_type(type));
105 if (sit != side_table.end()) {
106 auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
107 for (; sit != hi_sit; ++sit) {
108 if (const auto side_number = (*sit)->side_number; side_number >= 0) {
109 const int brother_side_number = (*sit)->brother_side_number;
110 const int sense = (*sit)->sense;
111
112 data[side_number].getSense() = sense;
113 if (brother_side_number != -1)
114 data[brother_side_number].getSense() = sense;
115 }
116 }
117 }
119}
120
121// ** Order **
122
123template <typename ENTMULTIINDEX>
124static inline int getMaxOrder(const ENTMULTIINDEX &multi_index) {
125 int max_order = 0;
126 for (auto ent_field_weak_ptr : multi_index)
127 if (auto e = ent_field_weak_ptr.lock()) {
128 const int order = e->getMaxOrder();
129 max_order = (max_order < order) ? order : max_order;
130 }
131 return max_order;
132}
133
135 int max_order = 0;
136 for (auto e : getDataFieldEnts()) {
137 if (auto ptr = e.lock()) {
138 const int order = ptr->getMaxOrder();
139 max_order = (max_order < order) ? order : max_order;
140 }
141 }
142 return max_order;
143}
144
147}
148
151}
152
154 const EntityType type, const FieldSpace space,
155 boost::ptr_vector<EntitiesFieldData::EntData> &data) const {
157
158 auto set_order = [&]() {
160 auto &side_table = numeredEntFiniteElementPtr->getSideNumberTable();
161
162 for (unsigned int s = 0; s != data.size(); ++s)
163 data[s].getOrder() = 0;
164
165 const FieldEntity_vector_view &data_field_ent = getDataFieldEnts();
166
167 for (auto r = fieldsPtr->get<BitFieldId_space_mi_tag>().equal_range(space);
168 r.first != r.second; ++r.first) {
169
170 const auto field_bit_number = (*r.first)->getBitNumber();
171 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
172 field_bit_number, get_id_for_min_type(type));
173 auto lo = std::lower_bound(data_field_ent.begin(), data_field_ent.end(),
174 lo_uid, cmp_uid_lo);
175 if (lo != data_field_ent.end()) {
176 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
177 field_bit_number, get_id_for_max_type(type));
178 auto hi =
179 std::upper_bound(lo, data_field_ent.end(), hi_uid, cmp_uid_hi);
180 for (; lo != hi; ++lo) {
181
182 if (auto ptr = lo->lock()) {
183
184 auto &e = *ptr;
185 auto sit = side_table.find(e.getEnt());
186 if (sit != side_table.end()) {
187 auto &side = *sit;
188 if (const auto side_number = side->side_number; side_number >= 0) {
189 ApproximationOrder ent_order = e.getMaxOrder();
190 auto &dat = data[side_number];
191 dat.getOrder() =
192 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
193 }
194 } else
195 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
196 "Entity on side of the element not found");
197 }
198 }
199 }
200 }
201
203 };
204
205 auto set_order_on_brother = [&]() {
207 auto &side_table =
208 numeredEntFiniteElementPtr->getSideNumberTable().get<0>();
209 auto sit = side_table.lower_bound(get_id_for_min_type(type));
210 if (sit != side_table.end()) {
211 auto hi_sit = side_table.upper_bound(get_id_for_max_type(type));
212 for (; sit != hi_sit; ++sit) {
213 const int brother_side_number = (*sit)->brother_side_number;
214 if (brother_side_number != -1) {
215 const int side_number = (*sit)->side_number;
216 data[brother_side_number].getOrder() = data[side_number].getOrder();
217 }
218 }
219 }
221 };
222
223 CHKERR set_order();
224 CHKERR set_order_on_brother();
225
227}
228
229// ** Indices **
230
231template <typename EXTRACTOR>
233 const int bit_number, FieldEntity_vector_view &ents_field,
234 VectorInt &nodes_indices, VectorInt &local_nodes_indices,
235 EXTRACTOR &&extractor) const {
237
238 auto field_it = fieldsPtr->get<BitFieldId_mi_tag>().find(
239 BitFieldId().set(bit_number - 1));
240 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
241
242#ifndef NDEBUG
243 if ((*field_it)->getBitNumber() != bit_number)
244 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
245#endif
246 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
247 bit_number, get_id_for_min_type<MBVERTEX>());
248 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
249 cmp_uid_lo);
250 if (lo != ents_field.end()) {
251 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
252 bit_number, get_id_for_max_type<MBVERTEX>());
253 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
254
255 const int num_nodes = getNumberOfNodes();
256 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
257 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
258
259 int nb_dofs = 0;
260 for (auto it = lo; it != hi; ++it) {
261 if (auto e = it->lock()) {
262 if (auto cache = extractor(e).lock()) {
263 if (cache->loHi[0] != cache->loHi[1]) {
264 nb_dofs += std::distance(cache->loHi[0], cache->loHi[1]);
265 break;
266 }
267 }
268 }
269 }
270
271 if (nb_dofs) {
272 nodes_indices.resize(max_nb_dofs, false);
273 local_nodes_indices.resize(max_nb_dofs, false);
274 } else {
275 nodes_indices.resize(0, false);
276 local_nodes_indices.resize(0, false);
277 }
278
279 if (nb_dofs != max_nb_dofs) {
280 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
281 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(), -1);
282 }
283
284 for (auto it = lo; it != hi; ++it) {
285 if (auto e = it->lock()) {
286 auto side_ptr = e->getSideNumberPtr();
287 if (const auto side_number = side_ptr->side_number; side_number >= 0) {
288 const auto brother_side_number = side_ptr->brother_side_number;
289 if (auto cache = extractor(e).lock()) {
290 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
291 auto &dof = **dit;
292 const int idx = dof.getPetscGlobalDofIdx();
293 const int local_idx = dof.getPetscLocalDofIdx();
294 const int pos =
295 side_number * nb_dofs_on_vert + dof.getDofCoeffIdx();
296 nodes_indices[pos] = idx;
297 local_nodes_indices[pos] = local_idx;
298 if (brother_side_number != -1) {
299 const int elem_idx = brother_side_number * nb_dofs_on_vert +
300 (*dit)->getDofCoeffIdx();
301 nodes_indices[elem_idx] = idx;
302 local_nodes_indices[elem_idx] = local_idx;
303 }
304 }
305 }
306 }
307 }
308 }
309 } else {
310 nodes_indices.resize(0, false);
311 local_nodes_indices.resize(0, false);
312 }
313
314 } else {
315 nodes_indices.resize(0, false);
316 local_nodes_indices.resize(0, false);
317 }
318
320}
321
324 const int bit_number) const {
325
326 struct Extractor {
327 boost::weak_ptr<EntityCacheNumeredDofs>
328 operator()(boost::shared_ptr<FieldEntity> &e) {
329 return e->entityCacheRowDofs;
330 }
331 };
332
333 return getNodesIndices(bit_number, getRowFieldEnts(),
334 data.dataOnEntities[MBVERTEX][0].getIndices(),
335 data.dataOnEntities[MBVERTEX][0].getLocalIndices(),
336 Extractor());
337}
338
341 const int bit_number) const {
342
343 struct Extractor {
344 boost::weak_ptr<EntityCacheNumeredDofs>
345 operator()(boost::shared_ptr<FieldEntity> &e) {
346 return e->entityCacheColDofs;
347 }
348 };
349
350 return getNodesIndices(bit_number, getColFieldEnts(),
351 data.dataOnEntities[MBVERTEX][0].getIndices(),
352 data.dataOnEntities[MBVERTEX][0].getLocalIndices(),
353 Extractor());
354}
355
356template <typename EXTRACTOR>
358 EntitiesFieldData &data, const int bit_number,
359 FieldEntity_vector_view &ents_field, const EntityType type_lo,
360 const EntityType type_hi, EXTRACTOR &&extractor) const {
362
363 for (EntityType t = type_lo; t != type_hi; ++t) {
364 for (auto &dat : data.dataOnEntities[t]) {
365 dat.getIndices().resize(0, false);
366 dat.getLocalIndices().resize(0, false);
367 }
368 }
369
370 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
371 bit_number, get_id_for_min_type(type_lo));
372 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
373 cmp_uid_lo);
374 if (lo != ents_field.end()) {
375 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
376 bit_number, get_id_for_max_type(type_hi));
377 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid, cmp_uid_hi);
378 if (lo != hi) {
379
380 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
381
382 for (auto it = lo; it != hi; ++it)
383 if (auto e = it->lock()) {
384
385 const EntityType type = e->getEntType();
386 auto side_ptr = e->getSideNumberPtr();
387 if (const auto side = side_ptr->side_number; side >= 0) {
388 const auto nb_dofs_on_ent = e->getNbDofsOnEnt();
389 const auto brother_side = side_ptr->brother_side_number;
390 auto &dat = data.dataOnEntities[type][side];
391 auto &ent_field_indices = dat.getIndices();
392 auto &ent_field_local_indices = dat.getLocalIndices();
393
394 ent_field_indices.resize(nb_dofs_on_ent, false);
395 ent_field_local_indices.resize(nb_dofs_on_ent, false);
396 std::fill(ent_field_indices.data().begin(),
397 ent_field_indices.data().end(), -1);
398 std::fill(ent_field_local_indices.data().begin(),
399 ent_field_local_indices.data().end(), -1);
400
401 if (auto cache = extractor(e).lock()) {
402 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
403 const int idx = (*dit)->getEntDofIdx();
404 ent_field_indices[idx] = (*dit)->getPetscGlobalDofIdx();
405 ent_field_local_indices[idx] = (*dit)->getPetscLocalDofIdx();
406 }
407 }
408
409 if (brother_side != -1) {
410 auto &dat_brother = data.dataOnEntities[type][brother_side];
411 dat_brother.getIndices().resize(nb_dofs_on_ent, false);
412 dat_brother.getLocalIndices().resize(nb_dofs_on_ent, false);
413 noalias(dat_brother.getIndices()) = dat.getIndices();
414 noalias(dat_brother.getLocalIndices()) = dat.getLocalIndices();
415 }
416 }
417 }
418 }
419 }
420
422}
423
425 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
426 const EntityType type_hi) const {
427
428 struct Extractor {
429 boost::weak_ptr<EntityCacheNumeredDofs>
430 operator()(boost::shared_ptr<FieldEntity> &e) {
431 return e->entityCacheRowDofs;
432 }
433 };
434
435 return getEntityIndices(data, bit_number, getRowFieldEnts(), type_lo, type_hi,
436 Extractor());
437}
438
440 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
441 const EntityType type_hi) const {
442
443 struct Extractor {
444 boost::weak_ptr<EntityCacheNumeredDofs>
445 operator()(boost::shared_ptr<FieldEntity> &e) {
446 return e->entityCacheColDofs;
447 }
448 };
449
450 return getEntityIndices(data, bit_number, getColFieldEnts(), type_lo, type_hi,
451 Extractor());
452}
453
455 const int bit_number, boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
456 VectorInt &indices) const {
458 auto dit = dofs->get<Unique_mi_tag>().lower_bound(
460 auto hi_dit = dofs->get<Unique_mi_tag>().upper_bound(
462 indices.resize(std::distance(dit, hi_dit));
463 for (; dit != hi_dit; dit++) {
464 int idx = (*dit)->getPetscGlobalDofIdx();
465 indices[(*dit)->getDofCoeffIdx()] = idx;
466 }
468}
469
472 const int bit_number) const {
474#ifndef NDEBUG
475 if (data.dataOnEntities[MBENTITYSET].size() == 0) {
476 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
477 "data.dataOnEntities[MBENTITYSET] is empty");
478 }
479#endif
481 data.dataOnEntities[MBENTITYSET][0].getIndices());
483}
484
487 const int bit_number) const {
489#ifndef NDEBUG
490 if (data.dataOnEntities[MBENTITYSET].size() == 0) {
491 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
492 "data.dataOnEntities[MBENTITYSET] is empty");
493 }
494#endif
496 data.dataOnEntities[MBENTITYSET][0].getIndices());
498}
499
500// ** Indices from problem **
501
503 const std::string &field_name, const NumeredDofEntity_multiIndex &dofs,
504 VectorInt &nodes_indices) const {
506
507 const Field *field_struture = mField.get_field_structure(field_name);
508 if (field_struture->getSpace() == H1) {
509
510 const int num_nodes = getNumberOfNodes();
511 nodes_indices.resize(field_struture->getNbOfCoeffs() * num_nodes, false);
512 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
513
514 auto &side_table = const_cast<SideNumber_multiIndex &>(
515 numeredEntFiniteElementPtr->getSideNumberTable());
516 auto siit = side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 0));
517 auto hi_siit =
518 side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 10000));
519
520 int nn = 0;
521 for (; siit != hi_siit; siit++, nn++) {
522 if (siit->get()->side_number >= 0) {
523 auto bit_number = mField.get_field_bit_number(field_name);
524 const EntityHandle ent = siit->get()->ent;
525 auto dit = dofs.get<Unique_mi_tag>().lower_bound(
527 auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
529 for (; dit != hi_dit; dit++) {
530 nodes_indices[siit->get()->side_number * (*dit)->getNbOfCoeffs() +
531 (*dit)->getDofCoeffIdx()] =
532 (*dit)->getPetscGlobalDofIdx();
533 }
534 }
535 }
536 } else {
537 nodes_indices.resize(0, false);
538 }
539
541}
542
544 const std::string &field_name, const NumeredDofEntity_multiIndex &dofs,
545 EntityType type, int side_number, VectorInt &indices) const {
547
548 indices.resize(0);
549
550 auto &side_table = const_cast<SideNumber_multiIndex &>(
551 numeredEntFiniteElementPtr->getSideNumberTable());
552 auto siit =
553 side_table.get<1>().lower_bound(boost::make_tuple(type, side_number));
554 auto hi_siit =
555 side_table.get<1>().upper_bound(boost::make_tuple(type, side_number));
556
557 for (; siit != hi_siit; siit++) {
558 if (siit->get()->side_number >= 0) {
559
560 const EntityHandle ent = siit->get()->ent;
561 auto bit_number = mField.get_field_bit_number(field_name);
562 auto dit = dofs.get<Unique_mi_tag>().lower_bound(
564 auto hi_dit = dofs.get<Unique_mi_tag>().upper_bound(
566 indices.resize(std::distance(dit, hi_dit));
567 for (; dit != hi_dit; dit++) {
568
569 indices[(*dit)->getEntDofIdx()] = (*dit)->getPetscGlobalDofIdx();
570 }
571 }
572 }
573
575}
576
578 const std::string &field_name, VectorInt &nodes_indices) const {
580 nodes_indices);
581}
582
585 EntityType type, int side_number,
586 VectorInt &indices) const {
588 type, side_number, indices);
589}
590
592 const std::string &field_name, VectorInt &nodes_indices) const {
594 nodes_indices);
595}
596
599 EntityType type, int side_number,
600 VectorInt &indices) const {
602 type, side_number, indices);
603}
604
605// ** Data **
606
608ForcesAndSurcesCore::getBitRefLevelOnData() {
610
611 // for (int s = H1; s != LASTSPACE; ++s) {
612 // dataOnElement[s]->dataOnEntities[MBENTITYSET].resize(0);
613 // }
614
615 for (auto &data : dataOnElement) {
616 if (data) {
617 for (auto &dat : data->dataOnEntities) {
618 for (auto &ent_dat : dat) {
619 ent_dat.getEntDataBitRefLevel().reset();
620 }
621 }
622 }
623 }
624
625 auto &field_ents = getDataFieldEnts();
626 for (auto it : field_ents) {
627 if (auto e = it.lock()) {
628 const FieldSpace space = e->getSpace();
629 if (space > NOFIELD) {
630 const EntityType type = e->getEntType();
631 const signed char side =
632 type == MBVERTEX ? 0 : e->getSideNumberPtr()->side_number;
633 if (side >= 0) {
634 if (auto &data = dataOnElement[space]) {
635 auto &dat = data->dataOnEntities[type][side];
636 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
637 }
638 }
639 } else {
640 if (auto &data = dataOnElement[NOFIELD]) {
641 auto &dat = data->dataOnEntities[MBENTITYSET][0];
642 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
643 }
644 }
645 }
646 }
647
649};
650
653 const int bit_number) const {
654
655 auto get_nodes_field_data = [&](VectorDouble &nodes_data,
656 VectorFieldEntities &field_entities,
657 VectorDofs &nodes_dofs, FieldSpace &space,
659 VectorInt &bb_node_order) {
661
662 nodes_data.resize(0, false);
663 nodes_dofs.resize(0, false);
664 field_entities.resize(0, false);
665
666 auto field_it =
667 fieldsPtr->get<BitFieldId_mi_tag>().find(BitFEId().set(bit_number - 1));
668 if (field_it != fieldsPtr->get<BitFieldId_mi_tag>().end()) {
669
670#ifndef NDEBUG
671 if ((*field_it)->getBitNumber() != bit_number)
672 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong bit number");
673#endif
674 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
675 space = (*field_it)->getSpace();
676 base = (*field_it)->getApproxBase();
677
678 auto &field_ents = getDataFieldEnts();
679 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
680 bit_number, get_id_for_min_type<MBVERTEX>());
681 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
682 cmp_uid_lo);
683 if (lo != field_ents.end()) {
684 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
685 bit_number, get_id_for_max_type<MBVERTEX>());
686 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
687 if (lo != hi) {
688
689 int nb_dofs = 0;
690 for (auto it = lo; it != hi; ++it) {
691 if (auto e = it->lock()) {
692 nb_dofs += e->getNbDofsOnEnt();
693 }
694 }
695
696 if (nb_dofs) {
697
698 const int num_nodes = getNumberOfNodes();
699 bb_node_order.resize(num_nodes, false);
700 bb_node_order.clear();
701 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
702 nodes_data.resize(max_nb_dofs, false);
703 nodes_dofs.resize(max_nb_dofs, false);
704 field_entities.resize(num_nodes, false);
705 std::fill(nodes_data.begin(), nodes_data.end(), 0);
706 std::fill(nodes_dofs.begin(), nodes_dofs.end(), nullptr);
707 std::fill(field_entities.begin(), field_entities.end(), nullptr);
708
709 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
710
711 for (auto it = lo; it != hi; ++it) {
712 if (auto e = it->lock()) {
713 const auto &sn = e->getSideNumberPtr();
714 // Some field entities on skeleton can have negative side
715 // number
716 if (const auto side_number = sn->side_number; side_number >= 0) {
717 const int brother_side_number = sn->brother_side_number;
718
719 field_entities[side_number] = e.get();
720 if (brother_side_number != -1) {
721 brother_ents_vec.emplace_back(e);
722 field_entities[side_number] = field_entities[side_number];
723 }
724
725 bb_node_order[side_number] = e->getMaxOrder();
726 int pos = side_number * nb_dofs_on_vert;
727 auto ent_filed_data_vec = e->getEntFieldData();
728 if (auto cache = e->entityCacheDataDofs.lock()) {
729 for (auto dit = cache->loHi[0]; dit != cache->loHi[1];
730 ++dit) {
731 const auto dof_idx = (*dit)->getEntDofIdx();
732 nodes_data[pos + dof_idx] = ent_filed_data_vec[dof_idx];
733 nodes_dofs[pos + dof_idx] =
734 reinterpret_cast<FEDofEntity *>((*dit).get());
735 }
736 }
737 }
738 }
739 }
740
741 for (auto &it : brother_ents_vec) {
742 if (const auto e = it.lock()) {
743 const auto &sn = e->getSideNumberPtr();
744 const int side_number = sn->side_number;
745 const int brother_side_number = sn->brother_side_number;
746 bb_node_order[brother_side_number] = bb_node_order[side_number];
747 int pos = side_number * nb_dofs_on_vert;
748 int brother_pos = brother_side_number * nb_dofs_on_vert;
749 for (int ii = 0; ii != nb_dofs_on_vert; ++ii) {
750 nodes_data[brother_pos] = nodes_data[pos];
751 nodes_dofs[brother_pos] = nodes_dofs[pos];
752 ++pos;
753 ++brother_pos;
754 }
755 }
756 }
757 }
758 }
759 }
760 }
761
763 };
764
765 return get_nodes_field_data(
766 data.dataOnEntities[MBVERTEX][0].getFieldData(),
767 data.dataOnEntities[MBVERTEX][0].getFieldEntities(),
768 data.dataOnEntities[MBVERTEX][0].getFieldDofs(),
769 data.dataOnEntities[MBVERTEX][0].getSpace(),
770 data.dataOnEntities[MBVERTEX][0].getBase(),
771 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder());
772}
773
775 EntitiesFieldData &data, const int bit_number, const EntityType type_lo,
776 const EntityType type_hi) const {
778 for (EntityType t = type_lo; t != type_hi; ++t) {
779 for (auto &dat : data.dataOnEntities[t]) {
780 dat.getOrder() = 0;
781 dat.getBase() = NOBASE;
782 dat.getSpace() = NOSPACE;
783 dat.getFieldData().resize(0, false);
784 dat.getFieldDofs().resize(0, false);
785 dat.getFieldEntities().resize(0, false);
786 }
787 }
788
789 auto &field_ents = getDataFieldEnts();
790 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
791 bit_number, get_id_for_min_type(type_lo));
792 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
793 cmp_uid_lo);
794 if (lo != field_ents.end()) {
795 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
796 bit_number, get_id_for_max_type(type_hi));
797 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
798 if (lo != hi) {
799
800 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
801
802 for (auto it = lo; it != hi; ++it)
803 if (auto e = it->lock()) {
804 auto side_ptr = e->getSideNumberPtr();
805 if (const auto side = side_ptr->side_number; side >= 0) {
806 const EntityType type = e->getEntType();
807 auto &dat = data.dataOnEntities[type][side];
808 auto &ent_field = dat.getFieldEntities();
809 auto &ent_field_dofs = dat.getFieldDofs();
810 auto &ent_field_data = dat.getFieldData();
811
812 const int brother_side = side_ptr->brother_side_number;
813 if (brother_side != -1)
814 brother_ents_vec.emplace_back(e);
815
816 dat.getBase() = e->getApproxBase();
817 dat.getSpace() = e->getSpace();
818 const int ent_order = e->getMaxOrder();
819 dat.getOrder() =
820 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
821
822 auto ent_data = e->getEntFieldData();
823 ent_field_data.resize(ent_data.size(), false);
824 noalias(ent_field_data) = ent_data;
825 ent_field_dofs.resize(ent_data.size(), false);
826 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(), nullptr);
827 ent_field.resize(1, false);
828 ent_field[0] = e.get();
829 if (auto cache = e->entityCacheDataDofs.lock()) {
830 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
831 ent_field_dofs[(*dit)->getEntDofIdx()] =
832 reinterpret_cast<FEDofEntity *>((*dit).get());
833 }
834 }
835 }
836 }
837
838 for (auto &it : brother_ents_vec) {
839 if (const auto e = it.lock()) {
840 const EntityType type = e->getEntType();
841 const int side = e->getSideNumberPtr()->side_number;
842 const int brother_side = e->getSideNumberPtr()->brother_side_number;
843 auto &dat = data.dataOnEntities[type][side];
844 auto &dat_brother = data.dataOnEntities[type][brother_side];
845 dat_brother.getBase() = dat.getBase();
846 dat_brother.getSpace() = dat.getSpace();
847 dat_brother.getOrder() = dat.getOrder();
848 dat_brother.getFieldData() = dat.getFieldData();
849 dat_brother.getFieldDofs() = dat.getFieldDofs();
850 dat_brother.getFieldEntities() = dat.getFieldEntities();
851 }
852 }
853 }
854 }
855
857}
858
860 const int bit_number, VectorDouble &ent_field_data,
861 VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const {
862
864
865 ent_field_data.resize(0, false);
866 ent_field_dofs.resize(0, false);
867 ent_field.resize(0, false);
868
869 auto &field_ents = getDataFieldEnts();
870 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
871 bit_number, get_id_for_min_type<MBVERTEX>());
872 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
873 cmp_uid_lo);
874 if (lo != field_ents.end()) {
875
876 ent_field.resize(field_ents.size(), false);
877 std::fill(ent_field.begin(), ent_field.end(), nullptr);
878
879 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
880 bit_number, get_id_for_max_type<MBENTITYSET>());
881 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
882 if (lo != hi) {
883
884 int side = 0;
885 for (auto it = lo; it != hi; ++it, ++side)
886 if (auto e = it->lock()) {
887
888 const auto size = e->getNbDofsOnEnt();
889 ent_field_data.resize(size, false);
890 ent_field_dofs.resize(size, false);
891 ent_field[side] = e.get();
892 noalias(ent_field_data) = e->getEntFieldData();
893
894 if (auto cache = e->entityCacheDataDofs.lock()) {
895 for (auto dit = cache->loHi[0]; dit != cache->loHi[1]; ++dit) {
896 ent_field_dofs[(*dit)->getEntDofIdx()] =
897 reinterpret_cast<FEDofEntity *>((*dit).get());
898 }
899 }
900 }
901 }
902 }
903
905}
906
909 const int bit_number) const {
911 if (data.dataOnEntities[MBENTITYSET].size() == 0)
912 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
913 "No space to insert data");
914
916 bit_number, data.dataOnEntities[MBENTITYSET][0].getFieldData(),
917 data.dataOnEntities[MBENTITYSET][0].getFieldDofs(),
918 data.dataOnEntities[MBENTITYSET][0].getFieldEntities());
920}
921
922// ** Face **
923
927 auto &face_nodes = data.facesNodes;
928 auto &face_nodes_order = data.facesNodesOrder;
929 auto &side_table = const_cast<SideNumber_multiIndex &>(
930 numeredEntFiniteElementPtr->getSideNumberTable());
931 const auto ent = numeredEntFiniteElementPtr->getEnt();
932 const auto type = numeredEntFiniteElementPtr->getEntType();
933 const auto nb_faces = CN::NumSubEntities(type, 2);
934 const EntityHandle *conn_ele;
935 int num_nodes_ele;
936 CHKERR mField.get_moab().get_connectivity(ent, conn_ele, num_nodes_ele, true);
937 auto side_ptr_it = side_table.get<1>().lower_bound(
938 boost::make_tuple(CN::TypeDimensionMap[2].first, 0));
939 auto hi_side_ptr_it = side_table.get<1>().upper_bound(
940 boost::make_tuple(CN::TypeDimensionMap[2].second, 100));
941
942 for (; side_ptr_it != hi_side_ptr_it; ++side_ptr_it) {
943 const auto side = (*side_ptr_it)->side_number;
944 const auto sense = (*side_ptr_it)->sense;
945 const auto offset = (*side_ptr_it)->offset;
946 const auto face_ent = (*side_ptr_it)->ent;
947
948 EntityType face_type;
949 int nb_nodes_face;
950 auto face_indices =
951 CN::SubEntityVertexIndices(type, 2, side, face_type, nb_nodes_face);
952 face_nodes.resize(nb_faces, nb_nodes_face);
953 face_nodes_order.resize(nb_faces, nb_nodes_face);
954
955 if (sense == 1)
956 for (int n = 0; n != nb_nodes_face; ++n)
957 face_nodes_order(side, n) = (n + offset) % nb_nodes_face;
958 else
959 for (int n = 0; n != nb_nodes_face; ++n)
960 face_nodes_order(side, n) =
961 (nb_nodes_face - (n - offset) % nb_nodes_face) % nb_nodes_face;
962
963 for (int n = 0; n != nb_nodes_face; ++n)
964 face_nodes(side, n) = face_indices[face_nodes_order(side, n)];
965
966#ifndef NDEBUG
967 auto check = [&]() {
969 const EntityHandle *conn_face;
970 // int nb_nodes_face;
971 CHKERR mField.get_moab().get_connectivity(face_ent, conn_face,
972 nb_nodes_face, true);
973 face_nodes.resize(nb_faces, nb_nodes_face);
974 for (int nn = 0; nn != nb_nodes_face; ++nn) {
975 if (face_nodes(side, nn) !=
976 std::distance(
977 conn_ele,
978 std::find(conn_ele, &conn_ele[num_nodes_ele], conn_face[nn]))) {
979 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
980 "Wrong face numeration");
981 }
982 }
984 };
985 CHKERR check();
986#endif
987 }
988
990}
991
992// ** Space and Base **
993
995 EntitiesFieldData &data) const {
997
998 if (nInTheLoop == 0) {
999 data.sPace.reset();
1000 data.bAse.reset();
1001 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
1002 data.spacesOnEntities[t].reset();
1003 data.basesOnEntities[t].reset();
1004 }
1005 for (int s = 0; s != LASTSPACE; ++s) {
1006 data.basesOnSpaces[s].reset();
1007 }
1008 }
1009
1010 auto fe_type = numeredEntFiniteElementPtr->getEntType();
1011
1012 if (getDataFieldEntsPtr())
1013 for (auto e : getDataFieldEnts()) {
1014 if (auto ptr = e.lock()) {
1015 // get data from entity
1016 const EntityType type = ptr->getEntType();
1018 const FieldSpace space = ptr->getSpace();
1019 const FieldApproximationBase approx = ptr->getApproxBase();
1020 // set data
1021 data.sPace.set(space);
1022 data.bAse.set(approx);
1023 data.spacesOnEntities[type].set(space);
1024 data.basesOnEntities[type].set(approx);
1025 data.basesOnSpaces[space].set(approx);
1026 }
1027 }
1028 }
1029 else
1030 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1031 "data fields ents not allocated on element");
1032
1034}
1035
1037 const FieldApproximationBase b) {
1039 if (dataOnElement[H1]->bAse.test(b)) {
1040 switch (static_cast<FieldApproximationBase>(b)) {
1041 case NOBASE:
1042 break;
1044 // BERNSTEIN_BEZIER_BASE is not hierarchical base
1045 break;
1050 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1051 "Functions genrating approximation base not defined");
1052
1053 for (int space = H1; space != LASTSPACE; ++space) {
1054 if (dataOnElement[H1]->sPace.test(space) &&
1055 dataOnElement[H1]->bAse.test(b) &&
1056 dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1057 CHKERR getElementPolynomialBase()->getValue(
1058 gaussPts,
1059 boost::make_shared<EntPolynomialBaseCtx>(
1060 *dataOnElement[space], static_cast<FieldSpace>(space),
1061 static_cast<FieldApproximationBase>(b), NOBASE));
1062 }
1063 }
1064 break;
1065 case USER_BASE:
1066 if (!getUserPolynomialBase())
1067 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1068 "Functions genrating user approximation base not defined");
1069
1070 for (int space = H1; space != LASTSPACE; ++space)
1071 if (dataOnElement[H1]->sPace.test(space) &&
1072 dataOnElement[H1]->bAse.test(b) &&
1073 dataOnElement[H1]->basesOnSpaces[space].test(b)) {
1074 CHKERR getUserPolynomialBase()->getValue(
1075 gaussPts,
1076 boost::make_shared<EntPolynomialBaseCtx>(
1077 *dataOnElement[space], static_cast<FieldSpace>(space),
1078 static_cast<FieldApproximationBase>(b), NOBASE));
1079 }
1080 break;
1081 default:
1083 "Base <%s> not yet implemented",
1085 }
1086 }
1088}
1089
1092 /// Use the some node base. Node base is usually used for construction other
1093 /// bases.
1094 for (int space = HCURL; space != LASTSPACE; ++space) {
1095 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE) =
1096 dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getNSharedPtr(NOBASE);
1097 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1098 NOBASE) =
1099 dataOnElement[H1]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1100 NOBASE);
1101 }
1102 for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
1104 static_cast<FieldApproximationBase>(b));
1105 }
1107}
1108
1112
1113 const auto ele_type = numeredEntFiniteElementPtr->getEntType();
1114
1115 auto get_nodal_base_data = [&](EntitiesFieldData &data, auto field_ptr) {
1117 auto &space = data.dataOnEntities[MBVERTEX][0].getSpace();
1118 auto &base = data.dataOnEntities[MBVERTEX][0].getBase();
1119 auto &bb_node_order = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder();
1120
1121 auto &field_ents = getDataFieldEnts();
1122 auto bit_number = field_ptr->getBitNumber();
1123 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1124 bit_number, get_id_for_min_type<MBVERTEX>());
1125 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1126 cmp_uid_lo);
1127 if (lo != field_ents.end()) {
1128 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1129 bit_number, get_id_for_max_type<MBVERTEX>());
1130 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1131 if (lo != hi) {
1132
1133 for (auto it = lo; it != hi; ++it)
1134 if (auto first_e = it->lock()) {
1135 space = first_e->getSpace();
1136 base = first_e->getApproxBase();
1137 const int num_nodes = getNumberOfNodes();
1138 bb_node_order.resize(num_nodes, false);
1139 bb_node_order.clear();
1140 const int nb_dof_idx = first_e->getNbOfCoeffs();
1141
1142 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1143
1144 for (; it != hi; ++it) {
1145 if (auto e = it->lock()) {
1146 const auto &sn = e->getSideNumberPtr();
1147 const int side_number = sn->side_number;
1148 const int brother_side_number = sn->brother_side_number;
1149 if (brother_side_number != -1)
1150 brother_ents_vec.emplace_back(e);
1151 bb_node_order[side_number] = e->getMaxOrder();
1152 }
1153 }
1154
1155 for (auto &it : brother_ents_vec) {
1156 if (const auto e = it.lock()) {
1157 const auto &sn = e->getSideNumberPtr();
1158 const int side_number = sn->side_number;
1159 const int brother_side_number = sn->brother_side_number;
1160 bb_node_order[brother_side_number] = bb_node_order[side_number];
1161 }
1162 }
1163
1164 break;
1165 }
1166 }
1167 }
1168
1170 };
1171
1172 auto get_entity_base_data = [&](EntitiesFieldData &data, auto field_ptr,
1173 const EntityType type_lo,
1174 const EntityType type_hi) {
1176 for (EntityType t = MBEDGE; t != MBPOLYHEDRON; ++t) {
1177 for (auto &dat : data.dataOnEntities[t]) {
1178 dat.getOrder() = 0;
1179 dat.getBase() = NOBASE;
1180 dat.getSpace() = NOSPACE;
1181 dat.getFieldData().resize(0, false);
1182 dat.getFieldDofs().resize(0, false);
1183 }
1184 }
1185
1186 auto &field_ents = getDataFieldEnts();
1187 auto bit_number = field_ptr->getBitNumber();
1188 const auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1189 bit_number, get_id_for_min_type(type_lo));
1190 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1191 cmp_uid_lo);
1192 if (lo != field_ents.end()) {
1193 const auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1194 bit_number, get_id_for_max_type(type_hi));
1195 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid, cmp_uid_hi);
1196
1197 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1198 for (; lo != hi; ++lo) {
1199 if (auto e = lo->lock()) {
1200 if (auto cache = e->entityCacheDataDofs.lock()) {
1201 if (cache->loHi[0] != cache->loHi[1]) {
1202 if (const auto side = e->getSideNumberPtr()->side_number;
1203 side >= 0) {
1204 const EntityType type = e->getEntType();
1205 auto &dat = data.dataOnEntities[type][side];
1206 const int brother_side =
1207 e->getSideNumberPtr()->brother_side_number;
1208 if (brother_side != -1)
1209 brother_ents_vec.emplace_back(e);
1210 dat.getBase() = e->getApproxBase();
1211 dat.getSpace() = e->getSpace();
1212 const auto ent_order = e->getMaxOrder();
1213 dat.getOrder() =
1214 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
1215 }
1216 }
1217 }
1218 }
1219 }
1220
1221 for (auto &ent_ptr : brother_ents_vec) {
1222 if (auto e = ent_ptr.lock()) {
1223 const EntityType type = e->getEntType();
1224 const int side = e->getSideNumberPtr()->side_number;
1225 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1226 auto &dat = data.dataOnEntities[type][side];
1227 auto &dat_brother = data.dataOnEntities[type][brother_side];
1228 dat_brother.getBase() = dat.getBase();
1229 dat_brother.getSpace() = dat.getSpace();
1230 dat_brother.getOrder() = dat.getOrder();
1231 }
1232 }
1233 }
1235 };
1236
1237 for (auto &e : getDataFieldEnts()) {
1238 if (auto ent_data_ptr = e.lock()) {
1239 if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1240 auto space = ent_data_ptr->getSpace();
1241 for (EntityType t = MBVERTEX; t != MBPOLYHEDRON; ++t) {
1242 for (auto &dat : (*dataOnElement[space]).dataOnEntities[t]) {
1243 for (auto &ptr : dat.getBBAlphaIndicesByOrderArray())
1244 ptr.reset();
1245 for (auto &ptr : dat.getBBNByOrderArray())
1246 ptr.reset();
1247 for (auto &ptr : dat.getBBDiffNByOrderArray())
1248 ptr.reset();
1249 }
1250 }
1251 }
1252 }
1253 }
1254
1255 auto check_space = [&](const auto space) {
1256 switch (space) {
1257 case H1:
1258 for (auto t = MBVERTEX; t <= ele_type; ++t) {
1259 if (dataOnElement[H1]->spacesOnEntities[t].test(H1))
1260 return true;
1261 }
1262 return false;
1263 case HCURL:
1264 for (auto t = MBEDGE; t <= ele_type; ++t) {
1265 if (dataOnElement[HCURL]->spacesOnEntities[t].test(HCURL))
1266 return true;
1267 }
1268 return false;
1269 case HDIV:
1270 for (auto t = MBTRI; t <= ele_type; ++t) {
1271 if (dataOnElement[HDIV]->spacesOnEntities[t].test(HDIV))
1272 return true;
1273 }
1274 return false;
1275 case L2:
1276 return dataOnElement[L2]->spacesOnEntities[ele_type].test(L2);
1277 break;
1278 default:
1279 THROW_MESSAGE("Not implemented");
1280 }
1281 };
1282
1283 std::set<string> fields_list;
1284 for (auto &e : getDataFieldEnts()) {
1285 if (auto ent_data_ptr = e.lock()) {
1286 if (ent_data_ptr->getApproxBase() == AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1287 auto field_name = ent_data_ptr->getName();
1288 if (fields_list.find(field_name) == fields_list.end()) {
1289 auto field_ptr = ent_data_ptr->getFieldRawPtr();
1290 auto space = ent_data_ptr->getSpace();
1291 CHKERR get_nodal_base_data(*dataOnElement[space], field_ptr);
1292 CHKERR get_entity_base_data(*dataOnElement[space], field_ptr, MBEDGE,
1293 ele_type);
1294 if (check_space(space)) {
1295 CHKERR getElementPolynomialBase()->getValue(
1296 gaussPts, boost::make_shared<EntPolynomialBaseCtx>(
1297 *dataOnElement[space], field_name,
1298 static_cast<FieldSpace>(space),
1300 fields_list.insert(field_name);
1301 }
1302 }
1303 }
1304 }
1305 }
1306
1308};
1309
1312
1313 // Data on elements for proper spaces
1314 for (int space = H1; space != LASTSPACE; ++space) {
1315 dataOnElement[space]->setElementType(type);
1316 derivedDataOnElement[space]->setElementType(type);
1317 }
1318
1320}
1321
1322#define FUNCTION_NAME_WITH_OP_NAME(OP) \
1323 std::ostringstream ss; \
1324 ss << "(Calling user data operator " \
1325 << boost::typeindex::type_id_runtime(OP).pretty_name() << " rowField " \
1326 << (OP).rowFieldName << " colField " << (OP).colFieldName << ") "
1327
1328#define CATCH_OP_ERRORS(OP) \
1329 catch (MoFEMExceptionInitial const &ex) { \
1330 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1331 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1332 ex.errorCode, PETSC_ERROR_INITIAL, ex.what()); \
1333 } \
1334 catch (MoFEMExceptionRepeat const &ex) { \
1335 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1336 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1337 ex.errorCode, PETSC_ERROR_REPEAT, " "); \
1338 } \
1339 catch (MoFEMException const &ex) { \
1340 FUNCTION_NAME_WITH_OP_NAME(OP) << ex.errorMessage; \
1341 SETERRQ(PETSC_COMM_SELF, ex.errorCode, ss.str().c_str()); \
1342 } \
1343 catch (std::exception const &ex) { \
1344 std::string message("Error: " + std::string(ex.what()) + " at " + \
1345 boost::lexical_cast<std::string>(__LINE__) + " : " + \
1346 std::string(__FILE__) + " in " + \
1347 std::string(PETSC_FUNCTION_NAME)); \
1348 FUNCTION_NAME_WITH_OP_NAME(OP) << message; \
1349 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str()); \
1350 }
1351
1354
1355 using UDO = UserDataOperator;
1356
1357 std::array<std::string, 2> field_name;
1358 std::array<Field *, 3> field_struture;
1359 std::array<int, 2>
1360 field_id; // note the this is field bit number (nor field bit)
1361 std::array<FieldSpace, 2> space;
1362 std::array<FieldApproximationBase, 2> base;
1363
1364 constexpr std::array<UDO::OpType, 2> types = {UDO::OPROW, UDO::OPCOL};
1365 std::array<int, 2> last_eval_field_id = {0, 0};
1366
1367 std::array<boost::shared_ptr<EntitiesFieldData>, 2> op_data;
1368
1369 auto swap_bases = [&](auto &op) {
1371 for (size_t ss = 0; ss != 2; ++ss) {
1372 if (op.getOpType() & types[ss] || op.getOpType() & UDO::OPROWCOL) {
1373 switch (base[ss]) {
1375 CHKERR op_data[ss]->baseSwap(field_name[ss],
1377 default:
1378 break;
1379 }
1380 }
1381 }
1382
1384 };
1385
1386 const EntityType type = numeredEntFiniteElementPtr->getEntType();
1387
1388 // evaluate entity data only, no field specific data provided or known
1389 auto evaluate_op_space = [&](auto &op) {
1391
1392 // reseat all data which all field dependent
1393 dataOnElement[op.sPace]->resetFieldDependentData();
1394 std::fill(last_eval_field_id.begin(), last_eval_field_id.end(), 0);
1395
1396 switch (op.sPace) {
1397 case NOSPACE:
1398 try {
1399 CHKERR op.doWork(
1400 0, MBENTITYSET,
1401 dataOnElement[NOSPACE]->dataOnEntities[MBENTITYSET][0]);
1402 }
1403 CATCH_OP_ERRORS(op);
1404 break;
1405 case NOFIELD:
1406 case H1:
1407 case HCURL:
1408 case HDIV:
1409 case L2:
1410 try {
1411
1412 for (EntityType t = MBVERTEX; t != MBMAXTYPE; ++t) {
1413 for (auto &e : dataOnElement[op.sPace]->dataOnEntities[t]) {
1414 e.getSpace() = op.sPace;
1415 }
1416 }
1417
1418 CHKERR op.opRhs(*dataOnElement[op.sPace], false);
1419 }
1420 CATCH_OP_ERRORS(op);
1421 break;
1422 default:
1423 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1424 "Not implemented for this space", op.sPace);
1425 }
1426
1428 };
1429
1430 // set entity data
1431 auto set_op_entityties_data = [&](auto ss, auto &op) {
1433
1434#ifndef NDEBUG
1435 if ((op.getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
1437 .none()) {
1438 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1439 "no data field < %s > on finite element < %s >",
1440 field_name[ss].c_str(), getFEName().c_str());
1441 }
1442#endif
1443
1444 op_data[ss] =
1445 !ss ? dataOnElement[space[ss]] : derivedDataOnElement[space[ss]];
1446
1447 for (auto &data : op_data[ss]->dataOnEntities[MBENTITYSET]) {
1448 CHKERR data.resetFieldDependentData();
1449 }
1450
1451 auto get_data_for_nodes = [&]() {
1453 if (!ss)
1454 CHKERR getRowNodesIndices(*op_data[ss], field_id[ss]);
1455 else
1456 CHKERR getColNodesIndices(*op_data[ss], field_id[ss]);
1457 CHKERR getNodesFieldData(*op_data[ss], field_id[ss]);
1459 };
1460
1461 // get data on entities but not nodes
1462 auto get_data_for_entities = [&]() {
1464 CHKERR getEntityFieldData(*op_data[ss], field_id[ss], MBEDGE);
1465 if (!ss)
1466 CHKERR getEntityRowIndices(*op_data[ss], field_id[ss], MBEDGE);
1467 else
1468 CHKERR getEntityColIndices(*op_data[ss], field_id[ss], MBEDGE);
1470 };
1471
1472 auto get_data_for_meshset = [&]() {
1474 if (!ss) {
1475 CHKERR getNoFieldRowIndices(*op_data[ss], field_id[ss]);
1476 } else {
1477 CHKERR getNoFieldColIndices(*op_data[ss], field_id[ss]);
1478 }
1479 CHKERR getNoFieldFieldData(*op_data[ss], field_id[ss]);
1481 };
1482
1483 switch (space[ss]) {
1484 case H1:
1485 CHKERR get_data_for_nodes();
1486 case HCURL:
1487 case HDIV:
1488 CHKERR get_data_for_entities();
1489 break;
1490 case L2:
1491 switch (type) {
1492 case MBVERTEX:
1493 CHKERR get_data_for_nodes();
1494 break;
1495 default:
1496 CHKERR get_data_for_entities();
1497 }
1498 break;
1499 case NOFIELD:
1500 // if (!getNinTheLoop()) {
1501 // NOFIELD data are the same for each element, can be
1502 // retrieved only once
1503 CHKERR get_data_for_meshset();
1504 // }
1505 break;
1506 default:
1507 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1508 "not implemented for this space < %s >",
1509 FieldSpaceNames[space[ss]]);
1510 }
1512 };
1513
1514 // evalate operators with field data
1515 auto evaluate_op_for_fields = [&](auto &op) {
1517
1518 if (op.getOpType() & UDO::OPROW) {
1519 try {
1520 CHKERR op.opRhs(*op_data[0], false);
1521 }
1522 CATCH_OP_ERRORS(op);
1523 }
1524
1525 if (op.getOpType() & UDO::OPCOL) {
1526 try {
1527 CHKERR op.opRhs(*op_data[1], false);
1528 }
1529 CATCH_OP_ERRORS(op);
1530 }
1531
1532 if (op.getOpType() & UDO::OPROWCOL) {
1533 try {
1534 CHKERR op.opLhs(*op_data[0], *op_data[1]);
1535 }
1536 CATCH_OP_ERRORS(op);
1537 }
1539 };
1540
1541 // Collect bit ref level on all entities, fields and spaces
1543
1544 auto oit = opPtrVector.begin();
1545 auto hi_oit = opPtrVector.end();
1546
1547 // interate over all operators on element
1548 for (; oit != hi_oit; oit++) {
1549
1550 try {
1551
1552 CHKERR oit->setPtrFE(this);
1553
1554 if ((oit->opType & UDO::OPSPACE) == UDO::OPSPACE) {
1555
1556 // run operator for space but specific field
1557 CHKERR evaluate_op_space(*oit);
1558
1559 } else if (
1560
1561 (oit->opType & (UDO::OPROW | UDO::OPCOL | UDO::OPROWCOL)) ==
1562 oit->opType
1563
1564 ) {
1565
1566 field_name[0] = oit->rowFieldName;
1567 field_name[1] = oit->colFieldName;
1568
1569 for (size_t ss = 0; ss != 2; ss++) {
1570
1571#ifndef NDEBUG
1572 if (field_name[ss].empty()) {
1573 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1574 "Not set Field name in operator %d (0-row, 1-column) in "
1575 "operator %s",
1576 ss,
1577 (boost::typeindex::type_id_runtime(*oit).pretty_name())
1578 .c_str());
1579 }
1580#endif
1581
1582 field_struture[ss] = mField.get_field_structure(field_name[ss]);
1583 field_id[ss] = field_struture[ss]->getBitNumber();
1584 space[ss] = field_struture[ss]->getSpace();
1585 base[ss] = field_struture[ss]->getApproxBase();
1586 }
1587
1588 // not that if field name do not change between operators, entity field
1589 // data are nor rebuild
1590 for (size_t ss = 0; ss != 2; ss++) {
1591
1592 if (oit->getOpType() & types[ss] ||
1593 oit->getOpType() & UDO::OPROWCOL) {
1594 if (last_eval_field_id[ss] != field_id[ss]) {
1595 CHKERR set_op_entityties_data(ss, *oit);
1596 last_eval_field_id[ss] = field_id[ss];
1597 }
1598 }
1599 }
1600
1601 CHKERR swap_bases(*oit);
1602
1603 // run operator for given field or row, column or both
1604 CHKERR evaluate_op_for_fields(*oit);
1605
1606 CHKERR swap_bases(*oit);
1607
1608 } else {
1609
1610 for (int i = 0; i != 5; ++i)
1611 if (oit->opType & (1 << i))
1612 MOFEM_LOG("SELF", Sev::error) << UDO::OpTypeNames[i];
1613 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1614 "Impossible operator type");
1615 }
1616 }
1617 CATCH_OP_ERRORS(*oit);
1618 }
1620}
1621
1623 "OPROW", " OPCOL", "OPROWCOL", "OPSPACE", "OPLAST"};
1624
1626 const std::string field_name, const EntityType type, const int side,
1627 VectorInt &indices) const {
1629 if (ptrFE == NULL)
1630 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1631
1632 switch (type) {
1633 case MBVERTEX:
1635 break;
1636 default:
1638 }
1640}
1641
1643 const std::string field_name, const EntityType type, const int side,
1644 VectorInt &indices) const {
1646 if (ptrFE == NULL)
1647 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1648
1649 switch (type) {
1650 case MBVERTEX:
1651 CHKERR ptrFE->getProblemNodesColIndices(field_name, indices);
1652 break;
1653 default:
1654 CHKERR ptrFE->getProblemTypeColIndices(field_name, type, side, indices);
1655 }
1657}
1658
1662 sidePtrFE = const_cast<ForcesAndSourcesCore *>(side_fe_ptr);
1664}
1665
1667 const ForcesAndSourcesCore *refine_fe_ptr) {
1669 refinePtrFE = const_cast<ForcesAndSourcesCore *>(refine_fe_ptr);
1671}
1672
1674 const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t side_dim,
1675 const EntityHandle ent_for_side, const int verb,
1676 const LogManager::SeverityLevel sev) {
1678 const auto ent = ent_for_side ? ent_for_side : getFEEntityHandle();
1679 const auto *problem_ptr = getFEMethod()->problemPtr;
1680
1681 side_fe->feName = fe_name;
1682
1683 CHKERR side_fe->setSideFEPtr(ptrFE);
1684 CHKERR side_fe->copyBasicMethod(*getFEMethod());
1685 CHKERR side_fe->copyPetscData(*getFEMethod());
1686 CHKERR side_fe->copyKsp(*getFEMethod());
1687 CHKERR side_fe->copySnes(*getFEMethod());
1688 CHKERR side_fe->copyTs(*getFEMethod());
1689
1690 CHKERR side_fe->preProcess();
1691
1692 Range adjacent_ents;
1693 CHKERR ptrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
1694 ent, side_dim, adjacent_ents);
1695 auto &numered_fe = problem_ptr->numeredFiniteElementsPtr
1697
1698 int nn = 0;
1699 side_fe->loopSize = adjacent_ents.size();
1700 for (auto fe_ent : adjacent_ents) {
1701 auto miit = numered_fe.find(boost::make_tuple(fe_name, fe_ent));
1702 if (miit != numered_fe.end()) {
1703
1704 if (verb >= VERBOSE)
1705 MOFEM_LOG("SELF", sev) << "Side finite element "
1706 << "(" << nn << "): " << **miit;
1707
1708 side_fe->nInTheLoop = nn++;
1709 side_fe->numeredEntFiniteElementPtr = *miit;
1710 CHKERR (*side_fe)();
1711 }
1712 }
1713
1714 CHKERR side_fe->postProcess();
1715
1717}
1718
1720 const string &fe_name, ForcesAndSourcesCore *this_fe, const int verb,
1721 const LogManager::SeverityLevel sev) {
1723
1724 if (verb >= VERBOSE)
1725 MOFEM_LOG("SELF", sev) << "This finite element: "
1726 << *getNumeredEntFiniteElementPtr();
1727
1728 const auto *problem_ptr = getFEMethod()->problemPtr;
1729 this_fe->feName = fe_name;
1730
1731 CHKERR this_fe->setRefineFEPtr(ptrFE);
1732 CHKERR this_fe->copyBasicMethod(*getFEMethod());
1733 CHKERR this_fe->copyPetscData(*getFEMethod());
1734 CHKERR this_fe->copyKsp(*getFEMethod());
1735 CHKERR this_fe->copySnes(*getFEMethod());
1736 CHKERR this_fe->copyTs(*getFEMethod());
1737
1738 CHKERR this_fe->preProcess();
1739
1740 this_fe->nInTheLoop = getNinTheLoop();
1741 this_fe->numeredEntFiniteElementPtr = getNumeredEntFiniteElementPtr();
1742
1743 CHKERR (*this_fe)();
1744
1745 CHKERR this_fe->postProcess();
1746
1748}
1749
1751 const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb,
1752 const LogManager::SeverityLevel sev) {
1754
1755 const auto *problem_ptr = getFEMethod()->problemPtr;
1756 auto &numered_fe = problem_ptr->numeredFiniteElementsPtr
1758
1759 const auto parent_ent = getNumeredEntFiniteElementPtr()->getParentEnt();
1760 auto miit = numered_fe.find(boost::make_tuple(fe_name, parent_ent));
1761 if (miit != numered_fe.end()) {
1762
1763 if (verb >= VERBOSE)
1764 MOFEM_LOG("SELF", sev) << "Parent finite element: " << **miit;
1765
1766 parent_fe->feName = fe_name;
1767
1768 CHKERR parent_fe->setRefineFEPtr(ptrFE);
1769 CHKERR parent_fe->copyBasicMethod(*getFEMethod());
1770 CHKERR parent_fe->copyPetscData(*getFEMethod());
1771 CHKERR parent_fe->copyKsp(*getFEMethod());
1772 CHKERR parent_fe->copySnes(*getFEMethod());
1773 CHKERR parent_fe->copyTs(*getFEMethod());
1774
1775 CHKERR parent_fe->preProcess();
1776
1777 parent_fe->loopSize = 1;
1778 parent_fe->nInTheLoop = 0;
1779 parent_fe->numeredEntFiniteElementPtr = *miit;
1780 CHKERR (*parent_fe)();
1781
1782 CHKERR parent_fe->postProcess();
1783 }
1784
1786}
1787
1789 const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb,
1790 const LogManager::SeverityLevel sev) {
1792
1793 const auto *problem_ptr = getFEMethod()->problemPtr;
1794 auto &ref_ents = *getPtrFE()->mField.get_ref_ents();
1795 auto &numered_fe = problem_ptr->numeredFiniteElementsPtr
1797
1798 const auto parent_ent = getNumeredEntFiniteElementPtr()->getEnt();
1799 const auto parent_type = getNumeredEntFiniteElementPtr()->getEntType();
1800 auto range =
1801 ref_ents.get<Composite_ParentEnt_And_EntType_mi_tag>().equal_range(
1802 boost::make_tuple(parent_type, parent_ent));
1803
1804 if (auto size = std::distance(range.first, range.second)) {
1805
1806 std::vector<EntityHandle> childs_vec;
1807 childs_vec.reserve(size);
1808 for (; range.first != range.second; ++range.first)
1809 childs_vec.emplace_back((*range.first)->getEnt());
1810
1811 Range childs;
1812
1813 if ((childs_vec.back() - childs_vec.front() + 1) == size)
1814 childs = Range(childs_vec.front(), childs_vec.back());
1815 else
1816 childs.insert_list(childs_vec.begin(), childs_vec.end());
1817
1818 child_fe->feName = fe_name;
1819
1820 CHKERR child_fe->setRefineFEPtr(ptrFE);
1821 CHKERR child_fe->copyBasicMethod(*getFEMethod());
1822 CHKERR child_fe->copyPetscData(*getFEMethod());
1823 CHKERR child_fe->copyKsp(*getFEMethod());
1824 CHKERR child_fe->copySnes(*getFEMethod());
1825 CHKERR child_fe->copyTs(*getFEMethod());
1826
1827 CHKERR child_fe->preProcess();
1828
1829 int nn = 0;
1830 child_fe->loopSize = size;
1831
1832 for (auto p = childs.pair_begin(); p != childs.pair_end(); ++p) {
1833
1834 auto miit = numered_fe.lower_bound(boost::make_tuple(fe_name, p->first));
1835 auto hi_miit =
1836 numered_fe.upper_bound(boost::make_tuple(fe_name, p->second));
1837
1838 for (; miit != hi_miit; ++miit) {
1839
1840 if (verb >= VERBOSE)
1841 MOFEM_LOG("SELF", sev) << "Child finite element "
1842 << "(" << nn << "): " << **miit;
1843
1844 child_fe->nInTheLoop = nn++;
1845 child_fe->numeredEntFiniteElementPtr = *miit;
1846 CHKERR (*child_fe)();
1847 }
1848 }
1849
1850 CHKERR child_fe->postProcess();
1851 };
1852
1854}
1855
1856int ForcesAndSourcesCore::getRule(int order_row, int order_col,
1857 int order_data) {
1858 return getRuleHook ? getRuleHook(order_row, order_col, order_data)
1859 : getRule(order_data);
1860}
1861
1863 int order_data) {
1864 return setRuleHook ? setRuleHook(this, order_row, order_col, order_data)
1865 : setGaussPts(order_data);
1866}
1867
1869
1870/** \deprecated setGaussPts(int row_order, int col_order, int data order);
1871 */
1874 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Sorry, not implemented");
1876}
1877
1879 const char type,
1880 const bool symm)
1881 : DataOperator(symm), opType(type), sPace(space), ptrFE(nullptr) {}
1882
1884 const std::string &field_name, const char type, const bool symm)
1885 : DataOperator(symm), opType(type), rowFieldName(field_name),
1886 colFieldName(field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1887
1889 const std::string &row_field_name, const std::string &col_field_name,
1890 const char type, const bool symm)
1891 : DataOperator(symm), opType(type), rowFieldName(row_field_name),
1892 colFieldName(col_field_name), sPace(LASTSPACE), ptrFE(nullptr) {}
1893
1896 if (preProcessHook) {
1897 ierr = preProcessHook();
1898 CHKERRG(ierr);
1899 }
1901}
1904 if (operatorHook) {
1905 ierr = operatorHook();
1906 CHKERRG(ierr);
1907 } else {
1908#ifndef NDEBUG
1909 MOFEM_LOG("SELF", Sev::warning)
1910 << "No method operator() overloaded on element entity on finite "
1911 "element <"
1912 << boost::typeindex::type_id_runtime(*this).pretty_name() << ">";
1913#endif
1914 }
1916}
1919 if (postProcessHook) {
1921 CHKERRG(ierr);
1922 }
1924}
1925
1929 if (!(ptrFE = dynamic_cast<ForcesAndSourcesCore *>(ptr)))
1930 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1931 "User operator and finite element do not work together");
1933}
1934
1935} // namespace MoFEM
static Index< 'n', 3 > n
static Index< 'p', 3 > p
#define CATCH_OP_ERRORS(OP)
constexpr double a
@ VERBOSE
Definition: definitions.h:222
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ LASTBASE
Definition: definitions.h:82
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ USER_BASE
user implemented approximation base
Definition: definitions.h:81
@ NOBASE
Definition: definitions.h:72
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
Definition: definitions.h:102
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:97
@ H1
continuous field
Definition: definitions.h:98
@ NOSPACE
Definition: definitions.h:96
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
static const char *const FieldSpaceNames[]
Definition: definitions.h:105
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
static const char *const ApproximationBaseNames[]
Definition: definitions.h:85
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
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.
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, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
virtual Field * get_field_structure(const std::string &name)=0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
SeverityLevel
Severity levels.
Definition: LogManager.hpp:43
FTensor::Index< 'i', SPACE_DIM > i
BlockParamData * cache
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
Definition: Types.hpp:54
int ApproximationOrder
Approximation on the entity.
Definition: Types.hpp:37
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
uint128_t UId
Unique Id.
Definition: Types.hpp:42
UBlasVector< int > VectorInt
Definition: Types.hpp:78
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static auto cmp_uid_lo(const boost::weak_ptr< FieldEntity > &a, const UId &b)
EntityHandle get_id_for_max_type()
static auto cmp_uid_hi(const UId &b, const boost::weak_ptr< FieldEntity > &a)
EntityHandle get_id_for_min_type()
ublas::vector< FieldEntity *, FieldEntAllocator > VectorFieldEntities
std::vector< boost::weak_ptr< FieldEntity > > FieldEntity_vector_view
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
static int getMaxOrder(const ENTMULTIINDEX &multi_index)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
const double r
rate factor
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr auto field_name
int loopSize
local number oe methods to process
MoFEMErrorCode copyBasicMethod(const BasicMethod &basic)
Copy data from other base method to this base method.
boost::function< MoFEMErrorCode()> preProcessHook
Hook function for pre-processing.
const Field_multiIndex * fieldsPtr
raw pointer to fields container
int nInTheLoop
number currently of processed method
const Problem * problemPtr
raw pointer to problem
boost::function< MoFEMErrorCode()> operatorHook
Hook function for operator.
int getNinTheLoop() const
get number of evaluated element in the loop
boost::function< MoFEMErrorCode()> postProcessHook
Hook function for post-processing.
Managing BitRefLevels.
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual BitFieldId get_field_id(const std::string &name) const =0
Get field Id.
virtual MPI_Comm & get_comm() const =0
base operator to do operations at Gauss Pt. level
static bool getDefTypeMap(const EntityType fe_type, const EntityType ent_type)
Deprecated interface functions.
this class derive data form other data structure
Data on single entity (This is passed as argument to DataOperator::doWork)
data structure for finite element entity
std::bitset< LASTSPACE > sPace
spaces on element
MatrixInt facesNodes
nodes on finite element faces
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::bitset< LASTBASE > bAse
bases on element
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
std::array< std::bitset< LASTBASE >, LASTSPACE > basesOnSpaces
base on spaces
MatrixInt facesNodesOrder
order of face nodes on element
std::array< std::bitset< LASTBASE >, MBMAXTYPE > basesOnEntities
bases on entity types
keeps information about indexed dofs for the finite element
std::string feName
Name of finite element.
auto & getRowFieldEnts() const
auto getFEName() const
get finite element name
const FieldEntity_vector_view & getDataFieldEnts() const
boost::shared_ptr< FieldEntity_vector_view > & getDataFieldEntsPtr() const
auto getColDofsPtr() const
EntityHandle getFEEntityHandle() const
auto getNumberOfNodes() const
auto getRowDofsPtr() const
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
auto & getColFieldEnts() const
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static UId getLoBitNumberUId(const FieldBitNumber bit_number)
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Provide data structure for (tensor) field approximation.
FieldCoefficientsNumber getNbOfCoeffs() const
Get number of field coefficients.
FieldSpace getSpace() const
Get field approximation space.
FieldBitNumber getBitNumber() const
Get number of set bit in Field ID. Each field has uid, get getBitNumber get number of bit set for giv...
MoFEMErrorCode getProblemColIndices(const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
Get col indices.
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
MoFEMErrorCode loopThis(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over parent elements. This function calls finite element with is oper...
MoFEMErrorCode loopParent(const string &fe_name, ForcesAndSourcesCore *parent_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over parent elements. This function calls finite element with is oper...
MoFEMErrorCode loopSide(const string &fe_name, ForcesAndSourcesCore *side_fe, const size_t dim, const EntityHandle ent_for_side=0, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over elements on the side of face. This function calls finite element...
MoFEMErrorCode getProblemRowIndices(const std::string filed_name, const EntityType type, const int side, VectorInt &indices) const
Get row indices.
MoFEMErrorCode loopChildren(const string &fe_name, ForcesAndSourcesCore *child_fe, const int verb=QUIET, const LogManager::SeverityLevel sev=Sev::noisy)
User call this function to loop over parent elements. This function calls finite element with is oper...
virtual MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
MoFEMErrorCode getProblemNodesRowIndices(const std::string &field_name, VectorInt &nodes_indices) const
MoFEMErrorCode getEntityIndices(EntitiesFieldData &data, const int bit_number, FieldEntity_vector_view &ents_field, const EntityType type_lo, const EntityType type_hi, EXTRACTOR &&extractor) const
MoFEMErrorCode getEntityRowIndices(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getRowNodesIndices(EntitiesFieldData &data, const int bit_number) const
get row node indices from FENumeredDofEntity_multiIndex
int getMaxRowOrder() const
Get max order of approximation for field in rows.
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
MoFEMErrorCode getProblemTypeColIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
MoFEMErrorCode getNoFieldIndices(const int bit_number, boost::shared_ptr< FENumeredDofEntity_multiIndex > dofs, VectorInt &nodes_indices) const
get NoField indices
virtual MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > derivedDataOnElement
Entity data on element entity columns fields.
MoFEMErrorCode getNoFieldRowIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
MoFEMErrorCode getNodesIndices(const int bit_number, FieldEntity_vector_view &ents_field, VectorInt &nodes_indices, VectorInt &local_nodes_indices, EXTRACTOR &&extractor) const
get node indices
MoFEMErrorCode getEntityFieldData(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode setRefineFEPtr(const ForcesAndSourcesCore *refine_fe_ptr)
Set the pointer to face element refined.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
ForcesAndSourcesCore(Interface &m_field)
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
MoFEMErrorCode getProblemNodesColIndices(const std::string &field_name, VectorInt &nodes_indices) const
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode getColNodesIndices(EntitiesFieldData &data, const int bit_number) const
get col node indices from FENumeredDofEntity_multiIndex
MoFEMErrorCode setSideFEPtr(const ForcesAndSourcesCore *side_fe_ptr)
Set the pointer to face element on the side.
int getMaxColOrder() const
Get max order of approximation for field in columns.
boost::ptr_vector< UserDataOperator > opPtrVector
Vector of finite element users data operators.
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
GaussHookFun setRuleHook
Set function to calculate integration rule.
MoFEMErrorCode getBitRefLevelOnData()
MoFEMErrorCode getNodesFieldData(EntitiesFieldData &data, const int bit_number) const
Get data on nodes.
MoFEMErrorCode getEntityColIndices(EntitiesFieldData &data, const int bit_number, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
ForcesAndSourcesCore * refinePtrFE
Element to integrate parent or child.
MoFEMErrorCode getNoFieldColIndices(EntitiesFieldData &data, const int bit_number) const
get col NoField indices
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode getNoFieldFieldData(const int bit_number, VectorDouble &ent_field_data, VectorDofs &ent_field_dofs, VectorFieldEntities &ent_field) const
Get field data on nodes.
MoFEMErrorCode getProblemTypeRowIndices(const std::string &field_name, EntityType type, int side_number, VectorInt &indices) const
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
RuleHookFun getRuleHook
Hook to get rule.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode copyKsp(const KspMethod &ksp)
copy data form another method
Definition: LoopMethods.cpp:68
MoFEMErrorCode copyPetscData(const PetscData &petsc_data)
Definition: LoopMethods.cpp:44
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredColDofsPtr
store DOFs on columns for this problem
MoFEMErrorCode copySnes(const SnesMethod &snes)
Copy snes data.
Definition: LoopMethods.cpp:88
MoFEMErrorCode copyTs(const TSMethod &ts)
Copy TS solver data.
boost::shared_ptr< SideNumber > getSideNumberPtr() const
Get the Side number.