v0.15.0
Loading...
Searching...
No Matches
FieldCore.cpp
Go to the documentation of this file.
1/** \file FieldCore.cpp
2 * \brief Core interface methods for managing fields.
3 */
4
5
6#include <MoFEM.hpp>
7
8#define FieldCoreFunctionBegin \
9 MoFEMFunctionBegin; \
10 MOFEM_LOG_CHANNEL("WORLD"); \
11 MOFEM_LOG_CHANNEL("SYNC"); \
12 MOFEM_LOG_FUNCTION(); \
13 MOFEM_LOG_TAG("WORLD", "FieldCore"); \
14 MOFEM_LOG_TAG("SYNC", "FieldCore");
15
16namespace MoFEM {
17
18BitFieldId Core::get_field_id(const std::string &name) const {
19 auto &set = fIelds.get<FieldName_mi_tag>();
20 auto miit = set.find(name);
21 if (miit == set.end()) {
22 THROW_MESSAGE("field < " + name +
23 " > not in database (top tip: check spelling)");
24 }
25 return (*miit)->getId();
26}
27
28FieldBitNumber Core::get_field_bit_number(const std::string name) const {
29 auto &set = fIelds.get<FieldName_mi_tag>();
30 auto miit = set.find(name);
31 if (miit == set.end())
32 THROW_MESSAGE("field < " + name +
33 " > not in database (top tip: check spelling)");
34 return (*miit)->getBitNumber();
35}
36
37std::string Core::get_field_name(const BitFieldId id) const {
38 auto &set = fIelds.get<BitFieldId_mi_tag>();
39 auto miit = set.find(id);
40 if (miit == set.end())
41 THROW_MESSAGE("field not in database (top tip: check spelling)");
42 return (*miit)->getName();
43}
44
46 auto &set = fIelds.get<BitFieldId_mi_tag>();
47 auto miit = set.find(id);
48 if (miit == set.end())
49 THROW_MESSAGE("field not in database (top tip: check spelling)");
50 return (*miit)->meshSet;
51}
52
53EntityHandle Core::get_field_meshset(const std::string name) const {
54 return get_field_meshset(get_field_id(name));
55}
56
57bool Core::check_field(const std::string &name) const {
58 auto miit = fIelds.get<FieldName_mi_tag>().find(name);
59 if (miit == fIelds.get<FieldName_mi_tag>().end())
60 return false;
61 else
62 return true;
63}
64
65const Field *Core::get_field_structure(const std::string &name,
66 enum MoFEMTypes bh) const {
67 auto miit = fIelds.get<FieldName_mi_tag>().find(name);
68 if (miit == fIelds.get<FieldName_mi_tag>().end()) {
69 if (bh == MF_EXIST)
70 throw MoFEMException(
72 std::string("field < " + name +
73 " > not in database (top tip: check spelling)")
74 .c_str());
75 else
76 return nullptr;
77 }
78 return miit->get();
79}
80
82 int dim,
83 Range &ents) const {
84
86 EntityHandle meshset = get_field_meshset(name);
87 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, true);
89}
90
92 EntityType type,
93 Range &ents) const {
95 EntityHandle meshset = get_field_meshset(name);
96 CHKERR get_moab().get_entities_by_type(meshset, type, ents, true);
98}
99
101 Range &ents) const {
103 EntityHandle meshset = get_field_meshset(name);
104 CHKERR get_moab().get_entities_by_handle(meshset, ents, true);
106}
107
108MoFEMErrorCode Core::addField(const std::string &name, const FieldSpace space,
109 const FieldContinuity continuity,
110 const FieldApproximationBase base,
111 const FieldCoefficientsNumber nb_of_coefficients,
112 const TagType tag_type, const enum MoFEMTypes bh,
113 int verb) {
114 MOFEM_LOG_CHANNEL("WORLD");
115 MOFEM_LOG_TAG("WORLD", "FieldCore");
116 MOFEM_LOG_CHANNEL("SYNC");
117 MOFEM_LOG_TAG("SYNC", "FieldCore");
120
121 // Add field mesh set to partion meshset. In case of no elements
122 // on processor part, when mesh file is read, finite element meshset is
123 // prevented from deletion by moab reader.
124 auto add_meshset_to_partition = [&](auto meshset) {
126 const void *tag_vals[] = {&rAnk};
127 ParallelComm *pcomm = ParallelComm::get_pcomm(
128 &get_moab(), get_basic_entity_data_ptr()->pcommID);
129 Tag part_tag = pcomm->part_tag();
130 Range tagged_sets;
131 CHKERR get_moab().get_entities_by_type_and_tag(0, MBENTITYSET, &part_tag,
132 tag_vals, 1, tagged_sets,
133 moab::Interface::UNION);
134 for (auto s : tagged_sets)
135 CHKERR get_moab().add_entities(s, &meshset, 1);
137 };
138
139 auto create_tags = [&](auto meshset, auto id) {
141 CHKERR get_moab().tag_set_data(th_FieldId, &meshset, 1, &id);
142 // space
143 CHKERR get_moab().tag_set_data(th_FieldSpace, &meshset, 1, &space);
144 // continuity
145 CHKERR get_moab().tag_set_data(th_FieldContinuity, &meshset, 1,
146 &continuity);
147 // base
148 CHKERR get_moab().tag_set_data(th_FieldBase, &meshset, 1, &base);
149
150 // name
151 void const *tag_data[] = {name.c_str()};
152 int tag_sizes[1];
153 tag_sizes[0] = name.size();
154 CHKERR get_moab().tag_set_by_ptr(th_FieldName, &meshset, 1, tag_data,
155 tag_sizes);
156 // rank
157 Tag th_rank;
158 int def_rank = 1;
159 const std::string tag_rank_name = "_Field_Rank_" + name;
160 CHKERR get_moab().tag_get_handle(tag_rank_name.c_str(), 1, MB_TYPE_INTEGER,
161 th_rank, MB_TAG_CREAT | MB_TAG_SPARSE,
162 &def_rank);
163 CHKERR get_moab().tag_set_data(th_rank, &meshset, 1, &nb_of_coefficients);
164
165 Version file_ver;
166 ierr = UnknownInterface::getFileVersion(moab, file_ver);
167 CHK_THROW_MESSAGE(ierr, "Not known file version");
168 if (file_ver.majorVersion >= 0 && file_ver.minorVersion >= 12 &&
169 file_ver.buildVersion >= 1) {
170
171 // Change tag names comparing to older versions
172
173 const std::string name_data_prefix("_App_Data_");
174 void const *tag_prefix_data[] = {name_data_prefix.c_str()};
175 int tag_prefix_sizes[1];
176 tag_prefix_sizes[0] = name_data_prefix.size();
177 CHKERR get_moab().tag_set_by_ptr(th_FieldName_DataNamePrefix, &meshset, 1,
178 tag_prefix_data, tag_prefix_sizes);
179 Tag th_app_order, th_field_data, th_field_data_vert;
180
181 // order
182 ApproximationOrder def_approx_order = -1;
183 const std::string tag_approx_order_name = "_App_Order_" + name;
184 CHKERR get_moab().tag_get_handle(
185 tag_approx_order_name.c_str(), 1, MB_TYPE_INTEGER, th_app_order,
186 MB_TAG_CREAT | MB_TAG_SPARSE, &def_approx_order);
187
188 // data
189 std::string tag_data_name = name_data_prefix + name;
190 const int def_len = 0;
191 CHKERR get_moab().tag_get_handle(
192 tag_data_name.c_str(), def_len, MB_TYPE_DOUBLE, th_field_data,
193 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, NULL);
194 std::string tag_data_name_verts = name_data_prefix + name + "_V";
195 VectorDouble def_vert_data(nb_of_coefficients);
196 def_vert_data.clear();
197 CHKERR get_moab().tag_get_handle(
198 tag_data_name_verts.c_str(), nb_of_coefficients, MB_TYPE_DOUBLE,
199 th_field_data_vert, MB_TAG_CREAT | tag_type, &*def_vert_data.begin());
200 } else {
201 // Deprecated names
202 const std::string name_data_prefix("_App_Data");
203 void const *tag_prefix_data[] = {name_data_prefix.c_str()};
204 int tag_prefix_sizes[1];
205 tag_prefix_sizes[0] = name_data_prefix.size();
206 CHKERR get_moab().tag_set_by_ptr(th_FieldName_DataNamePrefix, &meshset, 1,
207 tag_prefix_data, tag_prefix_sizes);
208 Tag th_app_order, th_field_data, th_field_data_vert;
209
210 // order
211 ApproximationOrder def_approx_order = -1;
212 const std::string tag_approx_order_name = "_App_Order_" + name;
213 CHKERR get_moab().tag_get_handle(
214 tag_approx_order_name.c_str(), 1, MB_TYPE_INTEGER, th_app_order,
215 MB_TAG_CREAT | MB_TAG_SPARSE, &def_approx_order);
216
217 // data
218 std::string tag_data_name = name_data_prefix + name;
219 const int def_len = 0;
220 CHKERR get_moab().tag_get_handle(
221 tag_data_name.c_str(), def_len, MB_TYPE_DOUBLE, th_field_data,
222 MB_TAG_CREAT | MB_TAG_VARLEN | MB_TAG_SPARSE, NULL);
223 std::string tag_data_name_verts = name_data_prefix + name + "V";
224 VectorDouble def_vert_data(nb_of_coefficients);
225 def_vert_data.clear();
226 CHKERR get_moab().tag_get_handle(
227 tag_data_name_verts.c_str(), nb_of_coefficients, MB_TYPE_DOUBLE,
228 th_field_data_vert, MB_TAG_CREAT | tag_type, &*def_vert_data.begin());
229 }
230
232 };
233
234 if (verb == -1)
235 verb = verbose;
236 *buildMoFEM = 0;
237 auto fit = fIelds.get<FieldName_mi_tag>().find(name);
238 if (fit != fIelds.get<FieldName_mi_tag>().end()) {
239 if (bh == MF_EXCL)
240 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
241 "field is <%s> in database", name.c_str());
242
243 } else {
244
245 EntityHandle meshset;
246 CHKERR get_moab().create_meshset(MESHSET_SET, meshset);
247 CHKERR add_meshset_to_partition(meshset);
248
249 // id
250 int field_shift = 0;
251 for (;
252 fIelds.get<BitFieldId_mi_tag>().find(BitFieldId().set(field_shift)) !=
253 fIelds.get<BitFieldId_mi_tag>().end();
254 ++field_shift) {
255 if (field_shift == BITFEID_SIZE)
256 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
257 "Maximal number of fields exceeded");
258 }
259
260 auto id = BitFieldId().set(field_shift);
261 CHKERR create_tags(meshset, id);
262
263 auto p = fIelds.insert(boost::make_shared<Field>(moab, meshset));
264 if (verb > QUIET) {
265 MOFEM_LOG("WORLD", Sev::inform) << "Add field " << **p.first;
266 MOFEM_LOG("WORLD", Sev::noisy)
267 << "Field " << (*p.first)->getName() << " core value < "
268 << this->getValue() << " > field value ( "
269 << static_cast<int>((*p.first)->getBitNumber()) << " )";
270 }
271
272 if (!p.second)
273 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
274 "field not inserted %s (top tip, it could be already "
275 "there)",
276 Field(moab, meshset).getName().c_str());
277 }
278
279 MOFEM_LOG_CHANNEL("WORLD");
280 MOFEM_LOG_CHANNEL("SYNC");
282}
283
285 const std::string name, const FieldSpace space,
286 const FieldApproximationBase base,
287 const FieldCoefficientsNumber nb_of_coefficients,
288
289 std::vector<
290
291 std::pair<EntityType,
293
294 >>
295 list_dof_side_map,
296
297 const TagType tag_type, const enum MoFEMTypes bh, int verb) {
299 CHKERR this->addField(name, space, DISCONTINUOUS, base, nb_of_coefficients,
300 tag_type, bh, verb);
301
302 auto fit = fIelds.get<FieldName_mi_tag>().find(name);
303 if (fit == fIelds.get<FieldName_mi_tag>().end())
304 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
305 "field <%s> not in database", name.c_str());
306
307 for(auto &p : list_dof_side_map) {
308 auto &data_side_map = (*fit)->getDofSideMap()[p.first];
309 CHKERR p.second(data_side_map);
310 }
311
313}
314
315MoFEMErrorCode Core::add_field(const std::string name, const FieldSpace space,
316 const FieldApproximationBase base,
317 const FieldCoefficientsNumber nb_of_coefficients,
318 const TagType tag_type, const enum MoFEMTypes bh,
319 int verb) {
320 return this->addField(name, space, CONTINUOUS, base, nb_of_coefficients,
321 tag_type, bh, verb);
322}
323
325 const std::string &name, int verb) {
326
327 *buildMoFEM = 0;
329 if (verb == -1)
330 verb = verbose;
331
332 MOFEM_LOG_CHANNEL("SYNC");
333 MOFEM_LOG_TAG("SYNC", "FieldCore");
335
337 idm = get_field_meshset(name);
338 FieldSpace space;
339 CHKERR get_moab().tag_get_data(th_FieldSpace, &idm, 1, &space);
340 FieldContinuity continuity;
341 CHKERR get_moab().tag_get_data(th_FieldContinuity, &idm, 1, &continuity);
342
343 switch (continuity) {
344 case CONTINUOUS:
345 case DISCONTINUOUS:
346 break;
347 default:
348 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
349 "sorry, unknown continuity added to entity");
350 }
351
352 std::vector<int> nb_ents_on_dim(3, 0);
353 switch (space) {
354 case L2:
355 CHKERR get_moab().add_entities(idm, ents);
356 break;
357 case H1:
358 CHKERR get_moab().add_entities(idm, ents);
359 if (continuity == DISCONTINUOUS)
360 break;
361 for (int dd = 0; dd != dim; ++dd) {
362 Range adj_ents;
363 CHKERR get_moab().get_adjacencies(ents, dd, false, adj_ents,
364 moab::Interface::UNION);
365 if (dd == 0) {
366 Range topo_nodes;
367 CHKERR get_moab().get_connectivity(ents, topo_nodes, true);
368 Range mid_nodes;
369 CHKERR get_moab().get_connectivity(ents, mid_nodes, false);
370 mid_nodes = subtract(mid_nodes, topo_nodes);
371 adj_ents = subtract(adj_ents, mid_nodes);
372 }
373 CHKERR get_moab().add_entities(idm, adj_ents);
374 nb_ents_on_dim[dd] = adj_ents.size();
375 }
376 break;
377 case HCURL:
378 CHKERR get_moab().add_entities(idm, ents);
379 if (continuity == DISCONTINUOUS)
380 break;
381 for (int dd = 1; dd != dim; ++dd) {
382 Range adj_ents;
383 CHKERR get_moab().get_adjacencies(ents, dd, false, adj_ents,
384 moab::Interface::UNION);
385 CHKERR get_moab().add_entities(idm, adj_ents);
386 nb_ents_on_dim[dd] = adj_ents.size();
387 }
388 break;
389 case HDIV:
390 CHKERR get_moab().add_entities(idm, ents);
391 if (continuity == DISCONTINUOUS)
392 break;
393 if (dim > 2) {
394 Range adj_ents;
395 CHKERR get_moab().get_adjacencies(ents, 2, false, adj_ents,
396 moab::Interface::UNION);
397 CHKERR get_moab().add_entities(idm, adj_ents);
398 nb_ents_on_dim[2] = adj_ents.size();
399 }
400 break;
401 default:
402 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
403 "sorry, unknown space added to entity");
404 }
405
406 if (verb >= VERBOSE) {
407 MOFEM_LOG("SYNC", Sev::noisy) << "add entities to field " << name;
408 MOFEM_LOG("SYNC", Sev::noisy) << "\tnb. add ents " << ents.size();
409 MOFEM_LOG("SYNC", Sev::noisy) << "\tnb. add faces " << nb_ents_on_dim[2];
410 MOFEM_LOG("SYNC", Sev::noisy) << "\tnb. add edges " << nb_ents_on_dim[1];
411 MOFEM_LOG("SYNC", Sev::noisy) << "\tnb. add nodes " << nb_ents_on_dim[0];
412 }
414}
415
417 const std::string &name,
418 int verb) {
420 Range ents_dim = ents.subset_by_dimension(dim);
421 CHKERR addEntsToFieldByDim(ents_dim, dim, name, verb);
422 MOFEM_LOG_SYNCHRONISE(mofemComm);
424}
425
427 const EntityType type,
428 const std::string &name,
429 int verb) {
431 Range ents_type = ents.subset_by_type(type);
432 if (!ents_type.empty()) {
433 const int dim = get_moab().dimension_from_handle(ents_type[0]);
434 CHKERR addEntsToFieldByDim(ents_type, dim, name, verb);
435 }
436 MOFEM_LOG_SYNCHRONISE(mofemComm);
438}
439
441 const int dim,
442 const std::string &name,
443 const bool recursive, int verb) {
445 Range ents;
446 CHKERR get_moab().get_entities_by_dimension(meshset, dim, ents, recursive);
447 CHKERR addEntsToFieldByDim(ents, dim, name, verb);
448 MOFEM_LOG_SYNCHRONISE(mofemComm);
450}
451
453 const EntityType type,
454 const std::string &name,
455 const bool recursive, int verb) {
457 Range ents;
458 CHKERR get_moab().get_entities_by_type(meshset, type, ents, recursive);
459 if (!ents.empty()) {
460 const int dim = get_moab().dimension_from_handle(ents[0]);
461 CHKERR addEntsToFieldByDim(ents, dim, name, verb);
462 }
463 MOFEM_LOG_SYNCHRONISE(mofemComm);
465}
466
468 const double coords[],
469 int size, int verb) {
471
472 if (verb == DEFAULT_VERBOSITY)
473 verb = verbose;
474
475 Range verts;
476
477 auto create_vertices = [&]() {
479
480 vector<double *> arrays_coord;
481 EntityHandle startv = 0;
482 ReadUtilIface *iface;
483 CHKERR get_moab().query_interface(iface);
484 CHKERR iface->get_node_coords(3, size, 0, startv, arrays_coord);
485 verts.insert(startv, startv + size - 1);
486 for (int n = 0; n != size; ++n)
487 for (auto d : {0, 1, 2})
488 arrays_coord[d][n] = coords[3 * n + d];
489
491 };
492
493 auto add_verts_to_field = [&]() {
495 EntityHandle field_meshset = get_field_meshset(name);
496 CHKERR get_moab().add_entities(field_meshset, verts);
498 };
499
500 CHKERR create_vertices();
501 CHKERR add_verts_to_field();
502
504}
505
506MoFEMErrorCode Core::setFieldOrderImpl(boost::shared_ptr<Field> field_ptr,
507 const Range &ents,
509 int verb) {
511
512 if (verb == DEFAULT_VERBOSITY)
513 verb = verbose;
514 *buildMoFEM = 0;
515
516 const auto field_meshset = field_ptr->getMeshset();
517 const auto bit_number = field_ptr->getBitNumber();
518
519 // intersection with field meshset
520 Range ents_of_id_meshset;
521 CHKERR get_moab().get_entities_by_handle(field_meshset, ents_of_id_meshset,
522 false);
523 Range field_ents = intersect(ents, ents_of_id_meshset);
524 if (verb > QUIET) {
525 MOFEM_LOG_C("SYNC", Sev::noisy,
526 "change nb. of ents for order in the field <%s> %d",
527 field_ptr->getName().c_str(), field_ents.size(), ents.size());
528 }
529
530 // ent view by field id (in set all MoabEnts has the same FieldId)
531 auto eiit = entsFields.get<Unique_mi_tag>().lower_bound(
532 FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
534 if (eiit != entsFields.get<Unique_mi_tag>().end()) {
535 auto hi_eiit = entsFields.get<Unique_mi_tag>().upper_bound(
536 FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
537 std::copy(eiit, hi_eiit, std::back_inserter(ents_id_view));
538 }
539
540 if (verb > QUIET)
541 MOFEM_LOG_C("SYNC", Sev::noisy,
542 "current nb. of ents in the multi index field <%s> %d",
543 field_ptr->getName().c_str(), ents_id_view.size());
544
545 // loop over ents
546 int nb_ents_set_order_up = 0;
547 int nb_ents_set_order_down = 0;
548 int nb_ents_set_order_new = 0;
549
550 FieldEntity_change_order modify_order_no_size_change(order, false);
551 FieldEntity_change_order modify_order_size_change(order, true);
552
553 for (auto pit = field_ents.const_pair_begin();
554 pit != field_ents.const_pair_end(); pit++) {
555 EntityHandle first = pit->first;
556 EntityHandle second = pit->second;
557
558 const auto type = static_cast<EntityType>(first >> MB_ID_WIDTH);
559
560 // Sanity check
561 switch (field_ptr->getSpace()) {
562 case H1:
563 break;
564 case HCURL:
565 if (type == MBVERTEX)
566 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
567 "Hcurl space on vertices makes no sense");
568
569 break;
570 case HDIV:
571 if (type == MBVERTEX)
572 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
573 "Hdiv space on vertices makes no sense");
574
575 if (type == MBEDGE)
576 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
577 "Hdiv space on edges makes no sense");
578
579 break;
580 default:
581 break;
582 }
583
584 // Entity is in database, change order only if needed
585 Range ents_in_database;
586 auto vit = ents_id_view.get<1>().lower_bound(first);
587 auto hi_vit = ents_id_view.get<1>().upper_bound(second);
588 if (order >= 0) {
589 for (; vit != hi_vit; ++vit) {
590 ents_in_database.insert(vit->get()->getEnt());
591 // entity is in database and order is changed or reset
592 const ApproximationOrder old_approximation_order =
593 (*vit)->getMaxOrder();
594
595 if (old_approximation_order != order) {
596
597 FieldEntity_multiIndex::iterator miit =
598 entsFields.get<Unique_mi_tag>().find((*vit)->getLocalUniqueId());
599
600 if ((*miit)->getMaxOrder() < order)
601 nb_ents_set_order_up++;
602 if ((*miit)->getMaxOrder() > order)
603 nb_ents_set_order_down++;
604
605 // set dofs inactive if order is reduced, and set new order to entity
606 // if order is increased (note that dofs are not build if order is
607 // increased)
608
609 bool can_change_size = true;
610 auto dit = dofsField.get<Unique_mi_tag>().lower_bound(
612 (*miit)->getEnt()));
613 if (dit != dofsField.get<Unique_mi_tag>().end()) {
614 auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(
616 (*miit)->getEnt()));
617
618 if (dit != hi_dit)
619 can_change_size = false;
620 for (; dit != hi_dit; dit++) {
621 if ((*dit)->getDofOrder() > order) {
622 bool success = dofsField.modify(dofsField.project<0>(dit),
624 if (!success)
625 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
626 "modification unsuccessful");
627 }
628 }
629 }
630
631 bool success =
632 entsFields.modify(entsFields.project<0>(miit),
633 can_change_size ? modify_order_size_change
634 : modify_order_no_size_change);
635 if (!success)
636 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
637 "modification unsuccessful");
638 }
639 }
640 }
641
642 Range new_ents = subtract(Range(first, second), ents_in_database);
643 for (Range::const_pair_iterator pit = new_ents.const_pair_begin();
644 pit != new_ents.const_pair_end(); ++pit) {
645 const auto first = pit->first;
646 const auto second = pit->second;
647 const auto ent_type = get_moab().type_from_handle(first);
648
649 if (!field_ptr->getFieldOrderTable()[ent_type])
650 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
651 "Number of degrees of freedom for entity %s for %s space on "
652 "base %s can "
653 "not be deduced",
654 moab::CN::EntityTypeName(ent_type),
655 field_ptr->getSpaceName().c_str(),
656 field_ptr->getApproxBaseName().c_str());
657
658 auto get_nb_dofs_on_order = [&](const int order) {
659 return order >= 0 ? (field_ptr->getFieldOrderTable()[ent_type])(order)
660 : 0;
661 };
662 const int nb_dofs_on_order = get_nb_dofs_on_order(order);
663 if (nb_dofs_on_order || order == -1) {
664
665 const int field_rank = field_ptr->getNbOfCoeffs();
666 const int nb_dofs = nb_dofs_on_order * field_rank;
667
668 // Entity is not in database and order is changed or reset
669 auto miit_ref_ent =
670 refinedEntities.get<Ent_mi_tag>().lower_bound(first);
671
672 auto create_tags_for_max_order = [&](const Range &ents) {
674 if (order >= 0) {
675 std::vector<ApproximationOrder> o_vec(ents.size(), order);
676 CHKERR get_moab().tag_set_data(field_ptr->th_AppOrder, ents,
677 &*o_vec.begin());
678 }
680 };
681
682 auto create_tags_for_data = [&](const Range &ents) {
684 if (order >= 0) {
685
686 if (nb_dofs > 0) {
687 if (ent_type == MBVERTEX) {
688 std::vector<FieldData> d_vec(nb_dofs * ents.size(), 0);
689 CHKERR get_moab().tag_set_data(field_ptr->th_FieldDataVerts,
690 ents, &*d_vec.begin());
691 } else {
692 std::vector<int> tag_size(ents.size(), nb_dofs);
693 std::vector<FieldData> d_vec(nb_dofs, 0);
694 std::vector<void const *> d_vec_ptr(ents.size(),
695 &*d_vec.begin());
696 CHKERR get_moab().tag_set_by_ptr(field_ptr->th_FieldData, ents,
697 &*d_vec_ptr.begin(),
698 &*tag_size.begin());
699 }
700 }
701 }
702
704 };
705
706 auto get_ents_in_ref_ent = [&](auto miit_ref_ent) {
707 auto hi = refinedEntities.get<Ent_mi_tag>().upper_bound(second);
708 Range in;
709 for (; miit_ref_ent != hi; ++miit_ref_ent)
710 in.insert(miit_ref_ent->get()->getEnt());
711 return in;
712 };
713
714 auto get_ents_max_order = [&](const Range &ents) {
715 boost::shared_ptr<std::vector<const void *>> vec(
716 new std::vector<const void *>());
717 vec->resize(ents.size());
718 CHKERR get_moab().tag_get_by_ptr(field_ptr->th_AppOrder, ents,
719 &*vec->begin());
720 return vec;
721 };
722
723 auto get_ents_field_data_vector_adaptor =
724 [&](const Range &ents,
725 boost::shared_ptr<std::vector<const void *>> &ents_max_orders) {
726 // create shared pointer and reserve memory
727 boost::shared_ptr<std::vector<double *>> vec(
728 new std::vector<double *>());
729 vec->reserve(ents.size());
730
731 if (order >= 0 && nb_dofs == 0) {
732 // set empty vector adaptor
733 for (int i = 0; i != ents.size(); ++i)
734 vec->emplace_back(nullptr);
735 } else {
736 moab::ErrorCode rval;
737 std::vector<int> tag_size(ents.size());
738 std::vector<const void *> d_vec_ptr(ents.size());
739
740 // get tags data
741 if (ent_type == MBVERTEX)
742 rval = get_moab().tag_get_by_ptr(field_ptr->th_FieldDataVerts,
743 ents, &*d_vec_ptr.begin(),
744 &*tag_size.begin());
745 else
746 rval = get_moab().tag_get_by_ptr(field_ptr->th_FieldData,
747 ents, &*d_vec_ptr.begin(),
748 &*tag_size.begin());
749
750 auto cast = [](auto p) {
751 return const_cast<FieldData *const>(
752 static_cast<const FieldData *>(p));
753 };
754
755 // some of entities has tag not set or zero dofs on entity
756 if (rval == MB_SUCCESS) {
757 // all is ok, all entities has tag set
758 auto tit = d_vec_ptr.begin();
759 auto oit = ents_max_orders->begin();
760 for (auto sit = tag_size.begin(); sit != tag_size.end();
761 ++sit, ++tit, ++oit)
762 vec->emplace_back(cast(*tit));
763
764 } else {
765 // set empty vector adaptor
766 for (int i = 0; i != ents.size(); ++i)
767 vec->emplace_back(nullptr);
768
769 // check order on all entities, and if for that order non zero
770 // dofs are expected get pointer to tag data and reset vector
771 // adaptor
772 auto oit = ents_max_orders->begin();
773 auto dit = vec->begin();
774 for (auto eit = ents.begin(); eit != ents.end();
775 ++eit, ++oit, ++dit) {
776
777 const int ent_order =
778 *static_cast<const ApproximationOrder *>(*oit);
779 const int ent_nb_dofs = get_nb_dofs_on_order(ent_order);
780
781 if (ent_nb_dofs) {
782 int tag_size;
783 const void *ret_val;
784 if (ent_type == MBVERTEX) {
785 rval = get_moab().tag_get_by_ptr(
786 field_ptr->th_FieldDataVerts, &*eit, 1, &ret_val,
787 &tag_size);
788 } else {
789 rval = get_moab().tag_get_by_ptr(
790 field_ptr->th_FieldData, &*eit, 1, &ret_val,
791 &tag_size);
792 if (rval != MB_SUCCESS) {
793
794 const int set_tag_size[] = {ent_nb_dofs * field_rank};
795 std::array<FieldData, MAX_DOFS_ON_ENTITY> set_d_vec;
796 std::fill(set_d_vec.begin(),
797 &set_d_vec[set_tag_size[0]], 0);
798 const void *set_d_vec_ptr[] = {set_d_vec.data()};
799 CHKERR get_moab().tag_set_by_ptr(
800 field_ptr->th_FieldData, &*eit, 1, set_d_vec_ptr,
801 set_tag_size);
802 rval = get_moab().tag_get_by_ptr(
803 field_ptr->th_FieldData, &*eit, 1, &ret_val,
804 &tag_size);
805
806 if (rval != MB_SUCCESS) {
810 MOFEM_LOG("SELF", Sev::error)
811 << "Error is triggered in MOAB, field tag data "
812 "for same reason can not be for accessed.";
813 MOFEM_LOG("SELF", Sev::error)
814 << "Set order: " << order;
815 MOFEM_LOG("SELF", Sev::error)
816 << "Nb. dofs on entity for given order: "
817 << set_tag_size[0];
818 MOFEM_LOG("SELF", Sev::error)
819 << "Entity type: "
820 << moab::CN::EntityTypeName(ent_type);
821 MOFEM_LOG("SELF", Sev::error)
822 << "Field: " << *field_ptr;
823 CHKERR rval;
824 }
825 }
826 }
827 const_cast<FieldData *&>(*dit) = cast(ret_val);
828 }
829 }
830 }
831 }
832 return vec;
833 };
834
835 auto ents_in_ref_ent = get_ents_in_ref_ent(miit_ref_ent);
836
837 CHKERR create_tags_for_max_order(ents_in_ref_ent);
838 CHKERR create_tags_for_data(ents_in_ref_ent);
839 auto ents_max_order = get_ents_max_order(ents_in_ref_ent);
840 auto ent_field_data =
841 get_ents_field_data_vector_adaptor(ents_in_ref_ent, ents_max_order);
842
843 // reserve memory for field dofs
844 auto ents_array = boost::make_shared<std::vector<FieldEntity>>();
845 // Add sequence to field data structure. Note that entities are
846 // allocated once into vector. This vector is passed into sequence as a
847 // weak_ptr. Vector is destroyed at the point last entity inside that
848 // vector is destroyed.
849 ents_array->reserve(second - first + 1);
850 auto vit_max_order = ents_max_order->begin();
851 auto vit_field_data = ent_field_data->begin();
852 for (int i = 0; i != ents_in_ref_ent.size(); ++i) {
853
854 ents_array->emplace_back(
855 field_ptr, *miit_ref_ent,
856 boost::shared_ptr<double *const>(ent_field_data,
857 &*vit_field_data),
858 boost::shared_ptr<const int>(
859 ents_max_order, static_cast<const int *>(*vit_max_order)));
860 ++vit_max_order;
861 ++vit_field_data;
862 ++miit_ref_ent;
863 }
864 if (!ents_array->empty())
865 if ((*ents_array)[0].getFieldRawPtr() != field_ptr.get())
866 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
867 "Get field ent poiter and field pointer do not match for "
868 "field %s",
869 field_ptr->getName().c_str());
870 nb_ents_set_order_new += ents_array->size();
871
872 // Check if any of entities in the range has bit level but is not added
873 // to database. That generate data inconsistency and error.
874 if (ents_in_ref_ent.size() < (second - first + 1)) {
875 Range ents_not_in_database =
876 subtract(Range(first, second), ents_in_ref_ent);
877 std::vector<const void *> vec_bits(ents_not_in_database.size());
878 CHKERR get_moab().tag_get_by_ptr(
879 get_basic_entity_data_ptr()->th_RefBitLevel, ents_not_in_database,
880 &*vec_bits.begin());
881 auto cast = [](auto p) {
882 return static_cast<const BitRefLevel *>(p);
883 };
884 for (auto v : vec_bits)
885 if (cast(v)->any())
886 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
887 "Try to add entities which are not seeded or added to "
888 "database");
889 }
890
891 // Add entities to database
892 auto hint = entsFields.end();
893 for (auto &v : *ents_array)
894 hint = entsFields.emplace_hint(hint, ents_array, &v);
895 }
896 }
897 }
898
899 if (verb > QUIET) {
900 MOFEM_LOG_C("SYNC", Sev::noisy,
901 "nb. of entities in field <%s> for which order was "
902 "increased %d (order %d)",
903 field_ptr->getName().c_str(), nb_ents_set_order_up, order);
904 MOFEM_LOG_C("SYNC", Sev::noisy,
905 "nb. of entities in field <%s> for which order was "
906 "reduced %d (order %d)",
907 field_ptr->getName().c_str(), nb_ents_set_order_down, order);
909 "SYNC", Sev::noisy,
910 "nb. of entities in field <%s> for which order set %d (order %d)",
911 field_ptr->getName().c_str(), nb_ents_set_order_new, order);
912 }
913
914 if (verb > QUIET) {
915 auto eiit = entsFields.get<Unique_mi_tag>().lower_bound(
916 FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
917 auto hi_eiit = entsFields.get<Unique_mi_tag>().upper_bound(
918 FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
919 MOFEM_LOG_C("SYNC", Sev::noisy,
920 "nb. of ents in the multi index field <%s> %d",
921 field_ptr->getName().c_str(), std::distance(eiit, hi_eiit));
922 }
923
924 if (verb > QUIET)
925 MOFEM_LOG_SEVERITY_SYNC(mofemComm, Sev::noisy);
926
928}
929
931 const ApproximationOrder order, int verb) {
932 MOFEM_LOG_CHANNEL("WORLD");
933 MOFEM_LOG_TAG("WORLD", "FieldCore");
934 MOFEM_LOG_CHANNEL("SYNC");
935 MOFEM_LOG_TAG("SYNC", "FieldCore");
938
939 // check field & meshset
940 auto field_it = fIelds.get<BitFieldId_mi_tag>().find(id);
941 if (field_it == fIelds.get<BitFieldId_mi_tag>().end())
942 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "no field found");
943
944 MOFEM_LOG("WORLD", Sev::noisy)
945 << "Field " << (*field_it)->getName() << " core value < "
946 << this->getValue() << " > field value ( "
947 << static_cast<int>((*field_it)->getBitNumber()) << " )"
948 << " nb. of ents " << ents.size();
949
950 CHKERR this->setFieldOrderImpl(*field_it, ents, order, verb);
951
953}
954
956 const EntityType type, const BitFieldId id,
957 const ApproximationOrder order, int verb) {
959 if (verb == -1)
960 verb = verbose;
961 *buildMoFEM = 0;
962 Range ents;
963 CHKERR get_moab().get_entities_by_type(meshset, type, ents);
964 CHKERR this->set_field_order(ents, id, order, verb);
966}
968 const EntityType type,
969 const std::string &name,
970 const ApproximationOrder order, int verb) {
972 if (verb == -1)
973 verb = verbose;
974 *buildMoFEM = 0;
975 CHKERR this->set_field_order(meshset, type, get_field_id(name), order, verb);
977}
978MoFEMErrorCode Core::set_field_order(const Range &ents, const std::string &name,
979 const ApproximationOrder order, int verb) {
981 if (verb == -1)
982 verb = verbose;
983 *buildMoFEM = 0;
984 CHKERR this->set_field_order(ents, get_field_id(name), order, verb);
986}
988 const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type,
989 const BitFieldId id, const ApproximationOrder order, int verb) {
991 if (verb == -1)
992 verb = verbose;
993 *buildMoFEM = 0;
994 Range ents;
996 ents, verb);
997 CHKERR this->set_field_order(ents, id, order, verb);
999}
1000
1002 const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type,
1003 const std::string &name, const ApproximationOrder order, int verb) {
1005 if (verb == -1)
1006 verb = verbose;
1007 *buildMoFEM = 0;
1008 Range ents;
1010 ents, verb);
1011 CHKERR this->set_field_order(ents, get_field_id(name), order, verb);
1013}
1014
1016Core::buildFieldForNoFieldImpl(boost::shared_ptr<Field> field_ptr,
1017 std::map<EntityType, int> &dof_counter,
1018 int verb) {
1020
1021 const auto bit_number = field_ptr->getBitNumber();
1022
1023 // ents in the field meshset
1024 Range ents;
1025 CHKERR get_moab().get_entities_by_handle(field_ptr->getMeshset(), ents,
1026 false);
1027 if (verb > VERBOSE)
1028 MOFEM_LOG_C("SYNC", Sev::noisy, "Ents in field %s meshset %d\n",
1029 field_ptr->getName().c_str(), ents.size());
1030
1031 // ent view by field id (in set all MoabEnts has the same FieldId)
1032 auto eiit = entsFields.get<Unique_mi_tag>().lower_bound(
1033 FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
1035 if (eiit != entsFields.get<Unique_mi_tag>().end()) {
1036 auto hi_eiit = entsFields.get<Unique_mi_tag>().upper_bound(
1037 FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
1038 std::copy(eiit, hi_eiit, std::back_inserter(ents_id_view));
1039 }
1040
1041 boost::shared_ptr<const int> zero_order(new const int(0));
1042
1043 for (auto ent : ents) {
1044 // search if field meshset is in database
1045 auto ref_ent_it = refinedEntities.get<Ent_mi_tag>().find(ent);
1046 if (ref_ent_it == refinedEntities.get<Ent_mi_tag>().end())
1047 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1048 "Entity is not in MoFEM database, entities in field meshset need "
1049 "to be seeded (i.e. bit ref level add to them)");
1050
1051 auto add_dofs = [&](auto field_eit) {
1053 // create dofs on this entity (nb. of dofs is equal to rank)
1054 for (FieldCoefficientsNumber rank = 0; rank < field_ptr->getNbOfCoeffs();
1055 rank++) {
1056 // insert dof
1057 auto p = dofsField.insert(
1058 boost::make_shared<DofEntity>(field_eit, 0, rank, rank));
1059 if (p.second) {
1060 dof_counter[MBENTITYSET]++; // Count entities in the meshset
1061 } else
1062 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1063 "Dof expected to be created");
1064 }
1066 };
1067
1068 // create database entity
1069 auto field_ent_it = ents_id_view.get<1>().find(ent);
1070 if (field_ent_it == ents_id_view.get<1>().end()) {
1071
1072 auto p = entsFields.insert(
1073
1074 boost::make_shared<FieldEntity>(
1075 field_ptr, *ref_ent_it,
1077 *ref_ent_it),
1078 boost::shared_ptr<const int>(zero_order, zero_order.get()))
1079
1080 );
1081
1082 if ((*p.first)->getFieldRawPtr() != field_ptr.get())
1083 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1084 "Get field ent poiter and field pointer do not match for "
1085 "field %s",
1086 field_ptr->getName().c_str());
1087
1088 if (!p.second)
1089 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1090 "Entity should be created");
1091
1092 CHKERR add_dofs(*(p.first));
1093
1094 } else {
1095
1096 // If there are DOFs in that range is more pragmatic to remove them
1097 // rather than to find sub-ranges or make them inactive
1098 auto dit = dofsField.get<Unique_mi_tag>().lower_bound(
1100 auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(
1102 dofsField.get<Unique_mi_tag>().erase(dit, hi_dit);
1103 CHKERR add_dofs(*field_ent_it);
1104 }
1105 }
1106
1107 if (verb > VERBOSE) {
1108 auto lo_dof = dofsField.get<Unique_mi_tag>().lower_bound(
1109 FieldEntity::getLoBitNumberUId(field_ptr->getBitNumber()));
1110 auto hi_dof = dofsField.get<Unique_mi_tag>().upper_bound(
1111 FieldEntity::getHiBitNumberUId(field_ptr->getBitNumber()));
1112 for (; lo_dof != hi_dof; lo_dof++)
1113 MOFEM_LOG("SYNC", Sev::noisy) << **lo_dof;
1114 MOFEM_LOG_SEVERITY_SYNC(mofemComm, Sev::noisy);
1115 }
1116
1118}
1119
1122 std::map<EntityType, int> &dof_counter, int verb) {
1124
1125 if (verb == -1)
1126 verb = verbose;
1127
1128 // find fields
1129 auto field_it = fIelds.get<BitFieldId_mi_tag>().find(id);
1130 if (field_it == fIelds.get<BitFieldId_mi_tag>().end())
1131 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found");
1132
1133 if (verb > QUIET)
1134 MOFEM_LOG("WORLD", Sev::noisy)
1135 << "Field " << (*field_it)->getName() << " core value < "
1136 << this->getValue() << " > field value () "
1137 << (*field_it)->getBitNumber() << " )";
1138
1139 CHKERR this->buildFieldForNoFieldImpl(*field_it, dof_counter, verb);
1140
1142}
1143
1145 const BitFieldId id, std::map<EntityType, int> &dof_counter,
1146 std::map<EntityType, int> &inactive_dof_counter, int verb) {
1148 if (verb == -1)
1149 verb = verbose;
1150
1151 // Find field
1152 auto &set_id = fIelds.get<BitFieldId_mi_tag>();
1153 auto field_it = set_id.find(id);
1154 if (field_it == set_id.end()) {
1155 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field not found");
1156 }
1157 const int bit_number = field_it->get()->getBitNumber();
1158 const int rank = field_it->get()->getNbOfCoeffs();
1159
1160 // Ents in the field meshset
1161 Range ents_of_id_meshset;
1162 CHKERR get_moab().get_entities_by_handle((*field_it)->meshSet,
1163 ents_of_id_meshset, false);
1164 if (verb > VERY_NOISY) {
1165 MOFEM_LOG_C("SYNC", Sev::noisy, "Ents in field %s meshset %d",
1166 (*field_it)->getName().c_str(), ents_of_id_meshset.size());
1167 }
1168
1169 for (auto p_eit = ents_of_id_meshset.pair_begin();
1170 p_eit != ents_of_id_meshset.pair_end(); ++p_eit) {
1171
1172 const EntityHandle first = p_eit->first;
1173 const EntityHandle second = p_eit->second;
1174 const auto lo_uid =
1175 FieldEntity::getLoLocalEntityBitNumber(bit_number, first);
1176 const auto hi_uid =
1177 FieldEntity::getHiLocalEntityBitNumber(bit_number, second);
1178
1179 auto feit = entsFields.get<Unique_mi_tag>().lower_bound(lo_uid);
1180 if (feit == entsFields.get<Unique_mi_tag>().end())
1181 continue;
1182 auto hi_feit = entsFields.get<Unique_mi_tag>().upper_bound(hi_uid);
1183
1184 // If there are DOFs in that range is more pragmatic to remove them
1185 // rather than to find sub-ranges or make them inactive
1186 auto dit = dofsField.get<Unique_mi_tag>().lower_bound(lo_uid);
1187 auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(hi_uid);
1188 dofsField.get<Unique_mi_tag>().erase(dit, hi_dit);
1189
1190 // Add vertices DOFs by bulk
1191 boost::shared_ptr<std::vector<DofEntity>> dofs_array =
1192 boost::make_shared<std::vector<DofEntity>>(std::vector<DofEntity>());
1193 // Add Sequence of DOFs to sequence container as weak_ptr
1194 int nb_dofs_on_ents = 0;
1195 for (auto tmp_feit = feit; tmp_feit != hi_feit; ++tmp_feit) {
1196 nb_dofs_on_ents += rank * tmp_feit->get()->getOrderNbDofs(
1197 tmp_feit->get()->getMaxOrder());
1198 }
1199
1200 // Reserve memory
1201 dofs_array->reserve(nb_dofs_on_ents);
1202
1203 // Create DOFs
1204 for (; feit != hi_feit; ++feit) {
1205 // Create dofs instances and shared pointers
1206 int DD = 0;
1207 // Loop orders (loop until max entity order is set)
1208 for (int oo = 0; oo <= feit->get()->getMaxOrder(); ++oo) {
1209 // Loop nb. dofs at order oo
1210 for (int dd = 0; dd < feit->get()->getOrderNbDofsDiff(oo); ++dd) {
1211 // Loop rank
1212 for (int rr = 0; rr < rank; ++rr, ++DD) {
1213 dofs_array->emplace_back(*feit, oo, rr, DD);
1214 ++dof_counter[feit->get()->getEntType()];
1215 }
1216 }
1217 }
1218 if (DD > feit->get()->getNbDofsOnEnt()) {
1219 std::ostringstream ss;
1220 ss << "rank " << rAnk << " ";
1221 ss << **feit << std::endl;
1222 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1223 "Expected number of DOFs on entity not equal to number added "
1224 "to database (DD = %d != %d = "
1225 "feit->get()->getNbDofsOnEnt())\n"
1226 "%s",
1227 DD, feit->get()->getNbDofsOnEnt(), ss.str().c_str());
1228 }
1229 }
1230
1231 // Insert into Multi-Index container
1232 int dofs_field_size0 = dofsField.size();
1233 auto hint = dofsField.end();
1234 for (auto &v : *dofs_array)
1235 hint = dofsField.emplace_hint(hint, dofs_array, &v);
1236
1237 // Add Sequence of DOFs to sequence container as weak_ptr
1238 field_it->get()->getDofSequenceContainer().push_back(dofs_array);
1239
1240 // Check data consistency
1241 if (PetscUnlikely(static_cast<int>(dofs_array.use_count()) !=
1242 static_cast<int>(dofs_array->size() + 1))) {
1243 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1244 "Wrong use count %ld != %ld", dofs_array.use_count(),
1245 dofs_array->size() + 1);
1246 }
1247 if (dofs_field_size0 + dofs_array->size() != dofsField.size()) {
1248 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1249 "Wrong number of inserted DOFs %ld != %ld", dofs_array->size(),
1250 dofsField.size() - dofs_field_size0);
1251 }
1252 }
1254}
1255
1256MoFEMErrorCode Core::buildField(const boost::shared_ptr<Field> &field,
1257 int verb) {
1259 if (verb == -1)
1260 verb = verbose;
1261 if (verb > QUIET)
1262 MOFEM_LOG("SYNC", Sev::verbose) << "Build field " << field->getName();
1263
1264 std::map<EntityType, int> dof_counter;
1265 std::map<EntityType, int> inactive_dof_counter;
1266
1267 // Need to rebuild order table since number of dofs on each order when
1268 // field was created.
1269 if (field->getApproxBase() == USER_BASE)
1270 CHKERR field->rebuildDofsOrderMap();
1271
1272 switch (field->getSpace()) {
1273 case NOFIELD:
1274 CHKERR this->buildFieldForNoField(field->getId(), dof_counter, verb);
1275 break;
1276 case L2:
1277 case H1:
1278 case HCURL:
1279 case HDIV:
1280 CHKERR this->buildFieldForL2H1HcurlHdiv(field->getId(), dof_counter,
1281 inactive_dof_counter, verb);
1282 break;
1283 default:
1284 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
1285 }
1286
1287 if (verb > QUIET) {
1288 int nb_added_dofs = 0;
1289 int nb_inactive_added_dofs = 0;
1290 for (auto const &it : dof_counter) {
1291 MOFEM_LOG("SYNC", Sev::verbose)
1292 << "Nb. of dofs (" << moab::CN::EntityTypeName(it.first) << ") "
1293 << it.second << " (inactive " << inactive_dof_counter[it.first]
1294 << ")";
1295 nb_added_dofs += it.second;
1296 nb_inactive_added_dofs += inactive_dof_counter[it.first];
1297 }
1298 MOFEM_LOG("SYNC", Sev::verbose)
1299 << "Nb. added dofs " << nb_added_dofs << " (number of inactive dofs "
1300 << nb_inactive_added_dofs << " )";
1301 }
1303}
1304
1305MoFEMErrorCode Core::build_field(const std::string field_name, int verb) {
1307 auto field_it = fIelds.get<FieldName_mi_tag>().find(field_name);
1308 if (field_it == fIelds.get<FieldName_mi_tag>().end())
1309 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND, "Field < %s > not found",
1310 field_name.c_str());
1311
1312 CHKERR this->buildField(*field_it, verb);
1313 if (verb > QUIET)
1314 MOFEM_LOG_SYNCHRONISE(mofemComm);
1316}
1317
1320 if (verb == -1)
1321 verb = verbose;
1322
1323 for (auto field : fIelds.get<BitFieldId_mi_tag>())
1324 CHKERR this->buildField(field, verb);
1325
1326 *buildMoFEM = 1 << 0;
1327 if (verb > QUIET) {
1328 MOFEM_LOG("SYNC", Sev::verbose) << "Number of dofs " << dofsField.size();
1329 MOFEM_LOG_SYNCHRONISE(mofemComm);
1330 }
1331
1333}
1334
1338 auto dit = dofsField.get<Unique_mi_tag>().lower_bound(
1339 FieldEntity::getLoBitNumberUId(get_field_bit_number((field_name))));
1340 auto hi_dit = dofsField.get<Unique_mi_tag>().upper_bound(
1341 FieldEntity::getHiBitNumberUId(get_field_bit_number(field_name)));
1342 MOFEM_LOG("SYNC", Sev::inform) << "List DOFs:";
1343 for (; dit != hi_dit; dit++)
1344 MOFEM_LOG("SYNC", Sev::inform) << *dit;
1345
1346 MOFEM_LOG_SYNCHRONISE(mofemComm);
1348}
1349
1352 MOFEM_LOG("SYNC", Sev::inform) << "List Fields:";
1353 for (auto &miit : fIelds.get<BitFieldId_mi_tag>())
1354 MOFEM_LOG("SYNC", Sev::inform) << *miit;
1355
1356 MOFEM_LOG_SYNCHRONISE(mofemComm);
1358}
1359
1360FieldEntityByUId::iterator
1362 return entsFields.get<Unique_mi_tag>().lower_bound(
1363 FieldEntity::getLoBitNumberUId(get_field_bit_number(field_name)));
1364}
1365FieldEntityByUId::iterator
1367 return entsFields.get<Unique_mi_tag>().upper_bound(
1368 FieldEntity::getHiBitNumberUId(get_field_bit_number(field_name)));
1369}
1370DofEntityByUId::iterator
1371Core::get_dofs_by_name_begin(const std::string &field_name) const {
1372 return dofsField.get<Unique_mi_tag>().lower_bound(
1373 FieldEntity::getLoBitNumberUId(get_field_bit_number(field_name)));
1374}
1375DofEntityByUId::iterator
1376Core::get_dofs_by_name_end(const std::string &field_name) const {
1377 return dofsField.get<Unique_mi_tag>().upper_bound(
1378 FieldEntity::getHiBitNumberUId(get_field_bit_number(field_name)));
1379}
1380DofEntityByUId::iterator
1382 const EntityHandle ent) const {
1383 return dofsField.get<Unique_mi_tag>().lower_bound(
1385 ent));
1386}
1387DofEntityByUId::iterator
1389 const EntityHandle ent) const {
1390 return dofsField.get<Unique_mi_tag>().upper_bound(
1392 ent));
1393}
1394DofEntityByUId::iterator
1396 const EntityType type) const {
1397 return dofsField.get<Unique_mi_tag>().lower_bound(
1399 get_id_for_min_type(type)));
1400}
1401DofEntityByUId::iterator
1403 const EntityType type) const {
1404 return dofsField.get<Unique_mi_tag>().upper_bound(
1406 get_id_for_max_type(type)));
1407}
1409Core::check_number_of_ents_in_ents_field(const std::string &name) const {
1411 auto it = fIelds.get<FieldName_mi_tag>().find(name);
1412 if (it == fIelds.get<FieldName_mi_tag>().end()) {
1413 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1414 "field not found < %s >", name.c_str());
1415 }
1416 EntityHandle meshset = (*it)->getMeshset();
1417 int num_entities;
1418 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1419
1420 auto count_field_ents = [&]() {
1421 auto bit_number = (*it)->getBitNumber();
1422 auto low_eit = entsFields.get<Unique_mi_tag>().lower_bound(
1423 FieldEntity::getLoBitNumberUId(bit_number));
1424 auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
1425 FieldEntity::getHiBitNumberUId(bit_number));
1426 return std::distance(low_eit, hi_eit);
1427 };
1428
1429 if (count_field_ents() > (unsigned int)num_entities) {
1430 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1431 "not equal number of entities in meshset and field multiindex "
1432 "< %s >",
1433 name.c_str());
1434 }
1436}
1439 for (auto &it : fIelds.get<FieldName_mi_tag>()) {
1440 if (it->getSpace() == NOFIELD)
1441 continue; // FIXME: should be treated properly, not test is just
1442 // skipped for this NOFIELD space
1443 EntityHandle meshset = it->getMeshset();
1444 int num_entities;
1445 CHKERR get_moab().get_number_entities_by_handle(meshset, num_entities);
1446
1447 auto count_field_ents = [&]() {
1448 auto bit_number = it->getBitNumber();
1449 auto low_eit = entsFields.get<Unique_mi_tag>().lower_bound(
1450 FieldEntity::getLoBitNumberUId(bit_number));
1451 auto hi_eit = entsFields.get<Unique_mi_tag>().upper_bound(
1452 FieldEntity::getHiBitNumberUId(bit_number));
1453 return std::distance(low_eit, hi_eit);
1454 };
1455
1456 if (count_field_ents() > (unsigned int)num_entities) {
1457 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1458 "not equal number of entities in meshset and field "
1459 "multiindex < %s >",
1460 it->getName().c_str());
1461 }
1462 }
1464}
1465
1466} // namespace MoFEM
#define FieldCoreFunctionBegin
Definition FieldCore.cpp:8
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
@ QUIET
@ DEFAULT_VERBOSITY
@ VERBOSE
@ VERY_NOISY
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
@ MF_EXIST
@ MF_EXCL
FieldApproximationBase
approximation base
Definition definitions.h:58
@ USER_BASE
user implemented approximation base
Definition definitions.h:68
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
Definition definitions.h:82
@ L2
field with C-1 continuity
Definition definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition definitions.h:84
@ H1
continuous field
Definition definitions.h:85
@ HCURL
field with continuous tangents
Definition definitions.h:86
@ HDIV
field with continuous normal traction
Definition definitions.h:87
#define MB_ID_WIDTH
FieldContinuity
Field continuity.
Definition definitions.h:99
@ CONTINUOUS
Regular field.
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#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 BITFEID_SIZE
max number of finite elements
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
constexpr int order
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
#define MOFEM_LOG_FUNCTION()
Set scope.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
int ApproximationOrder
Approximation on the entity.
Definition Types.hpp:26
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition Types.hpp:42
int FieldCoefficientsNumber
Number of field coefficients.
Definition Types.hpp:27
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
char FieldBitNumber
Field bit number.
Definition Types.hpp:28
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
EntityHandle get_id_for_max_type()
EntityHandle get_id_for_min_type()
Field_multiIndex::index< FieldName_mi_tag >::type::iterator field_it
static int getMaxOrder(const ENTMULTIINDEX &multi_index)
multi_index_container< boost::shared_ptr< FieldEntity >, indexed_by< sequenced<>, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< FieldEntity::interface_type_RefEntity, EntityHandle, &FieldEntity::getEnt > > > > FieldEntity_multiIndex_ent_view
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition Common.hpp:12
constexpr auto field_name
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > > >, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Managing BitRefLevels.
std::string get_field_name(const BitFieldId id) const
get field name from id
Definition FieldCore.cpp:37
MoFEMErrorCode set_field_order_by_entity_type_and_bit_ref(const BitRefLevel &bit, const BitRefLevel &mask, const EntityType type, const BitFieldId id, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)
const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const
get field structure
Definition FieldCore.cpp:65
MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)
Add entities to field meshset.
virtual MoFEMErrorCode add_broken_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_coefficients, std::vector< std::pair< EntityType, std::function< MoFEMErrorCode(BaseFunction::DofsSideMap &)> > > list_dof_side_map, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
Add field.
bool check_field(const std::string &name) const
check if field is in database
Definition FieldCore.cpp:57
DofEntityByUId::iterator get_dofs_by_name_and_type_end(const std::string &field_name, const EntityType ent) const
get begin iterator of field dofs of given name end ent type(instead you can use IT_GET_DOFS_FIELD_BY_...
MoFEMErrorCode get_field_entities_by_dimension(const std::string name, int dim, Range &ents) const
get entities in the field by dimension
Definition FieldCore.cpp:81
MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)
Add entities to field meshset.
MoFEMErrorCode buildFieldForNoFieldImpl(boost::shared_ptr< Field > field_ptr, std::map< EntityType, int > &dof_counter, int verb)
EntityHandle get_field_meshset(const BitFieldId id) const
Definition FieldCore.cpp:45
MoFEMErrorCode create_vertices_and_add_to_field(const std::string name, const double coords[], int size, int verb=DEFAULT_VERBOSITY)
Create a vertices and add to field object.
MoFEMErrorCode set_field_order(const Range &ents, const BitFieldId id, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)
DofEntityByUId::iterator get_dofs_by_name_and_type_begin(const std::string &field_name, const EntityType type) const
get begin iterator of field dofs of given name and ent type (instead you can use IT_GET_DOFS_FIELD_BY...
FieldEntityByUId::iterator get_ent_field_by_name_begin(const std::string &field_name) const
get begin iterator of field ents of given name (instead you can use IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP...
MoFEMErrorCode addEntsToFieldByDim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)
FieldBitNumber get_field_bit_number(const std::string name) const
get field bit number
Definition FieldCore.cpp:28
BitFieldId get_field_id(const std::string &name) const
Get field Id.
Definition FieldCore.cpp:18
DofEntityByUId::iterator get_dofs_by_name_and_ent_begin(const std::string &field_name, const EntityHandle ent) const
get begin iterator of field dofs of given name and ent(instead you can use IT_GET_DOFS_FIELD_BY_NAME_...
DofEntityByUId::iterator get_dofs_by_name_and_ent_end(const std::string &field_name, const EntityHandle ent) const
get begin iterator of field dofs of given name and ent (instead you can use IT_GET_DOFS_FIELD_BY_NAME...
MoFEMErrorCode check_number_of_ents_in_ents_field() const
check data consistency in entitiesPtr
MoFEMErrorCode build_field(const std::string field_name, int verb=DEFAULT_VERBOSITY)
build field by name
DofEntityByUId::iterator get_dofs_by_name_begin(const std::string &field_name) const
get begin iterator of field dofs of given name (instead you can use IT_GET_DOFS_FIELD_BY_NAME_FOR_LOO...
MoFEMErrorCode addField(const std::string &name, const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, const FieldCoefficientsNumber nb_coefficients, const TagType tag_type, const enum MoFEMTypes bh, int verb)
Template for add_field.
MoFEMErrorCode setFieldOrderImpl(boost::shared_ptr< Field > field_ptr, const Range &ents, const ApproximationOrder order, int verb)
MoFEMErrorCode buildField(const boost::shared_ptr< Field > &field, int verb=DEFAULT_VERBOSITY)
MoFEMErrorCode list_dofs_by_field_name(const std::string &name) const
MoFEMErrorCode buildFieldForNoField(const BitFieldId id, std::map< EntityType, int > &dof_counter, int verb=DEFAULT_VERBOSITY)
DofEntityByUId::iterator get_dofs_by_name_end(const std::string &field_name) const
get begin iterator of field dofs of given name (instead you can use IT_GET_DOFS_FIELD_BY_NAME_FOR_LOO...
FieldEntityByUId::iterator get_ent_field_by_name_end(const std::string &field_name) const
get begin iterator of field dofs of given name (instead you can use IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP...
MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)
MoFEMErrorCode get_field_entities_by_handle(const std::string name, Range &ents) const
get entities in the field by handle
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)
Add field.
MoFEMErrorCode buildFieldForL2H1HcurlHdiv(const BitFieldId id, std::map< EntityType, int > &dof_counter, std::map< EntityType, int > &inactive_dof_counter, int verb=DEFAULT_VERBOSITY)
MoFEMErrorCode list_fields() const
list entities in the field
MoFEMErrorCode get_field_entities_by_type(const std::string name, EntityType type, Range &ents) const
get entities in the field by type
Definition FieldCore.cpp:91
const int getValue() const
Definition Core.hpp:39
structure to change FieldEntity order
static UId getHiBitNumberUId(const FieldBitNumber bit_number)
static boost::shared_ptr< FieldData *const > makeSharedFieldDataAdaptorPtr(const boost::shared_ptr< Field > &field_ptr, const boost::shared_ptr< RefEntity > &ref_ents_ptr)
Return shared pointer to entity field data vector adaptor.
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)
MultiIndex Tag for field name.
Provide data structure for (tensor) field approximation.
static MoFEMErrorCode getFileVersion(moab::Interface &moab, Version &version)
Get database major version.