v0.15.0
Loading...
Searching...
No Matches
MedInterface.cpp
Go to the documentation of this file.
1/** \file MedInterface.cpp
2 * \brief Med file interface interface
3 *
4 * Interface loading mesh and data on mesh directly to mofem & moab
5 *
6 */
7
8#ifdef WITH_MED
9
10extern "C" {
11#include <med.h>
12}
13
14#if (MED_MAJOR_NUM >= 3)
15// To avoid too many ifdefs below we use defines for the bits of the
16// API that did not change too much between MED2 and MED3. If we
17// remove MED2 support at some point, please remove these defines and
18// replace the symbols accordingly.
19#define med_geometrie_element med_geometry_type
20#define med_entite_maillage med_entity_type
21#define med_type_champ med_field_type
22#define MED_LECTURE MED_ACC_RDONLY
23#define MED_LECTURE_AJOUT MED_ACC_RDEXT
24#define MED_NOEUD MED_NODE
25#define MED_MAILLE MED_CELL
26#define MED_NOEUD_MAILLE MED_NODE_ELEMENT
27#define MED_NOPFL MED_NO_PROFILE
28#define MEDouvrir MEDfileOpen
29#define MEDfermer MEDfileClose
30#define MEDnChamp MEDfieldnComponent
31#define MEDnValProfil MEDprofileSizeByName
32#define MEDfichDesEcr MEDfileCommentWr
33#define MED_ARETE MED_EDGE
34#define MED_FACETTE MED_FACE
35#else
36#error "MED file is not MED_MAJOR_NUM == 3"
37#endif
38
39#include <MedInterface.hpp>
40#include <unordered_set>
41
42namespace MoFEM {
43
45MedInterface::query_interface(boost::typeindex::type_index type_index,
46 UnknownInterface **iface) const {
47 *iface = const_cast<MedInterface *>(this);
48 return 0;
49}
50
52 : cOre(const_cast<Core &>(core)), flgFile(PETSC_FALSE) {
53
54 if (!LogManager::checkIfChannelExist("MEDWORLD")) {
55 auto core_log = logging::core::get();
56 core_log->add_sink(
58 LogManager::setLog("MEDWORLD");
59 MOFEM_LOG_TAG("MEDWORLD", "MED");
60 }
61
62 MOFEM_LOG("MEDWORLD", Sev::noisy) << "Mashset manager created";
63}
64
66 Interface &m_field = cOre;
67 char mesh_file_name[255];
69 PetscOptionsBegin(m_field.get_comm(), "", "MED Interface", "none");
70 ierr = PetscOptionsString("-med_file", "med file name", "", "mesh.med",
71 mesh_file_name, 255, &flgFile);
73 PetscOptionsEnd();
74 if (flgFile == PETSC_TRUE) {
75 medFileName = std::string(mesh_file_name);
76 } else {
77 medFileName = std::string("mesh.med");
78 }
80}
81
82MoFEMErrorCode MedInterface::medGetFieldNames(const string &file, int verb) {
83 Interface &m_field = cOre;
85 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
86 if (fid < 0) {
87 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
88 "Unable to open file '%s'", file.c_str());
89 }
90 med_int num_fields = MEDnField(fid);
91 for (int index = 0; index < num_fields; index++) {
92
93 med_int num_comp = MEDfieldnComponent(fid, index + 1);
94 if (num_comp <= 0) {
95 SETERRQ(m_field.get_comm(), MOFEM_IMPOSSIBLE_CASE,
96 "Could not get number of components for MED field");
97 }
98
99 char name[MED_NAME_SIZE + 1], mesh_name[MED_NAME_SIZE + 1];
100 char dt_unit[MED_SNAME_SIZE + 1];
101 std::vector<char> comp_name(num_comp * MED_SNAME_SIZE + 1);
102 std::vector<char> comp_unit(num_comp * MED_SNAME_SIZE + 1);
103 med_int num_steps = 0;
104 med_type_champ type;
105 med_bool local_mesh;
106 if (MEDfieldInfo(fid, index + 1, name, mesh_name, &local_mesh, &type,
107 &comp_name[0], &comp_unit[0], dt_unit, &num_steps) < 0) {
108 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
109 "Could not get MED field info");
110 }
111
112 std::string field_name = std::string(name);
114 fieldNames[field_name].fieldName = field_name;
115 fieldNames[field_name].meshName = std::string(mesh_name);
116 fieldNames[field_name].localMesh = (local_mesh == MED_TRUE);
117 for (int ff = 0; ff != num_comp; ff++) {
118 fieldNames[field_name].componentNames.push_back(
119 std::string(&comp_name[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
120 fieldNames[field_name].componentUnits.push_back(
121 std::string(&comp_unit[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
122 }
123 fieldNames[field_name].dtUnit = std::string(&dt_unit[0]);
124 fieldNames[field_name].ncSteps = num_steps;
125
126 if (verb > 0)
127 MOFEM_LOG("MEDWORLD", Sev::inform) << fieldNames[name];
128 }
129 if (MEDfileClose(fid) < 0) {
130 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
131 "Unable to close file '%s'", file.c_str());
132 }
134}
135
144
145MoFEMErrorCode MedInterface::readMed(const string &file, int verb) {
146 Interface &m_field = cOre;
148
149 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
150 if (fid < 0) {
151 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
152 "Unable to open file '%s'", file.c_str());
153 }
154
155 med_int v[3], vf[3];
156 MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
157 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
158
159 if (verb > 0)
160 MOFEM_LOG_C("MEDWORLD", Sev::inform,
161 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
162 vf[1], vf[2], v[0], v[1], v[2]);
163
164 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
165 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
166 "Cannot read MED file older than V2.2");
167 }
168
169 for (int i = 0; i < MEDnMesh(fid); i++) {
170 char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
171 bzero(mesh_name, MED_NAME_SIZE);
172 bzero(mesh_desc, MED_COMMENT_SIZE);
173 med_int space_dim;
174 med_mesh_type mesh_type;
175 med_int mesh_dim, n_step;
176 char dt_unit[MED_SNAME_SIZE + 1];
177 char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
178 med_sorting_type sorting_type;
179 med_axis_type axis_type;
180 if (MEDmeshInfo(fid, i + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
181 mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
182 axis_name, axis_unit) < 0) {
183 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
184 "Unable to read mesh information");
185 }
186 meshNames.push_back(std::string(mesh_name));
187 if (verb > 0) {
188 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Check mesh %s nsteps %d", mesh_name,
189 n_step);
190 }
191 }
192
193 std::map<int, Range> family_elem_map;
194 std::map<string, Range> group_elem_map;
195
196 for (unsigned int ii = 0; ii != meshNames.size(); ii++) {
197 CHKERR readMesh(file, ii, family_elem_map, verb);
198 CHKERR readFamily(file, ii, family_elem_map, group_elem_map, verb);
199 CHKERR makeBlockSets(group_elem_map, verb);
200 }
201
202 if (MEDfileClose(fid) < 0) {
203 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
204 "Unable to close file '%s'", file.c_str());
205 }
206
208}
209
210static std::vector<med_geometrie_element>
211moab2med_element_type(const EntityType type) {
212
213 std::vector<med_geometrie_element> types;
214
215 switch (type) {
216 case MBEDGE:
217 types.push_back(MED_SEG2);
218 types.push_back(MED_SEG3);
219 break;
220 case MBTRI:
221 types.push_back(MED_TRIA3);
222 types.push_back(MED_TRIA6);
223 break;
224 case MBQUAD:
225 types.push_back(MED_QUAD4);
226 types.push_back(MED_QUAD8);
227 break;
228 case MBTET:
229 types.push_back(MED_TETRA4);
230 types.push_back(MED_TETRA10);
231 break;
232 case MBHEX:
233 types.push_back(MED_HEXA8);
234 types.push_back(MED_HEXA20);
235 break;
236 case MBPRISM:
237 types.push_back(MED_PENTA6);
238 break;
239 case MBPYRAMID:
240 types.push_back(MED_PYRA5);
241 break;
242 case MBVERTEX:
243 // Not Currently Used
244 //types.push_back(MED_POINT1);
245 break;
246 default:
247 break;
248 }
249
250 return types;
251}
252
253MoFEMErrorCode MedInterface::readMesh(const string &file, const int index,
254 std::map<int, Range> &family_elem_map,
255 int verb) {
256
257 Interface &m_field = cOre;
259
260 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
261 if (fid < 0) {
262 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
263 "Unable to open file '%s'", file.c_str());
264 }
265 med_int v[3], vf[3];
266 MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
267 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
268 if (verb > 1)
269 MOFEM_LOG_C("MEDWORLD", Sev::noisy,
270 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
271 vf[1], vf[2], v[0], v[1], v[2]);
272
273 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
274 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
275 "Cannot read MED file older than V2.2");
276 }
277
278 char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
279 bzero(mesh_name, MED_NAME_SIZE + 1);
280 bzero(mesh_desc, MED_COMMENT_SIZE + 1);
281 med_int space_dim;
282 med_mesh_type mesh_type;
283 med_int mesh_dim, n_step;
284 char dt_unit[MED_SNAME_SIZE + 1];
285 char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
286 med_sorting_type sorting_type;
287 med_axis_type axis_type;
288 if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
289 mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
290 axis_name, axis_unit) < 0) {
291 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
292 "Unable to read mesh information");
293 }
294 if (verb > 0)
295 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Reading mesh %s nsteps %d", mesh_name,
296 n_step);
297
298 switch (axis_type) {
299 case MED_CARTESIAN:
300 break;
301 case MED_SPHERICAL:
302 case MED_CYLINDRICAL:
303 default:
304 SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
305 "Curvilinear coordinate system implemented");
306 }
307 if (space_dim < 2) {
308 SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
309 "Not implemented for space dim %d", space_dim);
310 }
311
312 EntityHandle mesh_meshset;
313 {
314 MeshsetsManager *meshsets_manager_ptr;
315 CHKERR m_field.getInterface(meshsets_manager_ptr);
316 int max_id = 0;
318 max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
319 }
320 max_id++;
321 CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id,
322 std::string(mesh_name));
323 CubitMeshSet_multiIndex::index<
325 cit =
326 meshsets_manager_ptr->getMeshsetsMultindex()
328 .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
329 max_id++;
330 mesh_meshset = cit->getMeshset();
331 meshMeshsets.push_back(mesh_meshset);
332 }
333
334 med_bool change_of_coord, geo_transform;
335 med_int num_nodes = MEDmeshnEntity(
336 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
337 MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
338 if (num_nodes < 0) {
339 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
340 "Could not read number of MED nodes");
341 }
342 if (num_nodes == 0) {
343 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
344 "No nodes in MED mesh");
345 }
346 if (verb > 0)
347 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Read number of nodes %d", num_nodes);
348
349 std::vector<med_float> coord_med(space_dim * num_nodes);
350 if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
351 MED_NO_INTERLACE, &coord_med[0]) < 0) {
352 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
353 "Could not read MED node coordinates");
354 }
355
356 // Add vertices to moab
357 ReadUtilIface *iface;
358 vector<double *> arrays_coord;
359 EntityHandle startv;
360 CHKERR m_field.get_moab().query_interface(iface);
361 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
362 Range verts(startv, startv + num_nodes - 1);
363 std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
364 arrays_coord[0]);
365 std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
366 arrays_coord[1]);
367 if (space_dim == 2) {
368 std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
369 } else {
370 std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
371 arrays_coord[2]);
372 }
373 CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
374 family_elem_map.clear();
375
376 // get family for vertices
377 {
378 std::vector<med_int> nodes_tags(num_nodes, 0);
379 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
380 MED_NODE, MED_NO_GEOTYPE,
381 &nodes_tags[0]) < 0) {
382 nodes_tags.clear();
383 }
384 for (int i = 0; i < num_nodes; i++) {
385 family_elem_map[nodes_tags.empty() ? i : nodes_tags[i]].insert(verts[i]);
386 }
387 }
388
389 // read elements (loop over all possible MSH element types)
390 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
391
392 auto types = moab2med_element_type(ent_type);
393
394 for (auto type : types) {
395
396 // get number of cells
397 med_bool change_of_coord;
398 med_bool geo_transform;
399 med_int num_ele = MEDmeshnEntity(
400 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL, type,
401 MED_CONNECTIVITY, MED_NODAL, &change_of_coord, &geo_transform);
402
403 // get connectivity
404 if (num_ele <= 0)
405 continue;
406
407 int num_nod_per_ele = type % 100;
408 if (verb > 0)
409 MOFEM_LOG("MEDWORLD", Sev::inform)
410 << "Reading elements " << num_ele << " of type "
411 << moab::CN::EntityTypeName(ent_type) << " number of nodes "
412 << num_nod_per_ele;
413
414 std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
415 if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
416 MED_CELL, type, MED_NODAL,
417 MED_FULL_INTERLACE, &conn_med[0]) < 0) {
418 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
419 "Could not read MED elements");
420 }
421
422 Range ents;
423
424 if (ent_type != MBVERTEX) {
425 EntityHandle *conn_moab;
426 EntityHandle starte;
427 CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
428 starte, conn_moab);
429 switch (ent_type) {
430 // FIXME: Some connectivity could not work, need to run and test
431 case MBTET: {
432 int ii = 0;
433 for (int ee = 0; ee != num_ele; ee++) {
434 std::vector<EntityHandle> n(num_nod_per_ele);
435 for (int nn = 0; nn != num_nod_per_ele; nn++) {
436 n[nn] = verts[conn_med[ii + nn] - 1];
437 }
438
439 std::array<int, 10> nodes_map{
440
441 1, 0, 2, 3,
442
443 4, 6, 5, 8,
444 7, 9
445
446 };
447
448 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
449 conn_moab[ii] = n[nodes_map[nn]];
450 }
451 }
452 }
453 case MBHEX: {
454 int ii = 0;
455 for (int ee = 0; ee != num_ele; ee++) {
456 std::vector<EntityHandle> n(num_nod_per_ele);
457 for (int nn = 0; nn != num_nod_per_ele; nn++) {
458 n[nn] = verts[conn_med[ii + nn] - 1];
459 }
460
461 std::array<int, 20> nodes_map{
462 0, 3, 2, 1, 4, 7, 6, 5, // corner nodes
463 11, 10, 9, 8, // edges on bottom face
464 16, 19, 18, 17, // vertical edges
465 15, 14, 13, 12 // edges on top face
466
467 };
468
469 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
470 conn_moab[ii] = n[nodes_map[nn]];
471 }
472 }
473 } break;
474 default: {
475 int ii = 0;
476 for (int ee = 0; ee != num_ele; ee++) {
477 for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
478 conn_moab[ii] = verts[conn_med[ii] - 1];
479 }
480 }
481 }
482 }
483 CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
484 conn_moab);
485 ents = Range(starte, starte + num_ele - 1);
486 } else {
487 // This is special case when in med vertices are defined as elements
488 int ii = 0;
489 std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
490 for (int ee = 0; ee != num_ele; ++ee)
491 for (int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
492 conn_moab[ii] = verts[conn_med[ii] - 1];
493 ents.insert_list(conn_moab.begin(), conn_moab.end());
494 }
495
496 // Add elements to family meshset
497 CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
498
499 // get family for cells
500 {
501 std::vector<med_int> fam(num_ele, 0);
502 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
503 MED_CELL, type, &fam[0]) < 0) {
504 MOFEM_LOG("MEDWORLD", Sev::warning) <<
505 "No family number for elements: using 0 as default family number";
506 }
507 for (int j = 0; j < num_ele; j++) {
508 family_elem_map[fam[j]].insert(ents[j]);
509 }
510 }
511 }
512 }
513
514 if (MEDfileClose(fid) < 0) {
515 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
516 "Unable to close file '%s'", file.c_str());
517 }
518
520}
521
523MedInterface::readFamily(const string &file, const int index,
524 const std::map<int, Range> &family_elem_map,
525 std::map<string, Range> &group_elem_map, int verb) {
526 Interface &m_field = cOre;
528
529 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
530 if (fid < 0) {
531 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
532 "Unable to open file '%s'", file.c_str());
533 }
534 med_int v[3], vf[3];
535 MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
536 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
537
538 if (verb > 1) {
539 MOFEM_LOG_C("MEDWORLD", Sev::noisy,
540 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
541 vf[1], vf[2], v[0], v[1], v[2]);
542 }
543 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
544 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
545 "Cannot read MED file older than V2.2");
546 }
547
548 // clear groups
549 group_elem_map.clear();
550
551 med_int num_families = MEDnFamily(fid, meshNames[index].c_str());
552 if (num_families < 0) {
553 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
554 "Could not read MED families");
555 }
556 for (int i = 0; i < num_families; i++) {
557 med_int num_attrib =
558 (vf[0] == 2)
559 ? MEDnFamily23Attribute(fid, meshNames[index].c_str(), i + 1)
560 : 0;
561 med_int num_groups = MEDnFamilyGroup(fid, meshNames[index].c_str(), i + 1);
562 if (num_attrib < 0 || num_groups < 0) {
563 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
564 "Could not read MED groups or attributes");
565 }
566
567 std::vector<med_int> attribId(num_attrib + 1);
568 std::vector<med_int> attrib_val(num_attrib + 1);
569 std::vector<char> attrib_des(MED_COMMENT_SIZE * num_attrib + 1);
570 std::vector<char> group_names(MED_LNAME_SIZE * num_groups + 1);
571 char family_name[MED_NAME_SIZE + 1];
572 med_int family_num;
573
574 if (vf[0] == 2) { // MED2 file
575 if (MEDfamily23Info(fid, meshNames[index].c_str(), i + 1, family_name,
576 &attribId[0], &attrib_val[0], &attrib_des[0],
577 &family_num, &group_names[0]) < 0) {
578 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
579 "Could not read info for MED2 family %d", i + 1);
580 }
581 } else {
582 if (MEDfamilyInfo(fid, meshNames[index].c_str(), i + 1, family_name,
583 &family_num, &group_names[0]) < 0) {
584 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
585 "Could not read info for MED3 family %d", i + 1);
586 }
587 }
588
589 for (int g = 0; g != num_groups; g++) {
590 std::string name =
591 std::string(&group_names[MED_LNAME_SIZE * g], MED_LNAME_SIZE - 1);
592 name.resize(NAME_TAG_SIZE - 1);
593 if (family_elem_map.find(family_num) == family_elem_map.end()) {
595 "MEDWORLD", Sev::warning,
596 "Family %d not read, likely type of element is not added "
597 "to moab database. Currently only triangle, quad, tetrahedral and "
598 "hexahedral elements are read to moab database",
599 family_num);
600 } else {
601 group_elem_map[name].merge(family_elem_map.at(family_num));
602 }
603 }
604 }
605
606 if (MEDfileClose(fid) < 0) {
607 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
608 "Unable to close file '%s'", file.c_str());
609 }
610
612}
613
615MedInterface::makeBlockSets(const std::map<string, Range> &group_elem_map,
616 int verb) {
617
618 Interface &m_field = cOre;
620 MeshsetsManager *meshsets_manager_ptr;
621 CHKERR m_field.getInterface(meshsets_manager_ptr);
622
623 int max_id = 0;
625 max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
626 }
627 max_id++;
628
629 for (std::map<string, Range>::const_iterator git = group_elem_map.begin();
630 git != group_elem_map.end(); git++) {
631 CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id, git->first);
632 CubitMeshSet_multiIndex::index<
634 cit =
635 meshsets_manager_ptr->getMeshsetsMultindex()
637 .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
638 EntityHandle meshsets = cit->getMeshset();
639 if (!git->second.empty()) {
640 CHKERR m_field.get_moab().add_entities(meshsets, git->second);
641 }
642 max_id++;
643 }
644
646 MOFEM_LOG("MEDWORLD", Sev::verbose) << *cit;
647
649}
650
659
661 boost::shared_ptr<std::vector<const CubitMeshSets *>> &meshsets_ptr) {
662 Interface &m_field = cOre;
664
665 auto &meshsets_idx =
666 m_field.getInterface<MeshsetsManager>()->getMeshsetsMultindex();
667 std::vector<const CubitMeshSets *> meshsets;
668
669 for (auto &m : meshsets_idx) {
670 meshsets.push_back(&m);
671 }
672 // Sort meshsets based on meshsetId
673 std::sort(meshsets.begin(), meshsets.end(),
674 [](const CubitMeshSets *a, const CubitMeshSets *b) {
675 return a->getMeshsetId() < b->getMeshsetId();
676 });
677
678 meshsets_ptr =
679 boost::make_shared<std::vector<const CubitMeshSets *>>(meshsets);
681}
682
683MoFEMErrorCode MedInterface::writeMed(boost::shared_ptr<Range> range_ptr,
684 int verb) {
686 if (medFileName.empty()) {
688 }
689
690 // Get all meshsets
691 boost::shared_ptr<std::vector<const CubitMeshSets *>> meshsets_ptr;
692 getMeshsets(meshsets_ptr);
693
694 CHKERR writeMed(medFileName, meshsets_ptr, range_ptr, verb);
696}
697
699 const string &file,
700 boost::shared_ptr<std::vector<const CubitMeshSets *>> meshsets_ptr,
701 boost::shared_ptr<Range> write_range_ptr, int verb) {
702 Interface &m_field = cOre;
703 // Get the MOAB instance from the field
704 moab::Interface &moab = m_field.get_moab();
705
707 MOFEM_LOG("MEDWORLD", Sev::warning)
708 << "WRITE_MED IS EXPERIMENTAL, MAY CONTAIN BUGS, ALWAYS CHECK THE OUTPUT "
709 "FILE";
710 // Open a med file with the specified version
711 med_idt fid =
712 MEDfileVersionOpen((char *)file.c_str(), MED_ACC_CREAT, MED_MAJOR_NUM,
713 MED_MINOR_NUM, MED_RELEASE_NUM);
714 // Throw an error if the file could not be opened
715 if (fid < 0) {
716 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
717 "Cannot create MED file");
718 return 0;
719 }
720
721 // Write appropriate header
722 CHKERR MEDfichDesEcr(fid, (char *)"MED file generated by MoFEM");
723
724 // Generate empty mesh
725 char dtUnit[MED_SNAME_SIZE + 1] = ""; // time unit
726 char axisName[3 * MED_SNAME_SIZE + 1] = ""; // axis name
727 char axisUnit[3 * MED_SNAME_SIZE + 1] = ""; // axis unit
728
729 PetscBool is_cubit_meshset = PETSC_TRUE;
730 int med_mesh_name_id = 0;
731
732 // Create the mesh
733 char mesh_name_char[255];
734 std::string mesh_name = "Mesh";
735
736 // Get mesh_name from command line
737 PetscOptionsBegin(PETSC_COMM_WORLD, "", "MED mesh options", "");
738 CHKERR PetscOptionsString("-med_mesh_name", "get med mesh name", "",
739 mesh_name.c_str(), mesh_name_char, 255, PETSC_NULLPTR);
740 PetscOptionsEnd();
741
742 mesh_name = mesh_name_char;
743
744 // check if mesh_name is a meshset
745 for (auto &m : *meshsets_ptr) {
746 if (m->getName() == mesh_name) {
747 med_mesh_name_id = m->getMeshsetId();
748 is_cubit_meshset = PETSC_FALSE;
749 break;
750 }
751 }
752
753 // Get the maximum meshset id
754 int max_id = 0;
755 for (auto &m : *meshsets_ptr) {
756 if (m->getMeshsetId() == med_mesh_name_id)
757 mesh_name = m->getName();
758 max_id = (max_id < m->getMeshsetId()) ? m->getMeshsetId() : max_id;
759 }
760
761 // Create a mesh
762 CHKERR MEDmeshCr(fid, mesh_name.c_str(), 3, 3, MED_UNSTRUCTURED_MESH,
763 "Mesh created", dtUnit, MED_SORT_DTIT, MED_CARTESIAN,
764 axisName, axisUnit);
765
766 // Create a map of meshset names to family ids and shared meshsets
767 med_int family_id = 0;
768 std::map<std::vector<std::string>, std::tuple<med_int, std::vector<int>>>
769 shared_meshsets_map;
770 std::map<EntityHandle, med_int> entityHandle_family_map;
771
772 // initialise zero family id
773 shared_meshsets_map[std::vector<string>()] =
774 std::make_tuple(family_id, std::vector<int>());
775
776 // function to get meshset names
777 auto get_set_name = [&](const CubitMeshSets *iit) {
778 std::string bc_type_name;
779 if (iit->getBcTypeULong() & BLOCKSET) {
780
781 bc_type_name = iit->getName();
782 if (bc_type_name == "NoNameSet") {
783 bc_type_name = "BLOCKSET_";
784 bc_type_name += std::to_string(iit->getMeshsetId());
785 }
786 return bc_type_name;
787 } else if (iit->getBcTypeULong() & SIDESET ||
788 iit->getBcTypeULong() & NODESET) {
789
790 CubitBCType cubitBcType(iit->getBcTypeULong());
791
792 unsigned jj = 0;
793 while (1 << jj != LASTSET_BC) {
794 const CubitBCType jj_bc_type = 1 << jj;
795 if ((iit->getBcType() & jj_bc_type).any()) {
796
797 std::string temp_bc_type_name = string(CubitBCNames[jj + 1]);
798 if (temp_bc_type_name.length() > 10) {
799 temp_bc_type_name = temp_bc_type_name.substr(0, 4);
800 }
801 bc_type_name += temp_bc_type_name;
802 bc_type_name += "_";
803 }
804 ++jj;
805 }
806 bc_type_name += std::to_string(iit->getMeshsetId());
807 } else {
808 bc_type_name = "UnknownSet";
809 }
810 return bc_type_name;
811 };
812
813 // loop over all entities in the write range
814 // FIXME: This is a very slow implementation
815 // We should use a binary search tree
816 // instead of looping each entity in
817 // write_range_ptr to find shared meshsets
818 for (auto &entity : *write_range_ptr) {
819 // check if entity is shared with other meshsets
820 std::vector<int> shared_meshsets;
821 std::vector<string> shared_names;
822
823 // function to add shared meshsets to vectors
824 auto add_shared_meshset = [&](auto &other_meshset) {
825 std::string other_name = get_set_name(other_meshset);
826 if (std::find(shared_names.begin(), shared_names.end(), other_name) ==
827 shared_names.end()) {
828 shared_meshsets.push_back(other_meshset->getMeshsetId());
829 shared_names.push_back(other_name);
830 }
831 };
832
833 // loop over all meshsets provided
834 for (auto &other_meshset : *meshsets_ptr) {
835 // skip if meshset is the global meshset given by user
836 if (med_mesh_name_id == other_meshset->getMeshsetId())
837 continue;
838
839 Range other_entities;
840 EntityHandle other_set = other_meshset->getMeshset();
841 moab.get_entities_by_handle(other_set, other_entities);
842 if (other_entities.empty())
843 continue;
844
845 bool is_in_meshset = moab.contains_entities(other_set, &entity, 1);
846 if (is_in_meshset) {
847 // add shared meshset id to list
848 add_shared_meshset(other_meshset);
849 }
850 }
851
852 // check if shared meshset is already in map
853 auto it = shared_meshsets_map.find(shared_names);
854 if (it == shared_meshsets_map.end()) {
855 // create new family id
856 family_id++;
857 // create and insert new tuple into the map
858 it = shared_meshsets_map
859 .insert(
860 {shared_names, std::make_tuple(family_id, shared_meshsets)})
861 .first;
862 }
863 // assign family id to entity
864 entityHandle_family_map[entity] = std::get<0>(it->second);
865 }
866
867 // loop to create families based on shared meshsets map
868 for (const auto &it : shared_meshsets_map) {
869 // create family
870 std::string family_name = "F_";
871 const auto &[family_id, shared_meshsets] = it.second;
872 family_name += std::to_string(family_id);
873
874 const auto &shared_meshset_names = it.first;
875 std::string group_name;
876 for (const auto &name : shared_meshset_names) {
877 // get meshset name
878 std::string meshset_name = name;
879 meshset_name.resize(MED_LNAME_SIZE, ' ');
880 group_name += meshset_name;
881 }
882
883 // create family
884 CHKERR MEDfamilyCr(fid, mesh_name.c_str(), family_name.c_str(), family_id,
885 shared_meshset_names.size(), group_name.c_str());
886
887 MOFEM_LOG("MEDWORLD", Sev::inform)
888 << "Creating family " << family_name << " with id " << family_id
889 << " and " << shared_meshset_names.size() << " groups " << std::endl;
890 }
891
892 // write nodes from range
893 Range verts;
894 verts = write_range_ptr->subset_by_type(MBVERTEX);
895
896 // Prepare arrays to hold the coordinates
897 std::vector<med_float> coord_med(3 * verts.size());
898 std::vector<med_int> fam;
899 std::vector<med_int> tags;
900 std::map<EntityHandle, med_int> entityHandle_tag_map;
901 med_int med_node_num = 1;
902
903 // For each vertex, get its coordinates
904 for (Range::iterator it = verts.begin(); it != verts.end(); ++it) {
905 double coords[3];
906 moab.get_coords(&(*it), 1, coords);
907 coord_med[3 * (it - verts.begin())] = coords[0];
908 coord_med[3 * (it - verts.begin()) + 1] = coords[1];
909 coord_med[3 * (it - verts.begin()) + 2] = coords[2];
910 fam.push_back(entityHandle_family_map[*it]);
911 tags.push_back(*it);
912 // map entityHandle to new tag (node) number
913 entityHandle_tag_map[*it] = med_node_num;
914 med_node_num++;
915 }
916 // Write the coordinates to the MED file
917 CHKERR MEDmeshNodeWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT,
918 MED_UNDEF_DT, MED_FULL_INTERLACE, verts.size(),
919 &coord_med[0], MED_FALSE, "", MED_TRUE, &tags[0],
920 MED_TRUE, &fam[0]);
921
922 // get last vertex id
923 double last_tag = tags.back();
924 // loop to write elements to families
925 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
926 // get all entities of type ent_type
927 Range entities;
928 moab.get_entities_by_type(0, ent_type, entities, true);
929 Range ents_to_write;
930 ents_to_write = intersect(entities, *write_range_ptr);
931
932 if (ents_to_write.empty())
933 continue;
934
935 // vectors to write
936 std::vector<med_int> tag_number;
937 std::vector<med_int> family_number;
938 std::vector<med_int> connectivity;
939
940 // loop over all entities
941 for (auto &entity : ents_to_write) {
942 if (ent_type != MBVERTEX) {
943 last_tag++;
944 // get family id for entity
945 family_number.push_back(entityHandle_family_map[entity]);
946 // get tag number for entity
947 tag_number.push_back(last_tag);
948 // get connectivity for entity
949 std::vector<EntityHandle> conn;
950 moab.get_connectivity(&entity, 1, conn);
951 // add connectivity to vector from map of entityHandle to new tag number
952 for (auto &c : conn) {
953 connectivity.push_back(entityHandle_tag_map[c]);
954 }
955 } else {
956 continue;
957 }
958 }
959
960 auto get_med_element_type = [](EntityType ent_type) {
961 med_geometrie_element type;
962 switch (ent_type) {
963 case MBHEX:
964 type = MED_HEXA8;
965 break;
966 case MBTET:
967 type = MED_TETRA4;
968 break;
969 case MBQUAD:
970 type = MED_QUAD4;
971 break;
972 case MBTRI:
973 type = MED_TRIA3;
974 break;
975 case MBEDGE:
976 type = MED_SEG2;
977 break;
978 case MBPRISM:
979 type = MED_PENTA6;
980 break;
981 default:
982 type = MED_POINT1;
983 break;
984 }
985 return type;
986 };
987
988 med_geometrie_element med_type = get_med_element_type(ent_type);
989 if (ent_type == MBENTITYSET) {
990 continue;
991 }
992 CHKERR MEDmeshElementWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT, 0.,
993 MED_CELL, med_type, MED_NODAL, MED_FULL_INTERLACE,
994 family_number.size(), &connectivity[0], MED_FALSE,
995 nullptr, MED_TRUE, &tag_number[0], MED_TRUE,
996 &family_number[0]);
997
998 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Writing %i elements of type %i (%s) ",
999 family_number.size(), med_type,
1000 moab::CN::EntityTypeName(ent_type));
1001 }
1002
1003 // Close the MED file
1004 CHKERR MEDfermer(fid);
1006}
1007
1008MoFEMErrorCode MedInterface::readFields(const std::string &file_name,
1009 const std::string &field_name,
1010 const bool load_series,
1011 const int only_step, int verb) {
1012
1013 Interface &m_field = cOre;
1015 med_idt fid = MEDfileOpen((char *)file_name.c_str(), MED_LECTURE);
1016 if (fid < 0) {
1017 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1018 "Unable to open file '%s'", file_name.c_str());
1019 }
1020
1021 med_int num_comp = MEDfieldnComponentByName(fid, field_name.c_str());
1022 if (num_comp <= 0) {
1023 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1024 "Could not get number of components for MED field");
1025 }
1026
1027 char meshName[MED_NAME_SIZE + 1];
1028 char dtUnit[MED_SNAME_SIZE + 1];
1029 std::vector<char> compName(num_comp * MED_SNAME_SIZE + 1);
1030 std::vector<char> compUnit(num_comp * MED_SNAME_SIZE + 1);
1031 med_int numSteps = 0;
1032 med_type_champ type;
1033 med_bool localMesh;
1034 if (MEDfieldInfoByName(fid, field_name.c_str(), meshName, &localMesh, &type,
1035 &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
1036 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1037 "Could not get MED field info");
1038 }
1039
1040 // Get meshset
1041 MeshsetsManager *meshsets_manager_ptr;
1042 CHKERR m_field.getInterface(meshsets_manager_ptr);
1043 const CubitMeshSets *cubit_meshset_ptr;
1044 CHKERR meshsets_manager_ptr->getCubitMeshsetPtr(meshName, &cubit_meshset_ptr);
1045 EntityHandle meshset = cubit_meshset_ptr->getMeshset();
1046
1047 // Paraview can only plot field which have 1, 3, or 9 components. If field has
1048 // more that 9 comonents, it is stored on MOAB mesh but not viable in
1049 // ParaView.
1050 int num_comp_msh = (num_comp <= 1) ? 1
1051 : (num_comp <= 3) ? 3
1052 : (num_comp <= 9) ? 9
1053 : num_comp;
1054
1055 // Create tag to store nodal or cell values read form med file. Note that tag
1056 // has prefix MED to avoid collision with other tags.
1057 Tag th;
1058 std::string tag_name = "MED_" + field_name;
1059 {
1060 std::vector<double> def_val(num_comp_msh, 0);
1061 CHKERR m_field.get_moab().tag_get_handle(
1062 tag_name.c_str(), num_comp_msh, MB_TYPE_DOUBLE, th,
1063 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val[0]);
1064 }
1065
1066 // Warning! The ordering of the elements in the last two lists is
1067 // important: it should match the ordering of the MSH element types
1068 // (when elements are saved without tags, this governs the order
1069 // with which we implicitly index them in GModel::readMED)
1070 const med_entity_type entType[] = {MED_NODE, MED_CELL, MED_NODE_ELEMENT};
1071 const med_geometrie_element eleType[] = {
1072 MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
1073 MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_QUAD9, MED_TETRA10,
1074 MED_HEXA27, MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
1075 // const int nodesPerEle[] = {0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 9, 10, 27, 1, 8,
1076 // 20, 15, 13};
1077
1078 std::vector<std::pair<int, int>> pairs;
1079 for (unsigned int i = 0; i < sizeof(entType) / sizeof(entType[0]); i++) {
1080 for (unsigned int j = 0; j < sizeof(eleType) / sizeof(eleType[0]); j++) {
1081 if ((!i && !j) || j) {
1082 med_int n = numSteps;
1083 if (n > 0) {
1084 pairs.push_back(std::pair<int, int>(i, j));
1085 numSteps = std::max(numSteps, n);
1086 }
1087 if (!i && !j)
1088 break; // MED_NOEUD does not care about eleType
1089 }
1090 }
1091 }
1092
1093 if (numSteps < 1 || pairs.empty()) {
1094 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1095 "Nothing to import from MED file");
1096 }
1097
1098 if (load_series) {
1099 SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1100 }
1101
1102 if (numSteps == 1 && only_step > 0) {
1103 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1104 "Field has only one step. Cannot load step with -med_time_step, "
1105 "remove from command line.");
1106 } else if (numSteps == only_step) {
1107 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1108 "Cannot load final step %d. -med_time_step should be %d.",
1109 only_step, only_step - 1);
1110 }
1111
1112 for (int step = (only_step == -1) ? 0 : only_step; step < numSteps; step++) {
1113 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Loading step %d of field %s", step + 1,
1114 field_name.c_str());
1115 if (!load_series && only_step != step)
1116 continue;
1117
1118 // FIXME: in MED3 we might want to loop over all profiles instead
1119 // of relying of the default one
1120
1121 // FIXME: MED3 allows to store multi-step meshes; we should
1122
1123 for (unsigned int pair = 0; pair < pairs.size(); pair++) {
1124
1125 // get step info
1126 med_entite_maillage ent = entType[pairs[pair].first];
1127 med_geometrie_element ele = eleType[pairs[pair].second];
1128 med_int numdt, numit, ngauss;
1129 med_float dt;
1130 if (MEDfieldComputingStepInfo(fid, field_name.c_str(), step + 1, &numdt,
1131 &numit, &dt) < 0) {
1132 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1133 "Could not read step info");
1134 }
1135
1136 char locName[MED_NAME_SIZE + 1], profileName[MED_NAME_SIZE + 1];
1137
1138 // get number of values in the field (numVal takes the number of
1139 // Gauss points or the number of nodes per element into account,
1140 // but not the number of components)
1141 med_int profileSize;
1142 med_int numVal = MEDfieldnValueWithProfile(
1143 fid, field_name.c_str(), numdt, numit, ent, ele, 1,
1144 MED_COMPACT_STMODE, profileName, &profileSize, locName, &ngauss);
1145 numVal *= ngauss;
1146
1147 if (numVal <= 0)
1148 continue;
1149
1150 // read field data
1151 std::vector<double> val(numVal * num_comp);
1152 if (MEDfieldValueWithProfileRd(fid, field_name.c_str(), numdt, numit, ent,
1153 ele, MED_COMPACT_STMODE, profileName,
1154 MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
1155 (unsigned char *)&val[0]) < 0) {
1156 SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1157 "Could not read field values");
1158 }
1159
1160 switch (ent) {
1161 case MED_CELL: {
1162 EntityType ent_type = MBMAXTYPE;
1163 switch (ele) {
1164 case MED_TETRA4:
1165 case MED_TETRA10:
1166 ent_type = MBTET;
1167 break;
1168 case MED_HEXA8:
1169 case MED_HEXA20:
1170 ent_type = MBHEX;
1171 break;
1172 default:
1173 MOFEM_LOG_C("MEDWORLD", Sev::warning,
1174 "Not yet implemented for this cell %d", ele);
1175 }
1176 if (ent_type != MBMAXTYPE) {
1177 if (ngauss == 1) {
1178 Range ents;
1179 CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
1180 ents, true);
1181 double e_vals[num_comp_msh];
1182 bzero(e_vals, sizeof(double) * num_comp_msh);
1183 std::vector<double>::iterator vit = val.begin();
1184 for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
1185 for (int ii = 0; ii != num_comp; ii++, vit++) {
1186 e_vals[ii] = *vit;
1187 }
1188 CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
1189 }
1190 } else {
1191 Range ents;
1192 CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
1193 ents, true);
1194 if (ents.size() * ngauss * num_comp != val.size()) {
1195 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1196 "data inconsistency");
1197 }
1198 // FIXME simply average gauss values, far from perfect need fix
1199 double e_vals[num_comp_msh];
1200 std::vector<double>::iterator vit = val.begin();
1201 for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
1202 bzero(e_vals, sizeof(double) * num_comp_msh);
1203 for (int gg = 0; gg != ngauss; gg++) {
1204 for (int ii = 0; ii != num_comp; ii++, vit++) {
1205 e_vals[ii] += *vit / ngauss;
1206 }
1207 }
1208 CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
1209 }
1210 }
1211 }
1212 } break;
1213 case MED_NODE:
1214 case MED_NODE_ELEMENT: {
1215 EntityType ent_type = MBVERTEX;
1216 Range ents;
1217 CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type, ents,
1218 true);
1219 double e_vals[num_comp_msh];
1220 bzero(e_vals, sizeof(double) * num_comp_msh);
1221 std::vector<double>::iterator vit = val.begin();
1222 for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
1223 for (int ii = 0; ii != num_comp; ii++, vit++) {
1224 e_vals[ii] = *vit;
1225 }
1226 CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
1227 }
1228 } break;
1229 default:
1230 MOFEM_LOG_C("MEDWORLD", Sev::inform, "Entity type %d not implemented",
1231 ent);
1232 }
1233 }
1234 }
1235
1237}
1238
1239std::ostream &operator<<(std::ostream &os,
1240 const MedInterface::FieldData &field_data) {
1241 os << "field name: " << field_data.fieldName;
1242 os << " mesh name: " << field_data.meshName;
1243 os << " local mesh: " << ((field_data.localMesh) ? "true" : "false");
1244 os << std::endl;
1245 os << " componentNames:";
1246 for (unsigned int ff = 0; ff != field_data.componentNames.size(); ff++) {
1247 os << " " << field_data.componentNames[ff];
1248 }
1249 os << std::endl;
1250 os << " componentUnits:";
1251 for (unsigned int ff = 0; ff != field_data.componentUnits.size(); ff++) {
1252 os << " " << field_data.componentUnits[ff];
1253 }
1254 os << std::endl;
1255 os << " dtUnit: " << field_data.dtUnit << endl;
1256 os << " number of steps: " << field_data.ncSteps;
1257 return os;
1258}
1259
1260 } // namespace MoFEM
1261
1262#endif // WITH_MED
#define MOFEM_LOG_C(channel, severity, format,...)
Med file interface interface.
constexpr double a
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ NODESET
@ SIDESET
@ LASTSET_BC
@ BLOCKSET
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
static const char *const CubitBCNames[]
Names of types of sets and boundary conditions.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
add cubit meshset
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
double dt
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'j', 3 > j
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
double FieldData
Field data type.
Definition Types.hpp:25
std::bitset< 32 > CubitBCType
Definition Types.hpp:52
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
std::ostream & operator<<(std::ostream &os, const EntitiesFieldData::EntData &e)
static std::vector< med_geometrie_element > moab2med_element_type(const EntityType type)
constexpr auto field_name
constexpr double g
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition Core.hpp:82
this struct keeps basic methods for moab meshset about material and boundary conditions
EntityHandle getMeshset() const
get bc meshset
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
std::vector< std::string > componentNames
std::vector< std::string > componentUnits
Interface for load MED files.
MoFEMErrorCode readMed(const string &file, int verb=1)
read MED file
MoFEMErrorCode writeMed(boost::shared_ptr< Range > range_ptr=nullptr, int verb=1)
write MED file
std::vector< EntityHandle > meshMeshsets
meshset for each mesh
MoFEMErrorCode readFamily(const string &file, const int index, const std::map< int, Range > &family_elem_mapint, std::map< string, Range > &group_elem_map, int verb=1)
read family and groups
MoFEMErrorCode makeBlockSets(const std::map< string, Range > &group_elem_map, int verb=1)
make from groups meshsets
std::map< std::string, FieldData > fieldNames
std::string medFileName
MED file name.
MoFEMErrorCode getMeshsets(boost::shared_ptr< std::vector< const CubitMeshSets * > > &meshsets_ptr)
MoFEMErrorCode getFileNameFromCommandLine(int verb=1)
Get MED file name from command line.
MedInterface(const MoFEM::Core &core)
std::vector< std::string > meshNames
list of meshes in MED file
MoFEM::Core & cOre
core database
MoFEMErrorCode readFields(const std::string &file_name, const std::string &field_name, const bool load_series=false, const int only_step=-1, int verb=1)
MoFEMErrorCode readMesh(const string &file, const int index, std::map< int, Range > &family_elem_map, int verb=1)
read mesh from MED file
PetscBool flgFile
true if file name given in command line
MoFEMErrorCode medGetFieldNames(const string &file, int verb=1)
Get field names in MED file.
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Interface for managing meshsets containing materials and boundary conditions.
CubitMeshSet_multiIndex & getMeshsetsMultindex()
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.