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  types.push_back(MED_POINT1);
244  default:
245  break;
246  }
247 
248  return types;
249 }
250 
251 MoFEMErrorCode MedInterface::readMesh(const string &file, const int index,
252  std::map<int, Range> &family_elem_map,
253  int verb) {
254 
255  Interface &m_field = cOre;
256  MoFEMFunctionBegin;
257 
258  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
259  if (fid < 0) {
260  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
261  "Unable to open file '%s'", file.c_str());
262  }
263  med_int v[3], vf[3];
264  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
265  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
266  if (verb > 1)
267  MOFEM_LOG_C("MEDWORLD", Sev::noisy,
268  "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
269  vf[1], vf[2], v[0], v[1], v[2]);
270 
271  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
272  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
273  "Cannot read MED file older than V2.2");
274  }
275 
276  char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
277  bzero(mesh_name, MED_NAME_SIZE + 1);
278  bzero(mesh_desc, MED_COMMENT_SIZE + 1);
279  med_int space_dim;
280  med_mesh_type mesh_type;
281  med_int mesh_dim, n_step;
282  char dt_unit[MED_SNAME_SIZE + 1];
283  char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
284  med_sorting_type sorting_type;
285  med_axis_type axis_type;
286  if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
287  mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
288  axis_name, axis_unit) < 0) {
289  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
290  "Unable to read mesh information");
291  }
292  if (verb > 0)
293  MOFEM_LOG_C("MEDWORLD", Sev::inform, "Reading mesh %s nsteps %d", mesh_name,
294  n_step);
295 
296  switch (axis_type) {
297  case MED_CARTESIAN:
298  break;
299  case MED_SPHERICAL:
300  case MED_CYLINDRICAL:
301  default:
302  SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
303  "Curvilinear coordinate system implemented");
304  }
305  if (space_dim < 2) {
306  SETERRQ1(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
307  "Not implemented for space dim %d", space_dim);
308  }
309 
310  EntityHandle mesh_meshset;
311  {
312  MeshsetsManager *meshsets_manager_ptr;
313  CHKERR m_field.getInterface(meshsets_manager_ptr);
314  int max_id = 0;
315  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
316  max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
317  }
318  max_id++;
319  CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id,
320  std::string(mesh_name));
321  CubitMeshSet_multiIndex::index<
322  Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
323  cit =
324  meshsets_manager_ptr->getMeshsetsMultindex()
325  .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
326  .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
327  max_id++;
328  mesh_meshset = cit->getMeshset();
329  meshMeshsets.push_back(mesh_meshset);
330  }
331 
332  med_bool change_of_coord, geo_transform;
333  med_int num_nodes = MEDmeshnEntity(
334  fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
335  MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
336  if (num_nodes < 0) {
337  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
338  "Could not read number of MED nodes");
339  }
340  if (num_nodes == 0) {
341  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
342  "No nodes in MED mesh");
343  }
344  if (verb > 0)
345  MOFEM_LOG_C("MEDWORLD", Sev::inform, "Read number of nodes %d", num_nodes);
346 
347  std::vector<med_float> coord_med(space_dim * num_nodes);
348  if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
349  MED_NO_INTERLACE, &coord_med[0]) < 0) {
350  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
351  "Could not read MED node coordinates");
352  }
353 
354  // Add vertices to moab
355  ReadUtilIface *iface;
356  vector<double *> arrays_coord;
357  EntityHandle startv;
358  CHKERR m_field.get_moab().query_interface(iface);
359  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
360  Range verts(startv, startv + num_nodes - 1);
361  std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
362  arrays_coord[0]);
363  std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
364  arrays_coord[1]);
365  if (space_dim == 2) {
366  std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
367  } else {
368  std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
369  arrays_coord[2]);
370  }
371  CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
372  family_elem_map.clear();
373 
374  // get family for vertices
375  {
376  std::vector<med_int> nodes_tags(num_nodes, 0);
377  if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
378  MED_NODE, MED_NO_GEOTYPE,
379  &nodes_tags[0]) < 0) {
380  nodes_tags.clear();
381  // SETERRQ(
382  // m_field.get_comm(),
383  // MOFEM_OPERATION_UNSUCCESSFUL,
384  // "No family number for elements: using 0 as default family number"
385  // );
386  }
387  for (int i = 0; i < num_nodes; i++) {
388  // cerr << verts[i] << " " /*<< ele_tags[j] << " "*/ << nodes_tags[i] <<
389  // endl;
390  family_elem_map[nodes_tags.empty() ? i : nodes_tags[i]].insert(verts[i]);
391  }
392  }
393 
394  // read elements (loop over all possible MSH element types)
395  for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
396 
397  auto types = moab2med_element_type(ent_type);
398 
399  for (auto type : types) {
400 
401  // get number of cells
402  med_bool change_of_coord;
403  med_bool geo_transform;
404  med_int num_ele = MEDmeshnEntity(
405  fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL, type,
406  MED_CONNECTIVITY, MED_NODAL, &change_of_coord, &geo_transform);
407 
408  // get connectivity
409  if (num_ele <= 0)
410  continue;
411 
412  int num_nod_per_ele = type % 100;
413  if (verb > 0)
414  MOFEM_LOG("MEDWORLD", Sev::inform)
415  << "Reading elements " << num_ele << " of type "
416  << moab::CN::EntityTypeName(ent_type) << " number of nodes "
417  << num_nod_per_ele;
418 
419  std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
420  if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
421  MED_CELL, type, MED_NODAL,
422  MED_FULL_INTERLACE, &conn_med[0]) < 0) {
423  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
424  "Could not read MED elements");
425  }
426 
427  // cerr << "type " << ent_type << " ";
428  // cerr << "num_ele " << num_ele << " " << num_nod_per_ele << endl;;
429 
430  Range ents;
431 
432  if (ent_type != MBVERTEX) {
433  EntityHandle *conn_moab;
434  EntityHandle starte;
435  CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
436  starte, conn_moab);
437  switch (ent_type) {
438  // FIXME: Some connectivity could not work, need to run and test
439  case MBTET: {
440  int ii = 0;
441  for (int ee = 0; ee != num_ele; ee++) {
442  EntityHandle n[4];
443  for (int nn = 0; nn != num_nod_per_ele; nn++) {
444  // conn_moab[ii] = verts[conn_med[ii]-1];
445  n[nn] = verts[conn_med[ii + nn] - 1];
446  }
447 
448  std::array<int, 10> nodes_map{
449 
450  1, 0, 2, 3,
451 
452  4, 6, 5, 8,
453  7, 9
454 
455  };
456 
457  for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
458  conn_moab[ii] = n[nodes_map[nn]];
459  }
460  }
461  } break;
462  default: {
463  int ii = 0;
464  for (int ee = 0; ee != num_ele; ee++) {
465  for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
466  // cerr << conn_med[ii] << " ";
467  conn_moab[ii] = verts[conn_med[ii] - 1];
468  }
469  // cerr << endl;
470  }
471  }
472  }
473  CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
474  conn_moab);
475  ents = Range(starte, starte + num_ele - 1);
476  } else {
477  // This is special case when in med vertices are defined as elements
478  int ii = 0;
479  std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
480  for (int ee = 0; ee != num_ele; ++ee)
481  for (int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
482  conn_moab[ii] = verts[conn_med[ii] - 1];
483  ents.insert_list(conn_moab.begin(), conn_moab.end());
484  }
485 
486  // Add elements to family meshset
487  CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
488 
489  // get family for cells
490  {
491  std::vector<med_int> fam(num_ele, 0);
492  if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
493  MED_CELL, type, &fam[0]) < 0) {
494  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
495  "No family number for elements: using 0 as default family "
496  "number");
497  }
498  // std::vector<med_int> ele_tags(num_ele);
499  // if(MEDmeshEntityNumberRd(
500  // fid, mesh_name,MED_NO_DT,MED_NO_IT,MED_CELL,type,&ele_tags[0]) < 0
501  // ) {
502  // ele_tags.clear();
503  // }
504  for (int j = 0; j < num_ele; j++) {
505  // cerr << ents[j] << " " /*<< ele_tags[j] << " "*/ << fam[j] << endl;
506  family_elem_map[fam[j]].insert(ents[j]);
507  }
508  }
509  }
510  }
511 
512  if (MEDfileClose(fid) < 0) {
513  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
514  "Unable to close file '%s'", file.c_str());
515  }
516 
517  MoFEMFunctionReturn(0);
518 }
519 
520 MoFEMErrorCode
521 MedInterface::readFamily(const string &file, const int index,
522  const std::map<int, Range> &family_elem_map,
523  std::map<string, Range> &group_elem_map, int verb) {
524  //
525  Interface &m_field = cOre;
526  MoFEMFunctionBegin;
527 
528  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
529  if (fid < 0) {
530  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
531  "Unable to open file '%s'", file.c_str());
532  }
533  med_int v[3], vf[3];
534  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
535  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
536 
537  if (verb > 1) {
538  MOFEM_LOG_C("MEDWORLD", Sev::noisy,
539  "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
540  vf[1], vf[2], v[0], v[1], v[2]);
541  }
542  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
543  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
544  "Cannot read MED file older than V2.2");
545  }
546 
547  // clear groups
548  group_elem_map.clear();
549 
550  med_int num_families = MEDnFamily(fid, meshNames[index].c_str());
551  if (num_families < 0) {
552  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
553  "Could not read MED families");
554  }
555  for (int i = 0; i < num_families; i++) {
556  med_int num_attrib =
557  (vf[0] == 2)
558  ? MEDnFamily23Attribute(fid, meshNames[index].c_str(), i + 1)
559  : 0;
560  med_int num_groups = MEDnFamilyGroup(fid, meshNames[index].c_str(), i + 1);
561  if (num_attrib < 0 || num_groups < 0) {
562  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
563  "Could not read MED groups or attributes");
564  }
565 
566  std::vector<med_int> attribId(num_attrib + 1);
567  std::vector<med_int> attrib_val(num_attrib + 1);
568  std::vector<char> attrib_des(MED_COMMENT_SIZE * num_attrib + 1);
569  std::vector<char> group_names(MED_LNAME_SIZE * num_groups + 1);
570  char family_name[MED_NAME_SIZE + 1];
571  med_int family_num;
572 
573  if (vf[0] == 2) { // MED2 file
574  if (MEDfamily23Info(fid, meshNames[index].c_str(), i + 1, family_name,
575  &attribId[0], &attrib_val[0], &attrib_des[0],
576  &family_num, &group_names[0]) < 0) {
577  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
578  "Could not read info for MED2 family %d", i + 1);
579  }
580  } else {
581  if (MEDfamilyInfo(fid, meshNames[index].c_str(), i + 1, family_name,
582  &family_num, &group_names[0]) < 0) {
583  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
584  "Could not read info for MED3 family %d", i + 1);
585  }
586  }
587 
588  // cerr << family_name << " " << family_num << " " << num_groups << endl;
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()) {
594  MOFEM_LOG_C(
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  // cerr << string(&group_names[MED_LNAME_SIZE*g]) << endl;
603  }
604  }
605  }
606 
607  if (MEDfileClose(fid) < 0) {
608  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
609  "Unable to close file '%s'", file.c_str());
610  }
611 
612  MoFEMFunctionReturn(0);
613 }
614 
615 MoFEMErrorCode
616 MedInterface::makeBlockSets(const std::map<string, Range> &group_elem_map,
617  int verb) {
618 
619  Interface &m_field = cOre;
620  MoFEMFunctionBegin;
621  MeshsetsManager *meshsets_manager_ptr;
622  CHKERR m_field.getInterface(meshsets_manager_ptr);
623 
624  int max_id = 0;
625  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit)) {
626  max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
627  }
628  max_id++;
629 
630  for (std::map<string, Range>::const_iterator git = group_elem_map.begin();
631  git != group_elem_map.end(); git++) {
632  CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id, git->first);
633  CubitMeshSet_multiIndex::index<
634  Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
635  cit =
636  meshsets_manager_ptr->getMeshsetsMultindex()
637  .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
638  .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
639  EntityHandle meshsets = cit->getMeshset();
640  if (!git->second.empty()) {
641  CHKERR m_field.get_moab().add_entities(meshsets, git->second);
642  }
643  max_id++;
644  }
645 
646  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, cit))
647  MOFEM_LOG("MEDWORLD", Sev::verbose) << *cit;
648 
649  MoFEMFunctionReturn(0);
650 }
651 
652 MoFEMErrorCode MedInterface::readMed(int verb) {
653  MoFEMFunctionBegin;
654  if (medFileName.empty()) {
655  CHKERR getFileNameFromCommandLine(verb);
656  }
657  CHKERR readMed(medFileName, verb);
658  MoFEMFunctionReturn(0);
659 }
660 
661 MoFEMErrorCode MedInterface::getMeshsets(
662  boost::shared_ptr<std::vector<const CubitMeshSets *>> &meshsets_ptr) {
663  Interface &m_field = cOre;
664  MoFEMFunctionBegin;
665 
666  auto &meshsets_idx =
667  m_field.getInterface<MeshsetsManager>()->getMeshsetsMultindex();
668  std::vector<const CubitMeshSets *> meshsets;
669 
670  for (auto &m : meshsets_idx) {
671  meshsets.push_back(&m);
672  }
673  // Sort meshsets based on meshsetId
674  std::sort(meshsets.begin(), meshsets.end(),
675  [](const CubitMeshSets *a, const CubitMeshSets *b) {
676  return a->getMeshsetId() < b->getMeshsetId();
677  });
678 
679  meshsets_ptr =
680  boost::make_shared<std::vector<const CubitMeshSets *>>(meshsets);
681  MoFEMFunctionReturn(0);
682 }
683 
684 MoFEMErrorCode MedInterface::writeMed(boost::shared_ptr<Range> range_ptr,
685  int verb) {
686  MoFEMFunctionBegin;
687  if (medFileName.empty()) {
688  CHKERR getFileNameFromCommandLine(verb);
689  }
690 
691  // Get all meshsets
692  boost::shared_ptr<std::vector<const CubitMeshSets *>> meshsets_ptr;
693  getMeshsets(meshsets_ptr);
694 
695  CHKERR writeMed(medFileName, meshsets_ptr, range_ptr, verb);
696  MoFEMFunctionReturn(0);
697 }
698 
699 MoFEMErrorCode MedInterface::writeMed(
700  const string &file,
701  boost::shared_ptr<std::vector<const CubitMeshSets *>> meshsets_ptr,
702  boost::shared_ptr<Range> write_range_ptr, int verb) {
703  Interface &m_field = cOre;
704  // Get the MOAB instance from the field
705  moab::Interface &moab = m_field.get_moab();
706 
707  MoFEMFunctionBeginHot;
708  MOFEM_LOG("MEDWORLD", Sev::warning)
709  << "WRITE_MED IS EXPERIMENTAL, MAY CONTAIN BUGS, ALWAYS CHECK THE OUTPUT "
710  "FILE";
711  // Open a med file with the specified version
712  med_idt fid =
713  MEDfileVersionOpen((char *)file.c_str(), MED_ACC_CREAT, MED_MAJOR_NUM,
714  MED_MINOR_NUM, MED_RELEASE_NUM);
715  // Throw an error if the file could not be opened
716  if (fid < 0) {
717  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
718  "Cannot create MED file");
719  return 0;
720  }
721 
722  // Write appropriate header
723  CHKERR MEDfichDesEcr(fid, (char *)"MED file generated by MoFEM");
724 
725  // Generate empty mesh
726  char dtUnit[MED_SNAME_SIZE + 1] = ""; // time unit
727  char axisName[3 * MED_SNAME_SIZE + 1] = ""; // axis name
728  char axisUnit[3 * MED_SNAME_SIZE + 1] = ""; // axis unit
729 
730  PetscBool is_cubit_meshset = PETSC_TRUE;
731  int med_mesh_name_id = 0;
732 
733  // Create the mesh
734  char mesh_name_char[255];
735  std::string mesh_name = "Mesh";
736 
737  // Get mesh_name from command line
738  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "MED mesh options", "");
739  CHKERR PetscOptionsString("-med_mesh_name", "get med mesh name", "",
740  mesh_name.c_str(), mesh_name_char, 255, PETSC_NULL);
741  ierr = PetscOptionsEnd();
742  CHKERRG(ierr);
743 
744  mesh_name = mesh_name_char;
745 
746  // check if mesh_name is a meshset
747  for (auto &m : *meshsets_ptr) {
748  if (m->getName() == mesh_name) {
749  med_mesh_name_id = m->getMeshsetId();
750  is_cubit_meshset = PETSC_FALSE;
751  break;
752  }
753  }
754 
755  // Get the maximum meshset id
756  int max_id = 0;
757  for (auto &m : *meshsets_ptr) {
758  if (m->getMeshsetId() == med_mesh_name_id)
759  mesh_name = m->getName();
760  max_id = (max_id < m->getMeshsetId()) ? m->getMeshsetId() : max_id;
761  }
762 
763  // Create a mesh
764  CHKERR MEDmeshCr(fid, mesh_name.c_str(), 3, 3, MED_UNSTRUCTURED_MESH,
765  "Mesh created", dtUnit, MED_SORT_DTIT, MED_CARTESIAN,
766  axisName, axisUnit);
767 
768  // Create a map of meshset names to family ids and shared meshsets
769  med_int family_id = 0;
770  std::map<std::vector<std::string>, std::tuple<med_int, std::vector<int>>>
771  shared_meshsets_map;
772  std::map<EntityHandle, med_int> entityHandle_family_map;
773 
774  // initialise zero family id
775  shared_meshsets_map[std::vector<string>()] =
776  std::make_tuple(family_id, std::vector<int>());
777 
778  // function to get meshset names
779  auto get_set_name = [&](const CubitMeshSets *iit) {
780  std::string bc_type_name;
781  if (iit->getBcTypeULong() & BLOCKSET) {
782 
783  bc_type_name = iit->getName();
784  if (bc_type_name == "NoNameSet") {
785  bc_type_name = "BLOCKSET_NoNameSet_";
786  bc_type_name += std::to_string(iit->getMeshsetId());
787  }
788  return bc_type_name;
789  } else if (iit->getBcTypeULong() & SIDESET ||
790  iit->getBcTypeULong() & NODESET) {
791 
792  CubitBCType cubitBcType(iit->getBcTypeULong());
793 
794  unsigned jj = 0;
795  while (1 << jj != LASTSET_BC) {
796  const CubitBCType jj_bc_type = 1 << jj;
797  if ((iit->getBcType() & jj_bc_type).any()) {
798 
799  bc_type_name += string(CubitBCNames[jj + 1]);
800  bc_type_name += "_";
801  }
802  ++jj;
803  }
804  bc_type_name += std::to_string(iit->getMeshsetId());
805  } else {
806  bc_type_name = "UnknownSet";
807  }
808  return bc_type_name;
809  };
810 
811  // loop over all entities in the write range
812  // FIXME: This is a very slow implementation
813  // We should use a binary search tree
814  // instead of looping each entity in
815  // write_range_ptr to find shared meshsets
816  for (auto &entity : *write_range_ptr) {
817  // check if entity is shared with other meshsets
818  std::vector<int> shared_meshsets;
819  std::vector<string> shared_names;
820 
821  // function to add shared meshsets to vectors
822  auto add_shared_meshset = [&](auto &other_meshset) {
823  std::string other_name = get_set_name(other_meshset);
824  if (std::find(shared_names.begin(), shared_names.end(), other_name) ==
825  shared_names.end()) {
826  shared_meshsets.push_back(other_meshset->getMeshsetId());
827  shared_names.push_back(other_name);
828  }
829  };
830 
831  // loop over all meshsets provided
832  for (auto &other_meshset : *meshsets_ptr) {
833  // skip if meshset is the global meshset given by user
834  if (med_mesh_name_id == other_meshset->getMeshsetId())
835  continue;
836 
837  Range other_entities;
838  EntityHandle other_set = other_meshset->getMeshset();
839  CHKERR moab.get_entities_by_handle(other_set, other_entities);
840 
841  // get entity type
842  EntityType ent_type = moab.type_from_handle(entity);
843 
844  bool is_in_meshset = moab.contains_entities(other_set, &entity, 1);
845  if (is_in_meshset) {
846  if (ent_type == MBVERTEX) {
847  // add shared meshset id to list
848  // check if higher dimension entity is in meshset i.e not a nodeset
849  bool is_in_higher_dim = false;
850  Range entities_in_higher_dim;
851  for (int i = 1; i < 4; i++) {
852  CHKERR moab.get_entities_by_dimension(other_set, i,
853  entities_in_higher_dim);
854  if (!entities_in_higher_dim.empty()) {
855  is_in_higher_dim = true;
856  break;
857  }
858  }
859  if (is_in_higher_dim) {
860  continue;
861  }
862  // check if name is already added to shared_names
863  add_shared_meshset(other_meshset);
864  } else {
865  // add shared meshset id to list
866  add_shared_meshset(other_meshset);
867  }
868  }
869  }
870 
871  // check if shared meshset is already in map
872  auto it = shared_meshsets_map.find(shared_names);
873  if (it == shared_meshsets_map.end()) {
874  // create new family id
875  family_id++;
876  // create and insert new tuple into the map
877  it = shared_meshsets_map
878  .insert(
879  {shared_names, std::make_tuple(family_id, shared_meshsets)})
880  .first;
881  }
882  // assign family id to entity
883  entityHandle_family_map[entity] = std::get<0>(it->second);
884  }
885 
886  // loop to create families based on shared meshsets map
887  for (const auto &it : shared_meshsets_map) {
888  // create family
889  std::string family_name = "F_";
890  const auto &[family_id, shared_meshsets] = it.second;
891  family_name += std::to_string(family_id);
892 
893  const auto &shared_meshset_names = it.first;
894  std::string group_name;
895  for (const auto &name : shared_meshset_names) {
896  // get meshset name
897  std::string meshset_name = name;
898  meshset_name.resize(MED_LNAME_SIZE, ' ');
899  group_name += meshset_name;
900  }
901 
902  // create family
903  CHKERR MEDfamilyCr(fid, mesh_name.c_str(), family_name.c_str(), family_id,
904  shared_meshset_names.size(), group_name.c_str());
905 
906  MOFEM_LOG("MEDWORLD", Sev::inform)
907  << "Creating family " << family_name << " with id " << family_id
908  << " and " << shared_meshset_names.size() << " groups " << std::endl;
909  }
910 
911  // write nodes
912  Range verts;
913  moab.get_entities_by_type(0, MBVERTEX, verts);
914  // Prepare arrays to hold the coordinates
915  std::vector<med_float> coord_med(3 * verts.size());
916  std::vector<med_int> fam;
917  std::vector<med_int> tags;
918 
919  // For each vertex, get its coordinates
920  for (Range::iterator it = verts.begin(); it != verts.end(); ++it) {
921  double coords[3];
922  moab.get_coords(&(*it), 1, coords);
923  coord_med[3 * (it - verts.begin())] = coords[0];
924  coord_med[3 * (it - verts.begin()) + 1] = coords[1];
925  coord_med[3 * (it - verts.begin()) + 2] = coords[2];
926  fam.push_back(entityHandle_family_map[*it]);
927  tags.push_back(*it);
928  }
929  // Write the coordinates to the MED file
930  CHKERR MEDmeshNodeWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT,
931  MED_UNDEF_DT, MED_FULL_INTERLACE, verts.size(),
932  &coord_med[0], MED_FALSE, "", MED_TRUE, &tags[0],
933  MED_TRUE, &fam[0]);
934 
935  // loop to write elements to families
936  for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
937  // get all entities of type ent_type
938  Range entities;
939  moab.get_entities_by_type(0, ent_type, entities, true);
940  Range ents_to_write;
941  ents_to_write = intersect(entities, *write_range_ptr);
942 
943  if (ents_to_write.empty())
944  continue;
945 
946  // vectors to write
947  std::vector<med_int> tag_number;
948  std::vector<med_int> family_number;
949  std::vector<med_int> connectivity;
950 
951  // loop over all entities
952  for (auto &entity : ents_to_write) {
953  if (ent_type != MBVERTEX) {
954  // get family id for entity
955  family_number.push_back(entityHandle_family_map[entity]);
956  // get tag number for entity
957  tag_number.push_back(entity);
958  // get connectivity for entity
959  std::vector<EntityHandle> conn;
960  moab.get_connectivity(&entity, 1, conn);
961  for (auto &c : conn) {
962  connectivity.push_back(c);
963  }
964  } else {
965  continue;
966  }
967  }
968 
969  auto get_med_element_type = [](EntityType ent_type) {
970  med_geometrie_element type;
971  switch (ent_type) {
972  case MBHEX:
973  type = MED_HEXA8;
974  break;
975  case MBTET:
976  type = MED_TETRA4;
977  break;
978  case MBQUAD:
979  type = MED_QUAD4;
980  break;
981  case MBTRI:
982  type = MED_TRIA3;
983  break;
984  case MBEDGE:
985  type = MED_SEG2;
986  break;
987  case MBPRISM:
988  type = MED_PENTA6;
989  break;
990  default:
991  type = MED_POINT1;
992  break;
993  }
994  return type;
995  };
996 
997  med_geometrie_element med_type = get_med_element_type(ent_type);
998  if (ent_type == MBENTITYSET) {
999  continue;
1000  }
1001  CHKERR MEDmeshElementWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT, 0.,
1002  MED_CELL, med_type, MED_NODAL, MED_FULL_INTERLACE,
1003  family_number.size(), &connectivity[0], MED_FALSE,
1004  nullptr, MED_TRUE, &tag_number[0], MED_TRUE,
1005  &family_number[0]);
1006 
1007  MOFEM_LOG_C("MEDWORLD", Sev::inform, "Writing %i elements of type %i (%s) ",
1008  family_number.size(), med_type,
1009  moab::CN::EntityTypeName(ent_type));
1010  }
1011 
1012  // Close the MED file
1013  CHKERR MEDfermer(fid);
1014  MoFEMFunctionReturnHot(0);
1015 }
1016 
1017 MoFEMErrorCode MedInterface::readFields(const std::string &file_name,
1018  const std::string &field_name,
1019  const bool load_series,
1020  const int only_step, int verb) {
1021 
1022  Interface &m_field = cOre;
1023  MoFEMFunctionBeginHot;
1024  med_idt fid = MEDfileOpen((char *)file_name.c_str(), MED_LECTURE);
1025  if (fid < 0) {
1026  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1027  "Unable to open file '%s'", file_name.c_str());
1028  }
1029 
1030  med_int num_comp = MEDfieldnComponentByName(fid, field_name.c_str());
1031  if (num_comp <= 0) {
1032  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1033  "Could not get number of components for MED field");
1034  }
1035 
1036  char meshName[MED_NAME_SIZE + 1];
1037  char dtUnit[MED_SNAME_SIZE + 1];
1038  std::vector<char> compName(num_comp * MED_SNAME_SIZE + 1);
1039  std::vector<char> compUnit(num_comp * MED_SNAME_SIZE + 1);
1040  med_int numSteps = 0;
1041  med_type_champ type;
1042  med_bool localMesh;
1043  if (MEDfieldInfoByName(fid, field_name.c_str(), meshName, &localMesh, &type,
1044  &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
1045  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1046  "Could not get MED field info");
1047  }
1048 
1049  // Get meshset
1050  MeshsetsManager *meshsets_manager_ptr;
1051  CHKERR m_field.getInterface(meshsets_manager_ptr);
1052  const CubitMeshSets *cubit_meshset_ptr;
1053  CHKERR meshsets_manager_ptr->getCubitMeshsetPtr(meshName, &cubit_meshset_ptr);
1054  EntityHandle meshset = cubit_meshset_ptr->getMeshset();
1055 
1056  // Paraview can only plot field which have 1, 3, or 9 components. If field has
1057  // more that 9 comonents, it is stored on MOAB mesh but not viable in
1058  // ParaView.
1059  int num_comp_msh = (num_comp <= 1) ? 1
1060  : (num_comp <= 3) ? 3
1061  : (num_comp <= 9) ? 9
1062  : num_comp;
1063 
1064  // Create tag to store nodal or cell values read form med file. Note that tag
1065  // has prefix MED to avoid collision with other tags.
1066  Tag th;
1067  std::string tag_name = "MED_" + field_name;
1068  {
1069  std::vector<double> def_val(num_comp_msh, 0);
1070  CHKERR m_field.get_moab().tag_get_handle(
1071  tag_name.c_str(), num_comp_msh, MB_TYPE_DOUBLE, th,
1072  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val[0]);
1073  }
1074 
1075  // Warning! The ordering of the elements in the last two lists is
1076  // important: it should match the ordering of the MSH element types
1077  // (when elements are saved without tags, this governs the order
1078  // with which we implicitly index them in GModel::readMED)
1079  const med_entity_type entType[] = {MED_NODE, MED_CELL, MED_NODE_ELEMENT};
1080  const med_geometrie_element eleType[] = {
1081  MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
1082  MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_QUAD9, MED_TETRA10,
1083  MED_HEXA27, MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
1084  // const int nodesPerEle[] = {0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 9, 10, 27, 1, 8,
1085  // 20, 15, 13};
1086 
1087  std::vector<std::pair<int, int>> pairs;
1088  for (unsigned int i = 0; i < sizeof(entType) / sizeof(entType[0]); i++) {
1089  for (unsigned int j = 0; j < sizeof(eleType) / sizeof(eleType[0]); j++) {
1090  if ((!i && !j) || j) {
1091  med_int n = numSteps;
1092  if (n > 0) {
1093  pairs.push_back(std::pair<int, int>(i, j));
1094  numSteps = std::max(numSteps, n);
1095  }
1096  if (!i && !j)
1097  break; // MED_NOEUD does not care about eleType
1098  }
1099  }
1100  }
1101 
1102  if (numSteps < 1 || pairs.empty()) {
1103  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1104  "Nothing to import from MED file");
1105  }
1106 
1107  if (load_series) {
1108  SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1109  }
1110 
1111  for (int step = (only_step == -1) ? 0 : only_step; step < numSteps; step++) {
1112 
1113  if (!load_series && only_step != step)
1114  continue;
1115 
1116  // cerr << only_step << " " << step << endl;
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  // int mult = 1;
1151  // if(ent == MED_NODE_ELEMENT) {
1152  // mult = nodesPerEle[pairs[pair].second];
1153  //}
1154  // else if(ngauss != 1){
1155  // mult = ngauss;
1156  //}
1157 
1158  // read field data
1159  std::vector<double> val(numVal * num_comp);
1160  if (MEDfieldValueWithProfileRd(fid, field_name.c_str(), numdt, numit, ent,
1161  ele, MED_COMPACT_STMODE, profileName,
1162  MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
1163  (unsigned char *)&val[0]) < 0) {
1164  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
1165  "Could not read field values");
1166  }
1167 
1168  // if (verb > 2) {
1169  // // FIXME: This not looks ok for me
1170  // cerr << ent << " " << ele << endl;
1171  // cerr << string(meshName) << " : " << string(profileName) << " : "
1172  // << string(locName) << " : " << profileSize << " : " << ngauss
1173  // << endl;
1174  // }
1175 
1176  switch (ent) {
1177  case MED_CELL: {
1178  EntityType ent_type = MBMAXTYPE;
1179  switch (ele) {
1180  case MED_TETRA4:
1181  case MED_TETRA10:
1182  ent_type = MBTET;
1183  break;
1184  case MED_HEXA8:
1185  ent_type = MBHEX;
1186  break;
1187  default:
1188  MOFEM_LOG_C("MEDWORLD", Sev::warning,
1189  "Not yet implemented for this cell %d", ele);
1190  }
1191  if (ent_type != MBMAXTYPE) {
1192  if (ngauss == 1) {
1193  Range ents;
1194  CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
1195  ents, true);
1196  double e_vals[num_comp_msh];
1197  bzero(e_vals, sizeof(double) * num_comp_msh);
1198  std::vector<double>::iterator vit = val.begin();
1199  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
1200  for (int ii = 0; ii != num_comp; ii++, vit++) {
1201  e_vals[ii] = *vit;
1202  }
1203  CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
1204  }
1205  } else {
1206  Range ents;
1207  CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
1208  ents, true);
1209  if (ents.size() * ngauss * num_comp != val.size()) {
1210  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1211  "data inconsistency");
1212  }
1213  // FIXME simply average gauss values, far from perfect need fix
1214  double e_vals[num_comp_msh];
1215  std::vector<double>::iterator vit = val.begin();
1216  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
1217  bzero(e_vals, sizeof(double) * num_comp_msh);
1218  for (int gg = 0; gg != ngauss; gg++) {
1219  for (int ii = 0; ii != num_comp; ii++, vit++) {
1220  e_vals[ii] += *vit / ngauss;
1221  }
1222  }
1223  CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
1224  }
1225  }
1226  }
1227  // SETERRQ1(
1228  // m_field.get_comm(),
1229  // MOFEM_NOT_IMPLEMENTED,
1230  // "Not implemented for more gauss pts ngauss = %d",
1231  // ngauss
1232  // );
1233  } break;
1234  case MED_NODE:
1235  case MED_NODE_ELEMENT: {
1236  EntityType ent_type = MBVERTEX;
1237  Range ents;
1238  CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type, ents,
1239  true);
1240  double e_vals[num_comp_msh];
1241  bzero(e_vals, sizeof(double) * num_comp_msh);
1242  std::vector<double>::iterator vit = val.begin();
1243  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
1244  for (int ii = 0; ii != num_comp; ii++, vit++) {
1245  e_vals[ii] = *vit;
1246  }
1247  CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
1248  }
1249  } break;
1250  default:
1251  MOFEM_LOG_C("MEDWORLD", Sev::inform, "Entity type %d not implemented",
1252  ent);
1253  }
1254  }
1255  }
1256 
1257  MoFEMFunctionReturnHot(0);
1258 }
1259 
1260 std::ostream &operator<<(std::ostream &os,
1261  const MedInterface::FieldData &field_data) {
1262  os << "field name: " << field_data.fieldName;
1263  os << " mesh name: " << field_data.meshName;
1264  os << " local mesh: " << ((field_data.localMesh) ? "true" : "false");
1265  os << std::endl;
1266  // os << " field type: ";
1267  // switch (field_data.fieldType) {
1268  // case MED_FLOAT64: os << "MED_FLOAT64"; break;
1269  // case MED_INT32: os << "MED_INT32"; break;
1270  // case MED_INT64: os << "MED_INT64"; break;
1271  // case MED_INT: os << "MED_INT"; break;
1272  // };
1273  os << " componentNames:";
1274  for (unsigned int ff = 0; ff != field_data.componentNames.size(); ff++) {
1275  os << " " << field_data.componentNames[ff];
1276  }
1277  os << std::endl;
1278  os << " componentUnits:";
1279  for (unsigned int ff = 0; ff != field_data.componentUnits.size(); ff++) {
1280  os << " " << field_data.componentUnits[ff];
1281  }
1282  os << std::endl;
1283  os << " dtUnit: " << field_data.dtUnit << endl;
1284  os << " number of steps: " << field_data.ncSteps;
1285  return os;
1286 }
1287 
1288  } // namespace MoFEM
1289 
1290 #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:661
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:616
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:521
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:684
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:1017
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:251
_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