v0.9.2
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 /*
9  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12  * License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
16  */
17 
18 #ifdef WITH_MED
19 
20 extern "C" {
21 #include <med.h>
22 }
23 
24 #if (MED_MAJOR_NUM >= 3)
25 // To avoid too many ifdefs below we use defines for the bits of the
26 // API that did not change too much between MED2 and MED3. If we
27 // remove MED2 support at some point, please remove these defines and
28 // replace the symbols accordingly.
29 #define med_geometrie_element med_geometry_type
30 #define med_entite_maillage med_entity_type
31 #define med_type_champ med_field_type
32 #define MED_LECTURE MED_ACC_RDONLY
33 #define MED_LECTURE_AJOUT MED_ACC_RDEXT
34 #define MED_NOEUD MED_NODE
35 #define MED_MAILLE MED_CELL
36 #define MED_NOEUD_MAILLE MED_NODE_ELEMENT
37 #define MED_NOPFL MED_NO_PROFILE
38 #define MEDouvrir MEDfileOpen
39 #define MEDfermer MEDfileClose
40 #define MEDnChamp MEDfieldnComponent
41 #define MEDnValProfil MEDprofileSizeByName
42 #else
43 #error "MED file is not MED_MAJOR_NUM == 3"
44 #endif
45 
46 #include <MedInterface.hpp>
47 
48 #define MedFunctionBegin \
49  MoFEMFunctionBegin; \
50  MOFEM_LOG_CHANNEL("WORLD"); \
51  MOFEM_LOG_CHANNEL("SYNC"); \
52  MOFEM_LOG_FUNCTION(); \
53  MOFEM_LOG_TAG("SYNC", "MedInterface"); \
54  MOFEM_LOG_TAG("WORLD", "MedInterface")
55 
56 namespace MoFEM {
57 
59  UnknownInterface **iface) const {
61  *iface = NULL;
62  if (uuid == IDD_MOFEMMedInterface) {
63  *iface = const_cast<MedInterface *>(this);
65  }
66  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
68 }
69 
71  : cOre(const_cast<Core &>(core)), flgFile(PETSC_FALSE) {}
72 
74  Interface &m_field = cOre;
75  char mesh_file_name[255];
77  ierr = PetscOptionsBegin(m_field.get_comm(), "", "MED Interface", "none");
78  CHKERRG(ierr);
79  ierr = PetscOptionsString("-med_file", "med file name", "", "mesh.med",
80  mesh_file_name, 255, &flgFile);
81  CHKERRG(ierr);
82  ierr = PetscOptionsEnd();
83  CHKERRG(ierr);
84  if (flgFile == PETSC_TRUE) {
85  medFileName = std::string(mesh_file_name);
86  } else {
87  medFileName = std::string("mesh.med");
88  }
90 }
91 
92 MoFEMErrorCode MedInterface::medGetFieldNames(const string &file, int verb) {
93  Interface &m_field = cOre;
95  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
96  if (fid < 0) {
97  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
98  "Unable to open file '%s'", file.c_str());
99  }
100  med_int num_fields = MEDnField(fid);
101  for (int index = 0; index < num_fields; index++) {
102 
103  med_int num_comp = MEDfieldnComponent(fid, index + 1);
104  if (num_comp <= 0) {
105  SETERRQ(m_field.get_comm(), MOFEM_IMPOSIBLE_CASE,
106  "Could not get number of components for MED field");
107  }
108 
109  char name[MED_NAME_SIZE + 1], mesh_name[MED_NAME_SIZE + 1];
110  char dt_unit[MED_SNAME_SIZE + 1];
111  std::vector<char> comp_name(num_comp * MED_SNAME_SIZE + 1);
112  std::vector<char> comp_unit(num_comp * MED_SNAME_SIZE + 1);
113  med_int num_steps = 0;
114  med_type_champ type;
115  med_bool local_mesh;
116  if (MEDfieldInfo(fid, index + 1, name, mesh_name, &local_mesh, &type,
117  &comp_name[0], &comp_unit[0], dt_unit, &num_steps) < 0) {
118  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
119  "Could not get MED field info");
120  }
121 
122  std::string field_name = std::string(name);
123  fieldNames[field_name] = FieldData();
124  fieldNames[field_name].fieldName = field_name;
125  fieldNames[field_name].meshName = std::string(mesh_name);
126  fieldNames[field_name].localMesh = (local_mesh == MED_TRUE);
127  for (int ff = 0; ff != num_comp; ff++) {
128  fieldNames[field_name].componentNames.push_back(
129  std::string(&comp_name[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
130  fieldNames[field_name].componentUnits.push_back(
131  std::string(&comp_unit[ff * MED_SNAME_SIZE], MED_SNAME_SIZE));
132  }
133  fieldNames[field_name].dtUnit = std::string(&dt_unit[0]);
134  fieldNames[field_name].ncSteps = num_steps;
135 
136  if (verb > 0)
137  MOFEM_LOG("WORLD", Sev::inform) << fieldNames[name];
138 
139  }
140  if (MEDfileClose(fid) < 0) {
141  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
142  "Unable to close file '%s'", file.c_str());
143  }
145 }
146 
149  if (medFileName.empty()) {
151  }
154 }
155 
156 MoFEMErrorCode MedInterface::readMed(const string &file, int verb) {
157  Interface &m_field = cOre;
159 
160  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
161  if (fid < 0) {
162  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
163  "Unable to open file '%s'", file.c_str());
164  }
165 
166  med_int v[3], vf[3];
167  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
168  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
169 
170  if (verb > 0)
171  MOFEM_LOG_C("WORLD", Sev::inform,
172  "Reading MED file V%d.%d.%d using MED library V%d.%d.%d",
173  vf[0], vf[1], vf[2], v[0], v[1], v[2]);
174 
175  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
176  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
177  "Cannot read MED file older than V2.2");
178  }
179 
180  for (int i = 0; i < MEDnMesh(fid); i++) {
181  char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
182  bzero(mesh_name, MED_NAME_SIZE);
183  bzero(mesh_desc, MED_COMMENT_SIZE);
184  med_int space_dim;
185  med_mesh_type mesh_type;
186  med_int mesh_dim, n_step;
187  char dt_unit[MED_SNAME_SIZE + 1];
188  char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
189  med_sorting_type sorting_type;
190  med_axis_type axis_type;
191  if (MEDmeshInfo(fid, i + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
192  mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
193  axis_name, axis_unit) < 0) {
194  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
195  "Unable to read mesh information");
196  }
197  meshNames.push_back(std::string(mesh_name));
198  if (verb > 0) {
199  MOFEM_LOG_C("WORLD", Sev::inform, "Check mesh %s nsteps %d", mesh_name,
200  n_step);
201  }
202  }
203 
204  std::map<int, Range> family_elem_map;
205  std::map<string, Range> group_elem_map;
206 
207  for (unsigned int ii = 0; ii != meshNames.size(); ii++) {
208  CHKERR readMesh(file, ii, family_elem_map, verb);
209  CHKERR readFamily(file, ii, family_elem_map, group_elem_map, verb);
210  CHKERR makeBlockSets(group_elem_map, verb);
211  }
212 
213  if (MEDfileClose(fid) < 0) {
214  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
215  "Unable to close file '%s'", file.c_str());
216  }
217 
219 }
220 
221 static med_geometrie_element moab2med_element_type(const EntityType type) {
222  switch (type) {
223  case MBEDGE:
224  return MED_SEG2;
225  case MBTRI:
226  return MED_TRIA3;
227  case MBQUAD:
228  return MED_QUAD4;
229  case MBTET:
230  return MED_TETRA4;
231  case MBHEX:
232  return MED_HEXA8;
233  case MBPRISM:
234  return MED_PENTA6;
235  case MBPYRAMID:
236  return MED_PYRA5;
237  case MBVERTEX:
238  return MED_POINT1;
239  default:
240  return MED_NONE;
241  }
242 }
243 
244 MoFEMErrorCode MedInterface::readMesh(const string &file, const int index,
245  std::map<int, Range> &family_elem_map,
246  int verb) {
247 
248  Interface &m_field = cOre;
250 
251  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
252  if (fid < 0) {
253  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
254  "Unable to open file '%s'", file.c_str());
255  }
256  med_int v[3], vf[3];
257  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
258  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
259  if (verb > 1)
260  MOFEM_LOG_C("WORLD", Sev::noisy,
261  "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
262  vf[1], vf[2], v[0], v[1], v[2]);
263 
264  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
265  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
266  "Cannot read MED file older than V2.2");
267  }
268 
269  char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
270  bzero(mesh_name, MED_NAME_SIZE + 1);
271  bzero(mesh_desc, MED_COMMENT_SIZE + 1);
272  med_int space_dim;
273  med_mesh_type mesh_type;
274  med_int mesh_dim, n_step;
275  char dt_unit[MED_SNAME_SIZE + 1];
276  char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
277  med_sorting_type sorting_type;
278  med_axis_type axis_type;
279  if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
280  mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
281  axis_name, axis_unit) < 0) {
282  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
283  "Unable to read mesh information");
284  }
285  if (verb > 0)
286  MOFEM_LOG_C("WORLD", Sev::inform, "Reading mesh %s nsteps %d", mesh_name,
287  n_step);
288 
289  switch (axis_type) {
290  case MED_CARTESIAN:
291  break;
292  case MED_SPHERICAL:
293  case MED_CYLINDRICAL:
294  default:
295  SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
296  "Curvilinear coordinate system implemented");
297  }
298  if (space_dim < 2) {
299  SETERRQ1(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED,
300  "Not implemented for space dim %d", space_dim);
301  }
302 
303  EntityHandle mesh_meshset;
304  {
305  MeshsetsManager *meshsets_manager_ptr;
306  CHKERR m_field.getInterface(meshsets_manager_ptr);
307  int max_id = 0;
309  max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
310  }
311  max_id++;
312  CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id,
313  std::string(mesh_name));
314  CubitMeshSet_multiIndex::index<
315  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator cit;
316  cit =
317  meshsets_manager_ptr->getMeshsetsMultindex()
319  .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
320  max_id++;
321  mesh_meshset = cit->getMeshset();
322  meshMeshsets.push_back(mesh_meshset);
323  }
324 
325  med_bool change_of_coord, geo_transform;
326  med_int num_nodes = MEDmeshnEntity(
327  fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
328  MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
329  if (num_nodes < 0) {
330  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
331  "Could not read number of MED nodes");
332  }
333  if (num_nodes == 0) {
334  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
335  "No nodes in MED mesh");
336  }
337  if (verb > 0)
338  MOFEM_LOG_C("WORLD", Sev::inform, "Read number of nodes %d", num_nodes);
339 
340  std::vector<med_float> coord_med(space_dim * num_nodes);
341  if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
342  MED_NO_INTERLACE, &coord_med[0]) < 0) {
343  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
344  "Could not read MED node coordinates");
345  }
346 
347  // Add vertices to moab
348  ReadUtilIface *iface;
349  vector<double *> arrays_coord;
350  EntityHandle startv;
351  CHKERR m_field.get_moab().query_interface(iface);
352  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
353  Range verts(startv, startv + num_nodes - 1);
354  std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
355  arrays_coord[0]);
356  std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
357  arrays_coord[1]);
358  if (space_dim == 2) {
359  std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
360  } else {
361  std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
362  arrays_coord[2]);
363  }
364  CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
365  family_elem_map.clear();
366 
367  // get family for vertices
368  {
369  std::vector<med_int> nodes_tags(num_nodes, 0);
370  if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
371  MED_NODE, MED_NO_GEOTYPE,
372  &nodes_tags[0]) < 0) {
373  nodes_tags.clear();
374  // SETERRQ(
375  // m_field.get_comm(),
376  // MOFEM_OPERATION_UNSUCCESSFUL,
377  // "No family number for elements: using 0 as default family number"
378  // );
379  }
380  for (int i = 0; i < num_nodes; i++) {
381  // cerr << verts[i] << " " /*<< ele_tags[j] << " "*/ << nodes_tags[i] <<
382  // endl;
383  family_elem_map[nodes_tags.empty() ? i : nodes_tags[i]].insert(verts[i]);
384  }
385  }
386 
387  // read elements (loop over all possible MSH element types)
388  for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
389 
390  med_geometrie_element type = moab2med_element_type(ent_type);
391  if (type == MED_NONE)
392  continue;
393 
394  // get number of cells
395  med_bool change_of_coord;
396  med_bool geo_transform;
397  med_int num_ele = MEDmeshnEntity(
398  fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL, type, MED_CONNECTIVITY,
399  MED_NODAL, &change_of_coord, &geo_transform);
400 
401  // get connectivity
402  if (num_ele <= 0)
403  continue;
404  int num_nod_per_ele = type % 100;
405  std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
406  if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
407  MED_CELL, type, MED_NODAL,
408  MED_FULL_INTERLACE, &conn_med[0]) < 0) {
409  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
410  "Could not read MED elements");
411  }
412 
413  // cerr << "type " << ent_type << " ";
414  // cerr << "num_ele " << num_ele << " " << num_nod_per_ele << endl;;
415 
416  Range ents;
417 
418  if (ent_type != MBVERTEX) {
419  EntityHandle *conn_moab;
420  EntityHandle starte;
421  CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
422  starte, conn_moab);
423  switch (ent_type) {
424  // FIXME: Some connectivity could not work, need to run and test
425  case MBTET: {
426  int ii = 0;
427  for (int ee = 0; ee != num_ele; ee++) {
428  EntityHandle n[4];
429  for (int nn = 0; nn != num_nod_per_ele; nn++) {
430  // conn_moab[ii] = verts[conn_med[ii]-1];
431  n[nn] = verts[conn_med[ii + nn] - 1];
432  }
433  EntityHandle n0 = n[0];
434  n[0] = n[1];
435  n[1] = n0;
436  for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
437  conn_moab[ii] = n[nn];
438  }
439  }
440  } break;
441  default: {
442  int ii = 0;
443  for (int ee = 0; ee != num_ele; ee++) {
444  for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
445  // cerr << conn_med[ii] << " ";
446  conn_moab[ii] = verts[conn_med[ii] - 1];
447  }
448  // cerr << endl;
449  }
450  }
451  }
452  CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
453  conn_moab);
454  ents = Range(starte, starte + num_ele - 1);
455  } else {
456  // This is special case when in med vertices are defined as elements
457  int ii = 0;
458  std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
459  for (int ee = 0; ee != num_ele; ++ee)
460  for (int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
461  conn_moab[ii] = verts[conn_med[ii] - 1];
462  ents.insert_list(conn_moab.begin(), conn_moab.end());
463  }
464 
465  // Add elements to family meshset
466  CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
467 
468  // get family for cells
469  {
470  std::vector<med_int> fam(num_ele, 0);
471  if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
472  MED_CELL, type, &fam[0]) < 0) {
473  SETERRQ(
475  "No family number for elements: using 0 as default family number");
476  }
477  // std::vector<med_int> ele_tags(num_ele);
478  // if(MEDmeshEntityNumberRd(
479  // fid, mesh_name,MED_NO_DT,MED_NO_IT,MED_CELL,type,&ele_tags[0]) < 0
480  // ) {
481  // ele_tags.clear();
482  // }
483  for (int j = 0; j < num_ele; j++) {
484  // cerr << ents[j] << " " /*<< ele_tags[j] << " "*/ << fam[j] << endl;
485  family_elem_map[fam[j]].insert(ents[j]);
486  }
487  }
488  }
489 
490  if (MEDfileClose(fid) < 0) {
491  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
492  "Unable to close file '%s'", file.c_str());
493  }
494 
496 }
497 
499 MedInterface::readFamily(const string &file, const int index,
500  const std::map<int, Range> &family_elem_map,
501  std::map<string, Range> &group_elem_map, int verb) {
502  //
503  Interface &m_field = cOre;
505 
506  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
507  if (fid < 0) {
508  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
509  "Unable to open file '%s'", file.c_str());
510  }
511  med_int v[3], vf[3];
512  MEDlibraryNumVersion(&v[0], &v[1], &v[2]);
513  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
514 
515  if(verb>1) {
516  MOFEM_LOG_C("WORLD", Sev::noisy,
517  "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
518  vf[1], vf[2], v[0], v[1], v[2]);
519  }
520  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
521  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
522  "Cannot read MED file older than V2.2");
523  }
524 
525  // clear groups
526  group_elem_map.clear();
527 
528  med_int num_families = MEDnFamily(fid, meshNames[index].c_str());
529  if (num_families < 0) {
530  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
531  "Could not read MED families");
532  }
533  for (int i = 0; i < num_families; i++) {
534  med_int num_attrib =
535  (vf[0] == 2)
536  ? MEDnFamily23Attribute(fid, meshNames[index].c_str(), i + 1)
537  : 0;
538  med_int num_groups = MEDnFamilyGroup(fid, meshNames[index].c_str(), i + 1);
539  if (num_attrib < 0 || num_groups < 0) {
540  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
541  "Could not read MED groups or attributes");
542  }
543 
544  std::vector<med_int> attribId(num_attrib + 1);
545  std::vector<med_int> attrib_val(num_attrib + 1);
546  std::vector<char> attrib_des(MED_COMMENT_SIZE * num_attrib + 1);
547  std::vector<char> group_names(MED_LNAME_SIZE * num_groups + 1);
548  char family_name[MED_NAME_SIZE + 1];
549  med_int family_num;
550 
551  if (vf[0] == 2) { // MED2 file
552  if (MEDfamily23Info(fid, meshNames[index].c_str(), i + 1, family_name,
553  &attribId[0], &attrib_val[0], &attrib_des[0],
554  &family_num, &group_names[0]) < 0) {
555  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
556  "Could not read info for MED2 family %d", i + 1);
557  }
558  } else {
559  if (MEDfamilyInfo(fid, meshNames[index].c_str(), i + 1, family_name,
560  &family_num, &group_names[0]) < 0) {
561  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
562  "Could not read info for MED3 family %d", i + 1);
563  }
564  }
565 
566  // cerr << family_name << " " << family_num << " " << num_groups << endl;
567  for (int g = 0; g != num_groups; g++) {
568  std::string name =
569  std::string(&group_names[MED_LNAME_SIZE * g], MED_LNAME_SIZE - 1);
570  name.resize(NAME_TAG_SIZE - 1);
571  if (family_elem_map.find(family_num) == family_elem_map.end()) {
572  MOFEM_LOG_C(
573  "WORLD", Sev::warning,
574  "Family %d not read, likely type of element is not added "
575  "to moab database. Currently only triangle, quad, tetrahedral and "
576  "hexahedral elements are read to moab database",
577  family_num);
578  } else {
579  group_elem_map[name].merge(family_elem_map.at(family_num));
580  // cerr << string(&group_names[MED_LNAME_SIZE*g]) << endl;
581  }
582  }
583  }
584 
585  if (MEDfileClose(fid) < 0) {
586  SETERRQ1(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
587  "Unable to close file '%s'", file.c_str());
588  }
589 
591 }
592 
594 MedInterface::makeBlockSets(const std::map<string, Range> &group_elem_map,
595  int verb) {
596 
597  Interface &m_field = cOre;
599  MeshsetsManager *meshsets_manager_ptr;
600  CHKERR m_field.getInterface(meshsets_manager_ptr);
601 
602  int max_id = 0;
604  max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
605  }
606  max_id++;
607 
608  for (std::map<string, Range>::const_iterator git = group_elem_map.begin();
609  git != group_elem_map.end(); git++) {
610  CHKERR meshsets_manager_ptr->addMeshset(BLOCKSET, max_id, git->first);
611  CubitMeshSet_multiIndex::index<
612  Composite_Cubit_msId_And_MeshSetType_mi_tag>::type::iterator cit;
613  cit =
614  meshsets_manager_ptr->getMeshsetsMultindex()
616  .find(boost::make_tuple(max_id, CubitBCType(BLOCKSET).to_ulong()));
617  EntityHandle meshsets = cit->getMeshset();
618  if (!git->second.empty()) {
619  CHKERR m_field.get_moab().add_entities(meshsets, git->second);
620  }
621  max_id++;
622  }
623 
625  MOFEM_LOG("WORLD", Sev::verbose) << *cit;
626 
628 }
629 
632  if (medFileName.empty()) {
634  }
635  CHKERR readMed(medFileName, verb);
637 }
638 
639 MoFEMErrorCode MedInterface::writeMed(const string &file, int verb) {
640  Interface &m_field = cOre;
642  SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
644 }
645 
646 MoFEMErrorCode MedInterface::readFields(const std::string &file_name,
647  const std::string &field_name,
648  const bool load_series,
649  const int only_step, int verb) {
650 
651  Interface &m_field = cOre;
653  med_idt fid = MEDfileOpen((char *)file_name.c_str(), MED_LECTURE);
654  if (fid < 0) {
655  SETERRQ1(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
656  "Unable to open file '%s'", file_name.c_str());
657  }
658 
659  med_int num_comp = MEDfieldnComponentByName(fid, field_name.c_str());
660  if (num_comp <= 0) {
661  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
662  "Could not get number of components for MED field");
663  }
664 
665  char meshName[MED_NAME_SIZE + 1];
666  char dtUnit[MED_SNAME_SIZE + 1];
667  std::vector<char> compName(num_comp * MED_SNAME_SIZE + 1);
668  std::vector<char> compUnit(num_comp * MED_SNAME_SIZE + 1);
669  med_int numSteps = 0;
670  med_type_champ type;
671  med_bool localMesh;
672  if (MEDfieldInfoByName(fid, field_name.c_str(), meshName, &localMesh, &type,
673  &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
674  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
675  "Could not get MED field info");
676  }
677 
678  // Get meshset
679  MeshsetsManager *meshsets_manager_ptr;
680  CHKERR m_field.getInterface(meshsets_manager_ptr);
681  const CubitMeshSets *cubit_meshset_ptr;
682  CHKERR meshsets_manager_ptr->getCubitMeshsetPtr(meshName, &cubit_meshset_ptr);
683  EntityHandle meshset = cubit_meshset_ptr->getMeshset();
684 
685  int num_comp_msh = (num_comp <= 1)
686  ? 1
687  : (num_comp <= 3) ? 3 : (num_comp <= 9) ? 9 : num_comp;
688 
689  // Create tag to store nodal or cell values read form med file. Note that tag
690  // has prefix MED to avoid collison with other tags.
691  Tag th;
692  std::string tag_name = "MED_" + field_name;
693  {
694  std::vector<double> def_val(num_comp_msh, 0);
695  CHKERR m_field.get_moab().tag_get_handle(
696  tag_name.c_str(), num_comp_msh, MB_TYPE_DOUBLE, th,
697  MB_TAG_CREAT | MB_TAG_SPARSE, &def_val[0]);
698  }
699 
700  // Warning! The ordering of the elements in the last two lists is
701  // important: it should match the ordering of the MSH element types
702  // (when elements are saved without tags, this governs the order
703  // with which we implicitly index them in GModel::readMED)
704  const med_entity_type entType[] = {MED_NODE, MED_CELL, MED_NODE_ELEMENT};
705  const med_geometrie_element eleType[] = {
706  MED_NONE, MED_SEG2, MED_TRIA3, MED_QUAD4, MED_TETRA4, MED_HEXA8,
707  MED_PENTA6, MED_PYRA5, MED_SEG3, MED_TRIA6, MED_QUAD9, MED_TETRA10,
708  MED_HEXA27, MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
709  // const int nodesPerEle[] = {0, 2, 3, 4, 4, 8, 6, 5, 3, 6, 9, 10, 27, 1, 8,
710  // 20, 15, 13};
711 
712  std::vector<std::pair<int, int>> pairs;
713  for (unsigned int i = 0; i < sizeof(entType) / sizeof(entType[0]); i++) {
714  for (unsigned int j = 0; j < sizeof(eleType) / sizeof(eleType[0]); j++) {
715  if ((!i && !j) || j) {
716  med_int n = numSteps;
717  if (n > 0) {
718  pairs.push_back(std::pair<int, int>(i, j));
719  numSteps = std::max(numSteps, n);
720  }
721  if (!i && !j)
722  break; // MED_NOEUD does not care about eleType
723  }
724  }
725  }
726 
727  if (numSteps < 1 || pairs.empty()) {
728  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
729  "Nothing to import from MED file");
730  }
731 
732  if (load_series) {
733  SETERRQ(m_field.get_comm(), MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
734  }
735 
736  for (int step = (only_step == -1) ? 0 : only_step; step < numSteps; step++) {
737 
738  if (!load_series && only_step != step)
739  continue;
740 
741  // cerr << only_step << " " << step << endl;
742 
743  // FIXME: in MED3 we might want to loop over all profiles instead
744  // of relying of the default one
745 
746  // FIXME: MED3 allows to store multi-step meshes; we should
747 
748  for (unsigned int pair = 0; pair < pairs.size(); pair++) {
749 
750  // get step info
751  med_entite_maillage ent = entType[pairs[pair].first];
752  med_geometrie_element ele = eleType[pairs[pair].second];
753  med_int numdt, numit, ngauss;
754  med_float dt;
755  if (MEDfieldComputingStepInfo(fid, field_name.c_str(), step + 1, &numdt,
756  &numit, &dt) < 0) {
757  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
758  "Could not read step info");
759  }
760 
761  char locName[MED_NAME_SIZE + 1], profileName[MED_NAME_SIZE + 1];
762 
763  // get number of values in the field (numVal takes the number of
764  // Gauss points or the number of nodes per element into account,
765  // but not the number of components)
766  med_int profileSize;
767  med_int numVal = MEDfieldnValueWithProfile(
768  fid, field_name.c_str(), numdt, numit, ent, ele, 1,
769  MED_COMPACT_STMODE, profileName, &profileSize, locName, &ngauss);
770  numVal *= ngauss;
771 
772  if (numVal <= 0)
773  continue;
774 
775  // int mult = 1;
776  // if(ent == MED_NODE_ELEMENT) {
777  // mult = nodesPerEle[pairs[pair].second];
778  //}
779  // else if(ngauss != 1){
780  // mult = ngauss;
781  //}
782 
783  // read field data
784  std::vector<double> val(numVal * num_comp);
785  if (MEDfieldValueWithProfileRd(fid, field_name.c_str(), numdt, numit, ent,
786  ele, MED_COMPACT_STMODE, profileName,
787  MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
788  (unsigned char *)&val[0]) < 0) {
789  SETERRQ(m_field.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
790  "Could not read field values");
791  }
792 
793  // if (verb > 2) {
794  // // FIXME: This not looks ok for me
795  // cerr << ent << " " << ele << endl;
796  // cerr << string(meshName) << " : " << string(profileName) << " : "
797  // << string(locName) << " : " << profileSize << " : " << ngauss
798  // << endl;
799  // }
800 
801  switch (ent) {
802  case MED_CELL: {
803  EntityType ent_type = MBMAXTYPE;
804  switch (ele) {
805  case MED_TETRA4:
806  case MED_TETRA10:
807  ent_type = MBTET;
808  break;
809  case MED_HEXA8:
810  ent_type = MBHEX;
811  break;
812  default:
813  MOFEM_LOG_C("WORLD", Sev::warning,
814  "Not yet implemented for this cell %d", ele);
815  }
816  if (ent_type != MBMAXTYPE) {
817  if (ngauss == 1) {
818  Range ents;
819  CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
820  ents, true);
821  double e_vals[num_comp_msh];
822  bzero(e_vals, sizeof(double) * num_comp_msh);
823  std::vector<double>::iterator vit = val.begin();
824  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
825  for (int ii = 0; ii != num_comp; ii++, vit++) {
826  e_vals[ii] = *vit;
827  }
828  CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
829  }
830  } else {
831  Range ents;
832  CHKERR m_field.get_moab().get_entities_by_type(meshset, ent_type,
833  ents, true);
834  if (ents.size() * ngauss * num_comp != val.size()) {
835  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
836  "data inconsistency");
837  }
838  // FIXME simply average gauss values, far from perfect need fix
839  double e_vals[num_comp_msh];
840  std::vector<double>::iterator vit = val.begin();
841  for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
842  bzero(e_vals, sizeof(double) * num_comp_msh);
843  for (int gg = 0; gg != ngauss; gg++) {
844  for (int ii = 0; ii != num_comp; ii++, vit++) {
845  e_vals[ii] += *vit / ngauss;
846  }
847  }
848  CHKERR m_field.get_moab().tag_set_data(th, &*eit, 1, e_vals);
849  }
850  }
851  }
852  // SETERRQ1(
853  // m_field.get_comm(),
854  // MOFEM_NOT_IMPLEMENTED,
855  // "Not implemented for more gauss pts ngauss = %d",
856  // ngauss
857  // );
858  } break;
859  case MED_NODE:
860  case MED_NODE_ELEMENT:
861  default:
862  MOFEM_LOG_C("WORLD", Sev::inform, "Entity type %d not implemented",
863  ent);
864  }
865  }
866  }
867 
869 }
870 
871 std::ostream &operator<<(std::ostream &os,
872  const MedInterface::FieldData &field_data) {
873  os << "field name: " << field_data.fieldName;
874  os << " mesh name: " << field_data.meshName;
875  os << " local mesh: " << ((field_data.localMesh) ? "true" : "false");
876  os << std::endl;
877  // os << " field type: ";
878  // switch (field_data.fieldType) {
879  // case MED_FLOAT64: os << "MED_FLOAT64"; break;
880  // case MED_INT32: os << "MED_INT32"; break;
881  // case MED_INT64: os << "MED_INT64"; break;
882  // case MED_INT: os << "MED_INT"; break;
883  // };
884  os << " componentNames:";
885  for (unsigned int ff = 0; ff != field_data.componentNames.size(); ff++) {
886  os << " " << field_data.componentNames[ff];
887  }
888  os << std::endl;
889  os << " componentUnits:";
890  for (unsigned int ff = 0; ff != field_data.componentUnits.size(); ff++) {
891  os << " " << field_data.componentUnits[ff];
892  }
893  os << std::endl;
894  os << " dtUnit: " << field_data.dtUnit << endl;
895  os << " number of steps: " << field_data.ncSteps;
896  return os;
897 }
898 
899 } // namespace MoFEM
900 
901 #endif // WITH_MED
Deprecated interface functions.
MoFEM interface unique ID.
MoFEMErrorCode makeBlockSets(const std::map< string, Range > &group_elem_map, int verb=1)
make from groups meshsets
MoFEMErrorCode readMesh(const string &file, const int index, std::map< int, Range > &family_elem_map, int verb=1)
read mesh from MED file
virtual moab::Interface & get_moab()=0
std::map< std::string, FieldData > fieldNames
MoFEM::Core & cOre
core database
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:303
base class for all interface classes
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:306
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
this struct keeps basic methods for moab meshset about material and boundary conditions
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
Interface for managing meshsets containing materials and boundary conditions.
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
CubitMeshSet_multiIndex & getMeshsetsMultindex()
std::vector< std::string > componentNames
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
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
std::ostream & operator<<(std::ostream &os, const DataForcesAndSourcesCore::EntData &e)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
PetscBool flgFile
true if file name given in command line
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
add cubit meshset
MoFEMErrorCode getFileNameFromCommandLine(int verb=1)
Get MED file name from command line.
static med_geometrie_element moab2med_element_type(const EntityType type)
MoFEMErrorCode readMed(const string &file, int verb=1)
read MED file
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MedInterface(const MoFEM::Core &core)
std::vector< EntityHandle > meshMeshsets
meshset for each mesh
EntityHandle getMeshset() const
get bc meshset
#define MedFunctionBegin
std::vector< std::string > meshNames
list of meshes in MED file
std::vector< std::string > componentUnits
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
double FieldData
Field data type.
Definition: Types.hpp:36
std::bitset< 32 > CubitBCType
Definition: Types.hpp:62
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
virtual MPI_Comm & get_comm() const =0
Med file interface interface.
MoFEMErrorCode writeMed(const string &file, int verb=1)
write MED file
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)
static const MOFEMuuid IDD_MOFEMMedInterface
MoFEMErrorCode medGetFieldNames(const string &file, int verb=1)
Get field names in MED file.
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
std::string medFileName
MED file name.