19 if (
auto a_ptr =
a.lock()) {
20 if (a_ptr->getLocalUniqueId() < b)
30 if (
auto a_ptr =
a.lock()) {
31 if (b < a_ptr->getLocalUniqueId())
43 mField(m_field), getRuleHook(0), setRuleHook(0),
61 dataOnElement[
HCURL]),
67 dataNoField(*dataOnElement[
NOFIELD].get()),
68 dataH1(*dataOnElement[
H1].get()), dataHcurl(*dataOnElement[
HCURL].get()),
69 dataHdiv(*dataOnElement[
HDIV].get()), dataL2(*dataOnElement[
L2].get()),
70 lastEvaluatedElementEntityType(MBMAXTYPE), sidePtrFE(nullptr),
71 refinePtrFE(nullptr) {
86 boost::ptr_vector<EntitiesFieldData::EntData> &data)
const {
91 if (sit != side_table.end()) {
93 for (; sit != hi_sit; ++sit) {
94 if (
const auto side_number = (*sit)->side_number; side_number >= 0) {
95 const int brother_side_number = (*sit)->brother_side_number;
96 const int sense = (*sit)->sense;
98 data[side_number].getSense() = sense;
99 if (brother_side_number != -1)
100 data[brother_side_number].getSense() = sense;
109template <
typename ENTMULTIINDEX>
112 for (
auto ent_field_weak_ptr : multi_index)
113 if (
auto e = ent_field_weak_ptr.lock()) {
114 const int order = e->getMaxOrder();
115 max_order = (max_order <
order) ?
order : max_order;
123 if (
auto ptr = e.lock()) {
124 const int order = ptr->getMaxOrder();
125 max_order = (max_order <
order) ?
order : max_order;
141 boost::ptr_vector<EntitiesFieldData::EntData> &data)
const {
144 auto set_order = [&]() {
148 for (
unsigned int s = 0; s != data.size(); ++s)
149 data[s].getOrder() = 0;
154 r.first != r.second; ++r.first) {
156 const auto field_bit_number = (*r.first)->getBitNumber();
159 auto lo = std::lower_bound(data_field_ent.begin(), data_field_ent.end(),
161 if (lo != data_field_ent.end()) {
165 std::upper_bound(lo, data_field_ent.end(), hi_uid,
cmp_uid_hi);
166 for (; lo != hi; ++lo) {
168 if (
auto ptr = lo->lock()) {
171 auto sit = side_table.find(e.getEnt());
172 if (sit != side_table.end()) {
174 if (
const auto side_number = side->side_number;
177 auto &dat = data[side_number];
179 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
183 "Entity on side of the element not found");
192 auto set_order_on_brother = [&]() {
197 if (sit != side_table.end()) {
199 for (; sit != hi_sit; ++sit) {
200 const int brother_side_number = (*sit)->brother_side_number;
201 if (brother_side_number != -1) {
202 const int side_number = (*sit)->side_number;
203 data[brother_side_number].getOrder() = data[side_number].getOrder();
211 CHKERR set_order_on_brother();
218template <
typename EXTRACTOR>
222 EXTRACTOR &&extractor)
const {
230 if ((*field_it)->getBitNumber() != bit_number)
234 bit_number, get_id_for_min_type<MBVERTEX>());
235 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
237 if (lo != ents_field.end()) {
239 bit_number, get_id_for_max_type<MBVERTEX>());
240 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid,
cmp_uid_hi);
243 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
244 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
247 for (
auto it = lo; it != hi; ++it) {
248 if (
auto e = it->lock()) {
249 if (
auto cache = extractor(e).lock()) {
251 nb_dofs += std::distance(
cache->loHi[0],
cache->loHi[1]);
259 nodes_indices.resize(max_nb_dofs,
false);
260 local_nodes_indices.resize(max_nb_dofs,
false);
262 nodes_indices.resize(0,
false);
263 local_nodes_indices.resize(0,
false);
266 if (nb_dofs != max_nb_dofs) {
267 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
268 std::fill(local_nodes_indices.begin(), local_nodes_indices.end(), -1);
271 for (
auto it = lo; it != hi; ++it) {
272 if (
auto e = it->lock()) {
273 auto side_ptr = e->getSideNumberPtr();
274 if (
const auto side_number = side_ptr->side_number;
276 const auto brother_side_number = side_ptr->brother_side_number;
277 if (
auto cache = extractor(e).lock()) {
278 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1]; ++dit) {
280 const int idx = dof.getPetscGlobalDofIdx();
281 const int local_idx = dof.getPetscLocalDofIdx();
283 side_number * nb_dofs_on_vert + dof.getDofCoeffIdx();
284 nodes_indices[pos] = idx;
285 local_nodes_indices[pos] = local_idx;
286 if (brother_side_number != -1) {
287 const int elem_idx = brother_side_number * nb_dofs_on_vert +
288 (*dit)->getDofCoeffIdx();
289 nodes_indices[elem_idx] = idx;
290 local_nodes_indices[elem_idx] = local_idx;
298 nodes_indices.resize(0,
false);
299 local_nodes_indices.resize(0,
false);
303 nodes_indices.resize(0,
false);
304 local_nodes_indices.resize(0,
false);
312 const int bit_number)
const {
315 boost::weak_ptr<EntityCacheNumeredDofs>
316 operator()(boost::shared_ptr<FieldEntity> &e) {
317 return e->entityCacheRowDofs;
329 const int bit_number)
const {
332 boost::weak_ptr<EntityCacheNumeredDofs>
333 operator()(boost::shared_ptr<FieldEntity> &e) {
334 return e->entityCacheColDofs;
344template <
typename EXTRACTOR>
348 const EntityType type_hi, EXTRACTOR &&extractor)
const {
353 dat.getIndices().resize(0,
false);
354 dat.getLocalIndices().resize(0,
false);
360 auto lo = std::lower_bound(ents_field.begin(), ents_field.end(), lo_uid,
362 if (lo != ents_field.end()) {
365 auto hi = std::upper_bound(lo, ents_field.end(), hi_uid,
cmp_uid_hi);
368 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
370 for (
auto it = lo; it != hi; ++it)
371 if (
auto e = it->lock()) {
374 auto side_ptr = e->getSideNumberPtr();
375 if (
const auto side = side_ptr->side_number; side >= 0) {
376 const auto nb_dofs_on_ent = e->getNbDofsOnEnt();
377 const auto brother_side = side_ptr->brother_side_number;
379 auto &ent_field_indices = dat.getIndices();
380 auto &ent_field_local_indices = dat.getLocalIndices();
382 ent_field_indices.resize(nb_dofs_on_ent,
false);
383 ent_field_local_indices.resize(nb_dofs_on_ent,
false);
384 std::fill(ent_field_indices.data().begin(),
385 ent_field_indices.data().end(), -1);
386 std::fill(ent_field_local_indices.data().begin(),
387 ent_field_local_indices.data().end(), -1);
389 if (
auto cache = extractor(e).lock()) {
390 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1]; ++dit) {
391 const int idx = (*dit)->getEntDofIdx();
392 ent_field_indices[idx] = (*dit)->getPetscGlobalDofIdx();
393 ent_field_local_indices[idx] = (*dit)->getPetscLocalDofIdx();
397 if (brother_side != -1) {
399 dat_brother.getIndices().resize(nb_dofs_on_ent,
false);
400 dat_brother.getLocalIndices().resize(nb_dofs_on_ent,
false);
401 noalias(dat_brother.getIndices()) = dat.getIndices();
402 noalias(dat_brother.getLocalIndices()) = dat.getLocalIndices();
417 boost::weak_ptr<EntityCacheNumeredDofs>
418 operator()(boost::shared_ptr<FieldEntity> &e) {
419 return e->entityCacheRowDofs;
432 boost::weak_ptr<EntityCacheNumeredDofs>
433 operator()(boost::shared_ptr<FieldEntity> &e) {
434 return e->entityCacheColDofs;
443 const int bit_number, boost::shared_ptr<FENumeredDofEntity_multiIndex> dofs,
450 indices.resize(std::distance(dit, hi_dit));
451 for (; dit != hi_dit; dit++) {
452 int idx = (*dit)->getPetscGlobalDofIdx();
453 indices[(*dit)->getDofCoeffIdx()] = idx;
460 const int bit_number)
const {
465 "data.dataOnEntities[MBENTITYSET] is empty");
475 const int bit_number)
const {
480 "data.dataOnEntities[MBENTITYSET] is empty");
499 nodes_indices.resize(field_struture->
getNbOfCoeffs() * num_nodes,
false);
500 std::fill(nodes_indices.begin(), nodes_indices.end(), -1);
504 auto siit = side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 0));
506 side_table.get<1>().lower_bound(boost::make_tuple(MBVERTEX, 10000));
509 for (; siit != hi_siit; siit++, nn++) {
510 if (siit->get()->side_number >= 0) {
517 for (; dit != hi_dit; dit++) {
518 nodes_indices[siit->get()->side_number * (*dit)->getNbOfCoeffs() +
519 (*dit)->getDofCoeffIdx()] =
520 (*dit)->getPetscGlobalDofIdx();
525 nodes_indices.resize(0,
false);
541 side_table.get<1>().lower_bound(boost::make_tuple(type, side_number));
543 side_table.get<1>().upper_bound(boost::make_tuple(type, side_number));
545 for (; siit != hi_siit; siit++) {
546 if (siit->get()->side_number >= 0) {
554 indices.resize(std::distance(dit, hi_dit));
555 for (; dit != hi_dit; dit++) {
557 indices[(*dit)->getEntDofIdx()] = (*dit)->getPetscGlobalDofIdx();
576 type, side_number, indices);
590 type, side_number, indices);
604 for (
auto &dat : data->dataOnEntities) {
605 for (
auto &ent_dat : dat) {
606 ent_dat.getEntDataBitRefLevel().reset();
613 for (
auto it : field_ents) {
614 if (
auto e = it.lock()) {
618 const signed char side =
619 type == MBVERTEX ? 0 : e->getSideNumberPtr()->side_number;
622 auto &dat = data->dataOnEntities[type][side];
623 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
628 auto &dat = data->dataOnEntities[MBENTITYSET][0];
629 dat.getEntDataBitRefLevel() |= e->getBitRefLevel();
640 const int bit_number)
const {
642 auto get_nodes_field_data = [&](
VectorDouble &nodes_data,
649 nodes_data.resize(0,
false);
650 nodes_dofs.resize(0,
false);
651 field_entities.resize(0,
false);
658 if ((*field_it)->getBitNumber() != bit_number)
661 const int nb_dofs_on_vert = (*field_it)->getNbOfCoeffs();
662 space = (*field_it)->getSpace();
663 base = (*field_it)->getApproxBase();
667 bit_number, get_id_for_min_type<MBVERTEX>());
668 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
670 if (lo != field_ents.end()) {
672 bit_number, get_id_for_max_type<MBVERTEX>());
673 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
677 for (
auto it = lo; it != hi; ++it) {
678 if (
auto e = it->lock()) {
679 nb_dofs += e->getNbDofsOnEnt();
686 bb_node_order.resize(num_nodes,
false);
687 bb_node_order.clear();
688 const int max_nb_dofs = nb_dofs_on_vert * num_nodes;
689 nodes_data.resize(max_nb_dofs,
false);
690 nodes_dofs.resize(max_nb_dofs,
false);
691 field_entities.resize(num_nodes,
false);
692 std::fill(nodes_data.begin(), nodes_data.end(), 0);
693 std::fill(nodes_dofs.begin(), nodes_dofs.end(),
nullptr);
694 std::fill(field_entities.begin(), field_entities.end(),
nullptr);
696 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
698 for (
auto it = lo; it != hi; ++it) {
699 if (
auto e = it->lock()) {
700 const auto &sn = e->getSideNumberPtr();
703 if (
const auto side_number = sn->side_number;
705 const int brother_side_number = sn->brother_side_number;
707 field_entities[side_number] = e.get();
708 if (brother_side_number != -1) {
709 brother_ents_vec.emplace_back(e);
710 field_entities[side_number] = field_entities[side_number];
713 bb_node_order[side_number] = e->getMaxOrder();
714 int pos = side_number * nb_dofs_on_vert;
715 auto ent_filed_data_vec = e->getEntFieldData();
716 if (
auto cache = e->entityCacheDataDofs.lock()) {
717 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1];
719 const auto dof_idx = (*dit)->getEntDofIdx();
720 nodes_data[pos + dof_idx] = ent_filed_data_vec[dof_idx];
721 nodes_dofs[pos + dof_idx] =
729 for (
auto &it : brother_ents_vec) {
730 if (
const auto e = it.lock()) {
732 const int side_number = sn->side_number;
733 const int brother_side_number = sn->brother_side_number;
734 bb_node_order[brother_side_number] = bb_node_order[side_number];
735 int pos = side_number * nb_dofs_on_vert;
736 int brother_pos = brother_side_number * nb_dofs_on_vert;
737 for (
int ii = 0; ii != nb_dofs_on_vert; ++ii) {
738 nodes_data[brother_pos] = nodes_data[pos];
739 nodes_dofs[brother_pos] = nodes_dofs[pos];
753 return get_nodes_field_data(
771 dat.getFieldData().resize(0,
false);
772 dat.getFieldDofs().resize(0,
false);
773 dat.getFieldEntities().resize(0,
false);
780 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
782 if (lo != field_ents.end()) {
785 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
788 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
790 for (
auto it = lo; it != hi; ++it)
791 if (
auto e = it->lock()) {
792 auto side_ptr = e->getSideNumberPtr();
793 if (
const auto side = side_ptr->side_number; side >= 0) {
796 auto &ent_field = dat.getFieldEntities();
797 auto &ent_field_dofs = dat.getFieldDofs();
798 auto &ent_field_data = dat.getFieldData();
800 const int brother_side = side_ptr->brother_side_number;
801 if (brother_side != -1)
802 brother_ents_vec.emplace_back(e);
804 dat.getBase() = e->getApproxBase();
805 dat.getSpace() = e->getSpace();
806 const int ent_order = e->getMaxOrder();
808 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
810 auto ent_data = e->getEntFieldData();
811 ent_field_data.resize(ent_data.size(),
false);
812 noalias(ent_field_data) = ent_data;
813 ent_field_dofs.resize(ent_data.size(),
false);
814 std::fill(ent_field_dofs.begin(), ent_field_dofs.end(),
nullptr);
815 ent_field.resize(1,
false);
816 ent_field[0] = e.get();
817 if (
auto cache = e->entityCacheDataDofs.lock()) {
818 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1]; ++dit) {
819 ent_field_dofs[(*dit)->getEntDofIdx()] =
826 for (
auto &it : brother_ents_vec) {
827 if (
const auto e = it.lock()) {
829 const int side = e->getSideNumberPtr()->side_number;
830 const int brother_side = e->getSideNumberPtr()->brother_side_number;
833 dat_brother.getBase() = dat.getBase();
834 dat_brother.getSpace() = dat.getSpace();
835 dat_brother.getOrder() = dat.getOrder();
836 dat_brother.getFieldData() = dat.getFieldData();
837 dat_brother.getFieldDofs() = dat.getFieldDofs();
838 dat_brother.getFieldEntities() = dat.getFieldEntities();
853 ent_field_data.resize(0,
false);
854 ent_field_dofs.resize(0,
false);
855 ent_field.resize(0,
false);
859 bit_number, get_id_for_min_type<MBVERTEX>());
860 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
862 if (lo != field_ents.end()) {
864 ent_field.resize(field_ents.size(),
false);
865 std::fill(ent_field.begin(), ent_field.end(),
nullptr);
868 bit_number, get_id_for_max_type<MBENTITYSET>());
869 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
873 for (
auto it = lo; it != hi; ++it, ++side)
874 if (
auto e = it->lock()) {
876 const auto size = e->getNbDofsOnEnt();
877 ent_field_data.resize(size,
false);
878 ent_field_dofs.resize(size,
false);
879 ent_field[side] = e.get();
880 noalias(ent_field_data) = e->getEntFieldData();
882 if (
auto cache = e->entityCacheDataDofs.lock()) {
883 for (
auto dit =
cache->loHi[0]; dit !=
cache->loHi[1]; ++dit) {
884 ent_field_dofs[(*dit)->getEntDofIdx()] =
897 const int bit_number)
const {
901 "No space to insert data");
921 const auto nb_faces = CN::NumSubEntities(type, 2);
925 auto side_ptr_it = side_table.get<1>().lower_bound(
926 boost::make_tuple(CN::TypeDimensionMap[2].first, 0));
927 auto hi_side_ptr_it = side_table.get<1>().upper_bound(
928 boost::make_tuple(CN::TypeDimensionMap[2].second, 100));
930 for (; side_ptr_it != hi_side_ptr_it; ++side_ptr_it) {
931 const auto side = (*side_ptr_it)->side_number;
932 const auto sense = (*side_ptr_it)->sense;
933 const auto offset = (*side_ptr_it)->offset;
934 const auto face_ent = (*side_ptr_it)->ent;
939 CN::SubEntityVertexIndices(type, 2, side, face_type, nb_nodes_face);
940 face_nodes.resize(nb_faces, nb_nodes_face);
941 face_nodes_order.resize(nb_faces, nb_nodes_face);
944 for (
int n = 0;
n != nb_nodes_face; ++
n)
945 face_nodes_order(side,
n) = (
n + offset) % nb_nodes_face;
947 for (
int n = 0;
n != nb_nodes_face; ++
n)
948 face_nodes_order(side,
n) =
949 (nb_nodes_face - (
n - offset) % nb_nodes_face) % nb_nodes_face;
951 for (
int n = 0;
n != nb_nodes_face; ++
n)
952 face_nodes(side,
n) = face_indices[face_nodes_order(side,
n)];
960 nb_nodes_face,
true);
961 face_nodes.resize(nb_faces, nb_nodes_face);
962 for (
int nn = 0; nn != nb_nodes_face; ++nn) {
963 if (face_nodes(side, nn) !=
966 std::find(conn_ele, &conn_ele[num_nodes_ele], conn_face[nn]))) {
968 "Wrong face numeration");
1002 if (
auto ptr = e.lock()) {
1009 data.
sPace.set(space);
1010 data.
bAse.set(approx);
1019 "data fields ents not allocated on element");
1039 "Functions genrating approximation base not defined");
1041 for (
int space =
H1; space !=
LASTSPACE; ++space) {
1047 boost::make_shared<EntPolynomialBaseCtx>(
1056 "Functions generating user approximation base not defined");
1064 boost::make_shared<EntPolynomialBaseCtx>(
1071 "Base <%s> not yet implemented",
1085 dataOnElement[space]->dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(
1107 auto &bb_node_order = data.
dataOnEntities[MBVERTEX][0].getBBNodeOrder();
1110 auto bit_number = field_ptr->getBitNumber();
1112 bit_number, get_id_for_min_type<MBVERTEX>());
1113 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1115 if (lo != field_ents.end()) {
1117 bit_number, get_id_for_max_type<MBVERTEX>());
1118 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
1121 for (
auto it = lo; it != hi; ++it)
1122 if (
auto first_e = it->lock()) {
1123 space = first_e->getSpace();
1124 base = first_e->getApproxBase();
1126 bb_node_order.resize(num_nodes,
false);
1127 bb_node_order.clear();
1128 const int nb_dof_idx = first_e->getNbOfCoeffs();
1130 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1132 for (; it != hi; ++it) {
1133 if (
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 if (brother_side_number != -1)
1138 brother_ents_vec.emplace_back(e);
1139 bb_node_order[side_number] = e->getMaxOrder();
1143 for (
auto &it : brother_ents_vec) {
1144 if (
const auto e = it.lock()) {
1145 const auto &sn = e->getSideNumberPtr();
1146 const int side_number = sn->side_number;
1147 const int brother_side_number = sn->brother_side_number;
1148 bb_node_order[brother_side_number] = bb_node_order[side_number];
1169 dat.getFieldData().resize(0,
false);
1170 dat.getFieldDofs().resize(0,
false);
1175 auto bit_number = field_ptr->getBitNumber();
1178 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(), lo_uid,
1180 if (lo != field_ents.end()) {
1183 auto hi = std::upper_bound(lo, field_ents.end(), hi_uid,
cmp_uid_hi);
1185 std::vector<boost::weak_ptr<FieldEntity>> brother_ents_vec;
1186 for (; lo != hi; ++lo) {
1187 if (
auto e = lo->lock()) {
1188 if (
auto cache = e->entityCacheDataDofs.lock()) {
1190 if (
const auto side = e->getSideNumberPtr()->side_number;
1194 const int brother_side =
1195 e->getSideNumberPtr()->brother_side_number;
1196 if (brother_side != -1)
1197 brother_ents_vec.emplace_back(e);
1198 dat.getBase() = e->getApproxBase();
1199 dat.getSpace() = e->getSpace();
1200 const auto ent_order = e->getMaxOrder();
1202 dat.getOrder() > ent_order ? dat.getOrder() : ent_order;
1209 for (
auto &ent_ptr : brother_ents_vec) {
1210 if (
auto e = ent_ptr.lock()) {
1212 const int side = e->getSideNumberPtr()->side_number;
1213 const int brother_side = e->getSideNumberPtr()->brother_side_number;
1216 dat_brother.getBase() = dat.getBase();
1217 dat_brother.getSpace() = dat.getSpace();
1218 dat_brother.getOrder() = dat.getOrder();
1226 if (
auto ent_data_ptr = e.lock()) {
1228 auto space = ent_data_ptr->getSpace();
1231 for (
auto &ptr : dat.getBBAlphaIndicesByOrderArray())
1233 for (
auto &ptr : dat.getBBNByOrderArray())
1235 for (
auto &ptr : dat.getBBDiffNByOrderArray())
1243 auto check_space = [&](
const auto space) {
1246 for (
auto t = MBVERTEX;
t <= ele_type; ++
t) {
1252 for (
auto t = MBEDGE;
t <= ele_type; ++
t) {
1258 for (
auto t = MBTRI;
t <= ele_type; ++
t) {
1271 std::set<string> fields_list;
1273 if (
auto ent_data_ptr = e.lock()) {
1276 if (fields_list.find(
field_name) == fields_list.end()) {
1277 auto field_ptr = ent_data_ptr->getFieldRawPtr();
1278 auto space = ent_data_ptr->getSpace();
1282 if (check_space(space)) {
1284 gaussPts, boost::make_shared<EntPolynomialBaseCtx>(
1302 for (
int space =
H1; space !=
LASTSPACE; ++space) {
1310#define FUNCTION_NAME_WITH_OP_NAME(OP) \
1311 std::ostringstream ss; \
1312 ss << "(Calling user data operator " \
1313 << boost::typeindex::type_id_runtime(OP).pretty_name() << " rowField " \
1314 << (OP).rowFieldName << " colField " << (OP).colFieldName << ") "
1316#define CATCH_OP_ERRORS(OP) \
1317 catch (MoFEMExceptionInitial const &ex) { \
1318 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1319 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1320 ex.errorCode, PETSC_ERROR_INITIAL, ex.what()); \
1322 catch (MoFEMExceptionRepeat const &ex) { \
1323 FUNCTION_NAME_WITH_OP_NAME(OP) << PETSC_FUNCTION_NAME; \
1324 return PetscError(PETSC_COMM_SELF, ex.lINE, ss.str().c_str(), __FILE__, \
1325 ex.errorCode, PETSC_ERROR_REPEAT, " "); \
1327 catch (MoFEMException const &ex) { \
1328 FUNCTION_NAME_WITH_OP_NAME(OP) << ex.errorMessage; \
1329 SETERRQ(PETSC_COMM_SELF, ex.errorCode, ss.str().c_str()); \
1331 catch (std::exception const &ex) { \
1332 std::string message("Error: " + std::string(ex.what()) + " at " + \
1333 boost::lexical_cast<std::string>(__LINE__) + " : " + \
1334 std::string(__FILE__) + " in " + \
1335 std::string(PETSC_FUNCTION_NAME)); \
1336 FUNCTION_NAME_WITH_OP_NAME(OP) << message; \
1337 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str()); \
1346 std::array<const Field *, 3> field_struture;
1349 std::array<FieldSpace, 2> space;
1350 std::array<FieldApproximationBase, 2> base;
1352 constexpr std::array<UDO::OpType, 2> types = {UDO::OPROW, UDO::OPCOL};
1353 std::array<int, 2> last_eval_field_id = {0, 0};
1355 std::array<boost::shared_ptr<EntitiesFieldData>, 2> op_data;
1357 auto swap_bases = [&](
auto &op) {
1359 for (
size_t ss = 0; ss != 2; ++ss) {
1360 if (op.getOpType() & types[ss] || op.getOpType() & UDO::OPROWCOL) {
1377 auto evaluate_op_space = [&](
auto &op) {
1382 std::fill(last_eval_field_id.begin(), last_eval_field_id.end(), 0);
1402 e.getSpace() = op.sPace;
1412 "Not implemented for this space", op.sPace);
1419 auto set_op_entityties_data = [&](
auto ss,
auto &op) {
1423 if ((op.getNumeredEntFiniteElementPtr()->getBitFieldIdData() &
1427 "no data field < %s > on finite element < %s >",
1435 for (
auto &data : op_data[ss]->dataOnEntities[MBENTITYSET]) {
1436 CHKERR data.resetFieldDependentData();
1439 auto get_data_for_nodes = [&]() {
1450 auto get_data_for_entities = [&]() {
1460 auto get_data_for_meshset = [&]() {
1471 switch (space[ss]) {
1473 CHKERR get_data_for_nodes();
1476 CHKERR get_data_for_entities();
1481 CHKERR get_data_for_nodes();
1484 CHKERR get_data_for_entities();
1491 CHKERR get_data_for_meshset();
1496 "not implemented for this space < %s >",
1503 auto evaluate_op_for_fields = [&](
auto &op) {
1506 if (op.getOpType() & UDO::OPROW) {
1508 CHKERR op.opRhs(*op_data[0],
false);
1513 if (op.getOpType() & UDO::OPCOL) {
1515 CHKERR op.opRhs(*op_data[1],
false);
1520 if (op.getOpType() & UDO::OPROWCOL) {
1522 CHKERR op.opLhs(*op_data[0], *op_data[1]);
1536 for (; oit != hi_oit; oit++) {
1540 CHKERR oit->setPtrFE(
this);
1542 if ((oit->opType & UDO::OPSPACE) == UDO::OPSPACE) {
1545 CHKERR evaluate_op_space(*oit);
1549 (oit->opType & (UDO::OPROW | UDO::OPCOL | UDO::OPROWCOL)) ==
1557 for (
size_t ss = 0; ss != 2; ss++) {
1562 "Not set Field name in operator %d (0-row, 1-column) in "
1565 (boost::typeindex::type_id_runtime(*oit).pretty_name())
1572 space[ss] = field_struture[ss]->getSpace();
1573 base[ss] = field_struture[ss]->getApproxBase();
1578 for (
size_t ss = 0; ss != 2; ss++) {
1580 if (oit->getOpType() & types[ss] ||
1581 oit->getOpType() & UDO::OPROWCOL) {
1582 if (last_eval_field_id[ss] != field_id[ss]) {
1583 CHKERR set_op_entityties_data(ss, *oit);
1584 last_eval_field_id[ss] = field_id[ss];
1592 CHKERR evaluate_op_for_fields(*oit);
1598 for (
int i = 0;
i != 5; ++
i)
1599 if (oit->opType & (1 <<
i))
1600 MOFEM_LOG(
"SELF", Sev::error) << UDO::OpTypeNames[
i];
1602 "Impossible operator type");
1610const char *
const ForcesAndSourcesCore::UserDataOperator::OpTypeNames[] = {
1611 "OPROW",
" OPCOL",
"OPROWCOL",
"OPSPACE",
"OPLAST"};
1667 const auto *problem_ptr = getFEMethod()->problemPtr;
1669 side_fe->
feName = fe_name;
1680 Range adjacent_ents;
1682 ent, side_dim, adjacent_ents);
1684 problem_ptr->numeredFiniteElementsPtr->get<
Unique_mi_tag>();
1686 auto fe_miit = ptrFE->mField.get_finite_elements()
1689 if (fe_miit != ptrFE->mField.get_finite_elements()
1693 side_fe->
loopSize = adjacent_ents.size();
1694 for (
auto fe_ent : adjacent_ents) {
1696 fe_ent, (*fe_miit)->getFEUId()));
1697 if (miit != numered_fe.end()) {
1699 MOFEM_LOG(
"SELF", sev) <<
"Side finite element "
1700 <<
"(" << nn <<
"): " << **miit;
1719 MOFEM_LOG(
"SELF", sev) <<
"This finite element: "
1720 << *getNumeredEntFiniteElementPtr();
1722 const auto *problem_ptr = getFEMethod()->problemPtr;
1723 this_fe->
feName = fe_name;
1749 const auto *problem_ptr = getFEMethod()->problemPtr;
1751 problem_ptr->numeredFiniteElementsPtr->get<
Unique_mi_tag>();
1753 parent_fe->
feName = fe_name;
1763 auto fe_miit = ptrFE->mField.get_finite_elements()
1766 if (fe_miit != ptrFE->mField.get_finite_elements()
1769 const auto parent_ent = getNumeredEntFiniteElementPtr()->getParentEnt();
1771 parent_ent, (*fe_miit)->getFEUId()));
1772 if (miit != numered_fe.end()) {
1774 MOFEM_LOG(
"SELF", sev) <<
"Parent finite element: " << **miit;
1792 const auto *problem_ptr = getFEMethod()->problemPtr;
1793 auto &ref_ents = *getPtrFE()->mField.get_ref_ents();
1795 problem_ptr->numeredFiniteElementsPtr->get<
Unique_mi_tag>();
1797 const auto parent_ent = getNumeredEntFiniteElementPtr()->getEnt();
1798 const auto parent_type = getNumeredEntFiniteElementPtr()->getEntType();
1801 boost::make_tuple(parent_type, parent_ent));
1803 if (
auto size = std::distance(range.first, range.second)) {
1805 std::vector<EntityHandle> childs_vec;
1806 childs_vec.reserve(size);
1807 for (; range.first != range.second; ++range.first)
1808 childs_vec.emplace_back((*range.first)->getEnt());
1812 if ((childs_vec.back() - childs_vec.front() + 1) == size)
1813 childs =
Range(childs_vec.front(), childs_vec.back());
1815 childs.insert_list(childs_vec.begin(), childs_vec.end());
1817 child_fe->
feName = fe_name;
1831 auto fe_miit = ptrFE->mField.get_finite_elements()
1834 if (fe_miit != ptrFE->mField.get_finite_elements()
1838 for (
auto p = childs.pair_begin();
p != childs.pair_end(); ++
p) {
1842 p->first, (*fe_miit)->getFEUId()));
1845 p->second, (*fe_miit)->getFEUId()));
1847 for (; miit != hi_miit; ++miit) {
1850 MOFEM_LOG(
"SELF", sev) <<
"Child finite element "
1851 <<
"(" << nn <<
"): " << **miit;
1888ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
const FieldSpace space,
1891 :
DataOperator(symm), opType(type), sPace(space), ptrFE(nullptr) {}
1893ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1894 const std::string
field_name,
const char type,
const bool symm)
1898ForcesAndSourcesCore::UserDataOperator::UserDataOperator(
1899 const std::string row_field_name,
const std::string col_field_name,
1900 const char type,
const bool symm)
1901 :
DataOperator(symm), opType(type), rowFieldName(row_field_name),
1902 colFieldName(col_field_name), sPace(
LASTSPACE), ptrFE(nullptr) {}
1920 <<
"No method operator() overloaded on element entity on finite "
1922 << boost::typeindex::type_id_runtime(*this).pretty_name() <<
">";
1941 "User operator and finite element do not work together");
#define CATCH_OP_ERRORS(OP)
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ USER_BASE
user implemented approximation base
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
@ LASTSPACE
FieldSpace in [ 0, LASTSPACE )
@ NOFIELD
scalar or vector of scalars describe (no true field)
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
static const char *const FieldSpaceNames[]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
static const char *const ApproximationBaseNames[]
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FTensor::Index< 'n', SPACE_DIM > n
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 const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
#define MOFEM_LOG(channel, severity)
Log.
SeverityLevel
Severity levels.
FTensor::Index< 'i', SPACE_DIM > i
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITFEID_SIZE > BitFEId
Finite element Id.
int ApproximationOrder
Approximation on the entity.
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
UBlasVector< int > VectorInt
implementation of Data Operators for Forces and Sources
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
constexpr double t
plate stiffness
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.
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
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.
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...
ForcesAndSourcesCore * ptrFE
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
boost::ptr_deque< UserDataOperator > opPtrVector
Vector of finite element users data operators.
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.
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
MoFEMErrorCode copyPetscData(const PetscData &petsc_data)
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.
MoFEMErrorCode copyTs(const TSMethod &ts)
Copy TS solver data.
boost::shared_ptr< SideNumber > getSideNumberPtr() const