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