v0.15.0
Loading...
Searching...
No Matches
Namespaces | Macros
FieldCore.cpp File Reference

Core interface methods for managing fields. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Namespaces

namespace  MoFEM
 implementation of Data Operators for Forces and Sources
 

Macros

#define FieldCoreFunctionBegin
 

Detailed Description

Core interface methods for managing fields.

Definition in file FieldCore.cpp.

Macro Definition Documentation

◆ FieldCoreFunctionBegin

#define FieldCoreFunctionBegin
Value:
MoFEMFunctionBegin; \
MOFEM_LOG_CHANNEL("WORLD"); \
MOFEM_LOG_CHANNEL("SYNC"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("WORLD", "FieldCore"); \
MOFEM_LOG_TAG("SYNC", "FieldCore");
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359

Definition at line 8 of file FieldCore.cpp.

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