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