255                                                {
  256 
  259 
  260  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
  261  if (fid < 0) {
  263             "Unable to open file '%s'", file.c_str());
  264  }
  266  MEDlibraryNumVersion(&
v[0], &
v[1], &
v[2]);
 
  267  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
  268  if (verb > 1)
  270                "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
  271                vf[1], vf[2], 
v[0], 
v[1], 
v[2]);
 
  272 
  273  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
  275            "Cannot read MED file older than V2.2");
  276  }
  277 
  278  char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
  279  bzero(mesh_name, MED_NAME_SIZE + 1);
  280  bzero(mesh_desc, MED_COMMENT_SIZE + 1);
  281  med_int space_dim;
  282  med_mesh_type mesh_type;
  283  med_int mesh_dim, n_step;
  284  char dt_unit[MED_SNAME_SIZE + 1];
  285  char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
  286  med_sorting_type sorting_type;
  287  med_axis_type axis_type;
  288  if (MEDmeshInfo(fid, index + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
  289                  mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
  290                  axis_name, axis_unit) < 0) {
  292            "Unable to read mesh information");
  293  }
  294  if (verb > 0)
  295    MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Reading mesh %s nsteps %d", mesh_name,
 
  296                n_step);
  297 
  298  switch (axis_type) {
  299  case MED_CARTESIAN:
  300    break;
  301  case MED_SPHERICAL:
  302  case MED_CYLINDRICAL:
  303  default:
  305            "Curvilinear coordinate system implemented");
  306  }
  307  if (space_dim < 2) {
  309             "Not implemented for space dim %d", space_dim);
  310  }
  311 
  313  {
  314    MeshsetsManager *meshsets_manager_ptr;
  315    CHKERR m_field.getInterface(meshsets_manager_ptr);
 
  316    int max_id = 0;
  318      max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
  319    }
  320    max_id++;
  322                                            std::string(mesh_name));
  323    CubitMeshSet_multiIndex::index<
  324        Composite_Cubit_msId_And_MeshsetType_mi_tag>::type::iterator cit;
  325    cit =
  326        meshsets_manager_ptr->getMeshsetsMultindex()
  327            .get<Composite_Cubit_msId_And_MeshsetType_mi_tag>()
  329    max_id++;
  330    mesh_meshset = cit->getMeshset();
  332  }
  333 
  334  med_bool change_of_coord, geo_transform;
  335  med_int num_nodes = MEDmeshnEntity(
  336      fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
  337      MED_COORDINATE, MED_NO_CMODE, &change_of_coord, &geo_transform);
  338  if (num_nodes < 0) {
  340            "Could not read number of MED nodes");
  341  }
  342  if (num_nodes == 0) {
  344            "No nodes in MED mesh");
  345  }
  346  if (verb > 0)
  347    MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Read number of nodes %d", num_nodes);
 
  348 
  349  std::vector<med_float> coord_med(space_dim * num_nodes);
  350  if (MEDmeshNodeCoordinateRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
  351                              MED_NO_INTERLACE, &coord_med[0]) < 0) {
  353            "Could not read MED node coordinates");
  354  }
  355 
  356  
  357  ReadUtilIface *iface;
  358  vector<double *> arrays_coord;
  360  CHKERR m_field.get_moab().query_interface(iface);
 
  361  CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays_coord);
 
  362  Range verts(startv, startv + num_nodes - 1);
 
  363  std::copy(&coord_med[0 * num_nodes], &coord_med[1 * num_nodes],
  364            arrays_coord[0]);
  365  std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
  366            arrays_coord[1]);
  367  if (space_dim == 2) {
  368    std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
  369  } else {
  370    std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
  371              arrays_coord[2]);
  372  }
  373  CHKERR m_field.get_moab().add_entities(mesh_meshset, verts);
 
  374  family_elem_map.clear();
  375 
  376  
  377  {
  378    std::vector<med_int> nodes_tags(num_nodes, 0);
  379    if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
  380                                    MED_NODE, MED_NO_GEOTYPE,
  381                                    &nodes_tags[0]) < 0) {
  382      nodes_tags.clear();
  383    }
  384    for (
int i = 0; 
i < num_nodes; 
i++) {
 
  385      family_elem_map[nodes_tags.empty() ? 
i : nodes_tags[
i]].insert(verts[
i]);
 
  386    }
  387  }
  388 
  389  
  390  for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
  391 
  393 
  394    for (auto type : types) {
  395 
  396      
  397      med_bool change_of_coord;
  398      med_bool geo_transform;
  399      med_int num_ele = MEDmeshnEntity(
  400          fid, mesh_name, MED_NO_DT, MED_NO_IT, MED_CELL, type,
  401          MED_CONNECTIVITY, MED_NODAL, &change_of_coord, &geo_transform);
  402 
  403      
  404      if (num_ele <= 0)
  405        continue;
  406 
  407      int num_nod_per_ele = 
type % 100;
 
  408      if (verb > 0)
  410            << "Reading elements " << num_ele << " of type "
  411            << moab::CN::EntityTypeName(ent_type) << " number of nodes "
  412            << num_nod_per_ele;
  413 
  414      std::vector<med_int> conn_med(num_ele * num_nod_per_ele);
  415      if (MEDmeshElementConnectivityRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
  416                                       MED_CELL, type, MED_NODAL,
  417                                       MED_FULL_INTERLACE, &conn_med[0]) < 0) {
  419                "Could not read MED elements");
  420      }
  421 
  423 
  424      if (ent_type != MBVERTEX) {
  427        CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 0,
 
  428                                          starte, conn_moab);
  429        switch (ent_type) {
  430        
  431        case MBTET: {
  432          int ii = 0;
  433          for (int ee = 0; ee != num_ele; ee++) {
  434            std::vector<EntityHandle> 
n(num_nod_per_ele);
 
  435            for (int nn = 0; nn != num_nod_per_ele; nn++) {
  436              n[nn] = verts[conn_med[ii + nn] - 1];
 
  437            }
  438 
  439            std::array<int, 10> nodes_map{
  440 
  441                1, 0, 2, 3,
  442 
  443                4, 6, 5, 8,
  444                7, 9
  445 
  446            };
  447 
  448            for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
  449              conn_moab[ii] = 
n[nodes_map[nn]];
 
  450            }
  451          }
  452        }
  453        case MBHEX: {
  454          int ii = 0;
  455          for (int ee = 0; ee != num_ele; ee++) {
  456            std::vector<EntityHandle> 
n(num_nod_per_ele);
 
  457            for (int nn = 0; nn != num_nod_per_ele; nn++) {
  458              n[nn] = verts[conn_med[ii + nn] - 1];
 
  459            }
  460 
  461            std::array<int, 20> nodes_map{
  462                0,  3,  2,  1,  4, 7, 6, 5, 
  463                11, 10, 9,  8,              
  464                16, 19, 18, 17,             
  465                15, 14, 13, 12              
  466 
  467            };
  468 
  469            for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
  470              conn_moab[ii] = 
n[nodes_map[nn]];
 
  471            }
  472          }
  473        } break;
  474        default: {
  475          int ii = 0;
  476          for (int ee = 0; ee != num_ele; ee++) {
  477            for (int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
  478              conn_moab[ii] = verts[conn_med[ii] - 1];
  479            }
  480          }
  481        }
  482        }
  483        CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
 
  484                                         conn_moab);
  485        ents = 
Range(starte, starte + num_ele - 1);
 
  486      } else {
  487        
  488        int ii = 0;
  489        std::vector<EntityHandle> conn_moab(num_ele * num_nod_per_ele);
  490        for (int ee = 0; ee != num_ele; ++ee)
  491          for (int nn = 0; nn != num_nod_per_ele; ++nn, ++ii)
  492            conn_moab[ii] = verts[conn_med[ii] - 1];
  493        ents.insert_list(conn_moab.begin(), conn_moab.end());
  494      }
  495 
  496      
  497      CHKERR m_field.get_moab().add_entities(mesh_meshset, ents);
 
  498 
  499      
  500      {
  501        std::vector<med_int> fam(num_ele, 0);
  502        if (MEDmeshEntityFamilyNumberRd(fid, mesh_name, MED_NO_DT, MED_NO_IT,
  503                                        MED_CELL, type, &fam[0]) < 0) {
  505                  "No family number for elements: using 0 as default family number";
  506        }
  507        for (
int j = 0; 
j < num_ele; 
j++) {
 
  508          family_elem_map[fam[
j]].insert(ents[
j]);
 
  509        }
  510      }
  511    }
  512  }
  513 
  514  if (MEDfileClose(fid) < 0) {
  516             "Unable to close file '%s'", file.c_str());
  517  }
  518 
  520}
static std::vector< med_geometrie_element > moab2med_element_type(const EntityType type)
std::vector< EntityHandle > meshMeshsets
meshset for each mesh