254 med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
257 "Unable to open file '%s'", file.c_str());
260 MEDlibraryNumVersion(&
v[0], &
v[1], &
v[2]);
261 MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
264 "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
265 vf[1], vf[2],
v[0],
v[1],
v[2]);
267 if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
269 "Cannot read MED file older than V2.2");
272 char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
273 bzero(mesh_name, MED_NAME_SIZE + 1);
274 bzero(mesh_desc, MED_COMMENT_SIZE + 1);
276 med_mesh_type mesh_type;
277 med_int mesh_dim, n_step;
278 char dt_unit[MED_SNAME_SIZE + 1];
279 char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
280 med_sorting_type sorting_type;
281 med_axis_type axis_type;
282 if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
283 mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
284 axis_name, axis_unit) < 0) {
286 "Unable to read mesh information");
289 MOFEM_LOG_C(
"MEDWORLD", Sev::inform,
"Reading mesh %s nsteps %d", mesh_name,
296 case MED_CYLINDRICAL:
299 "Curvilinear coordinate system implemented");
303 "Not implemented for space dim %d", space_dim);
308 MeshsetsManager *meshsets_manager_ptr;
309 CHKERR m_field.getInterface(meshsets_manager_ptr);
312 max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
316 std::string(mesh_name));
317 CubitMeshSet_multiIndex::index<
318 Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
320 meshsets_manager_ptr->getMeshsetsMultindex()
321 .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
324 mesh_meshset = cit->getMeshset();
328 med_bool change_of_coord, geo_transform;
329 med_int num_nodes = MEDmeshnEntity(
330 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
331 MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
334 "Could not read number of MED nodes");
336 if (num_nodes == 0) {
338 "No nodes in MED mesh");
341 MOFEM_LOG_C(
"MEDWORLD", Sev::inform,
"Read number of nodes %d", num_nodes);
343 std::vector<med_float> coord_med(space_dim * num_nodes);
344 if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
345 MED_NO_INTERLACE, &coord_med[0]) < 0) {
347 "Could not read MED node coordinates");
351 ReadUtilIface *iface;
352 vector<double *> arrays_coord;
354 CHKERR m_field.get_moab().query_interface(iface);
355 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
356 Range verts(startv, startv + num_nodes - 1);
357 std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
359 std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
361 if (space_dim == 2) {
362 std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
364 std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
367 CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
368 family_elem_map.clear();
372 std::vector<med_int> nodes_tags(num_nodes, 0);
373 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
374 MED_NODE, MED_NO_GEOTYPE,
375 &nodes_tags[0]) < 0) {
383 for (
int i = 0;
i < num_nodes;
i++) {
386 family_elem_map[nodes_tags.empty() ?
i : nodes_tags[
i]].insert(verts[
i]);
391 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
395 for (
auto type : types) {
398 med_bool change_of_coord;
399 med_bool geo_transform;
400 med_int num_ele = MEDmeshnEntity(
401 fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL,
type,
402 MED_CONNECTIVITY, MED_NODAL, &change_of_coord, &geo_transform);
408 int num_nod_per_ele =
type % 100;
411 <<
"Reading elements " << num_ele <<
" of type "
412 << moab::CN::EntityTypeName(ent_type) <<
" number of nodes "
415 std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
416 if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
417 MED_CELL,
type, MED_NODAL,
418 MED_FULL_INTERLACE, &conn_med[0]) < 0) {
420 "Could not read MED elements");
428 if (ent_type != MBVERTEX) {
431 CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
437 for (
int ee = 0; ee != num_ele; ee++) {
439 for (
int nn = 0; nn != num_nod_per_ele; nn++) {
441 n[nn] = verts[conn_med[ii + nn] - 1];
444 std::array<int, 10> nodes_map{
453 for (
int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
454 conn_moab[ii] =
n[nodes_map[nn]];
460 for (
int ee = 0; ee != num_ele; ee++) {
461 for (
int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
463 conn_moab[ii] = verts[conn_med[ii] - 1];
469 CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
471 ents =
Range(starte, starte + num_ele - 1);
475 std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
476 for (
int ee = 0; ee != num_ele; ++ee)
477 for (
int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
478 conn_moab[ii] = verts[conn_med[ii] - 1];
479 ents.insert_list(conn_moab.begin(), conn_moab.end());
483 CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
487 std::vector<med_int> fam(num_ele, 0);
488 if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
489 MED_CELL,
type, &fam[0]) < 0) {
491 "No family number for elements: using 0 as default family "
500 for (
int j = 0;
j < num_ele;
j++) {
502 family_elem_map[fam[
j]].insert(ents[
j]);
508 if (MEDfileClose(fid) < 0) {
510 "Unable to close file '%s'", file.c_str());
Public Member Functions inherited from MoFEM::UnknownInterface