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