149  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
 
  152             "Unable to open file '%s'", file.c_str());
 
  156  MEDlibraryNumVersion(&
v[0], &
v[1], &
v[2]);
 
  157  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
 
  161                "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
 
  162                vf[1], vf[2], 
v[0], 
v[1], 
v[2]);
 
  164  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
 
  166            "Cannot read MED file older than V2.2");
 
  169  for (
int i = 0; 
i < MEDnMesh(fid); 
i++) {
 
  170    char mesh_name[MED_NAME_SIZE + 1], mesh_desc[MED_COMMENT_SIZE + 1];
 
  171    bzero(mesh_name, MED_NAME_SIZE);
 
  172    bzero(mesh_desc, MED_COMMENT_SIZE);
 
  174    med_mesh_type mesh_type;
 
  175    med_int mesh_dim, n_step;
 
  176    char dt_unit[MED_SNAME_SIZE + 1];
 
  177    char axis_name[3 * MED_SNAME_SIZE + 1], axis_unit[3 * MED_SNAME_SIZE + 1];
 
  178    med_sorting_type sorting_type;
 
  179    med_axis_type axis_type;
 
  180    if (MEDmeshInfo(fid, 
i + 1, mesh_name, &space_dim, &mesh_dim, &mesh_type,
 
  181                    mesh_desc, dt_unit, &sorting_type, &n_step, &axis_type,
 
  182                    axis_name, axis_unit) < 0) {
 
  184              "Unable to read mesh information");
 
  186    meshNames.push_back(std::string(mesh_name));
 
  188      MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Check mesh %s nsteps %d", mesh_name,
 
  193  std::map<int, Range> family_elem_map;
 
  194  std::map<string, Range> group_elem_map;
 
  196  for (
unsigned int ii = 0; ii != 
meshNames.size(); ii++) {
 
  202  if (MEDfileClose(fid) < 0) {
 
  204             "Unable to close file '%s'", file.c_str());
 
 
  254                                      std::map<int, Range> &family_elem_map,
 
  260  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
 
  263             "Unable to open file '%s'", file.c_str());
 
  266  MEDlibraryNumVersion(&
v[0], &
v[1], &
v[2]);
 
  267  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
 
  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]);
 
  273  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
 
  275            "Cannot read MED file older than V2.2");
 
  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);
 
  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");
 
  295    MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Reading mesh %s nsteps %d", mesh_name,
 
  302  case MED_CYLINDRICAL:
 
  305            "Curvilinear coordinate system implemented");
 
  309             "Not implemented for space dim %d", space_dim);
 
  318      max_id = (max_id < cit->getMeshsetId()) ? cit->getMeshsetId() : max_id;
 
  322                                            std::string(mesh_name));
 
  323    CubitMeshSet_multiIndex::index<
 
  330    mesh_meshset = cit->getMeshset();
 
  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);
 
  340            "Could not read number of MED nodes");
 
  342  if (num_nodes == 0) {
 
  344            "No nodes in MED mesh");
 
  347    MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Read number of nodes %d", num_nodes);
 
  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");
 
  357  ReadUtilIface *iface;
 
  358  vector<double *> arrays_coord;
 
  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],
 
  365  std::copy(&coord_med[1 * num_nodes], &coord_med[2 * num_nodes],
 
  367  if (space_dim == 2) {
 
  368    std::fill(arrays_coord[2], &arrays_coord[2][num_nodes], 0.);
 
  370    std::copy(&coord_med[2 * num_nodes], &coord_med[3 * num_nodes],
 
  374  family_elem_map.clear();
 
  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) {
 
  384    for (
int i = 0; 
i < num_nodes; 
i++) {
 
  385      family_elem_map[nodes_tags.empty() ? 
i : nodes_tags[
i]].insert(verts[
i]);
 
  390  for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
 
  394    for (
auto type : types) {
 
  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);
 
  407      int num_nod_per_ele = type % 100;
 
  410            << 
"Reading elements " << num_ele << 
" of type " 
  411            << moab::CN::EntityTypeName(ent_type) << 
" number of nodes " 
  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");
 
  424      if (ent_type != MBVERTEX) {
 
  427        CHKERR iface->get_element_connect(num_ele, num_nod_per_ele, ent_type, 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];
 
  439            std::array<int, 10> nodes_map{
 
  448            for (
int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
 
  449              conn_moab[ii] = 
n[nodes_map[nn]];
 
  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];
 
  461            std::array<int, 20> nodes_map{
 
  462                0,  3,  2,  1,  4, 7, 6, 5, 
 
  469            for (
int nn = 0; nn != num_nod_per_ele; nn++, ii++) {
 
  470              conn_moab[ii] = 
n[nodes_map[nn]];
 
  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];
 
  483        CHKERR iface->update_adjacencies(starte, num_ele, num_nod_per_ele,
 
  485        ents = 
Range(starte, starte + num_ele - 1);
 
  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());
 
  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";
 
  507        for (
int j = 0; 
j < num_ele; 
j++) {
 
  508          family_elem_map[fam[
j]].insert(ents[
j]);
 
  514  if (MEDfileClose(fid) < 0) {
 
  516             "Unable to close file '%s'", file.c_str());
 
 
  524                         const std::map<int, Range> &family_elem_map,
 
  525                         std::map<string, Range> &group_elem_map, 
int verb) {
 
  529  med_idt fid = MEDfileOpen(file.c_str(), MED_ACC_RDONLY);
 
  532             "Unable to open file '%s'", file.c_str());
 
  535  MEDlibraryNumVersion(&
v[0], &
v[1], &
v[2]);
 
  536  MEDfileNumVersionRd(fid, &vf[0], &vf[1], &vf[2]);
 
  540                "Reading MED file V%d.%d.%d using MED library V%d.%d.%d", vf[0],
 
  541                vf[1], vf[2], 
v[0], 
v[1], 
v[2]);
 
  543  if (vf[0] < 2 || (vf[0] == 2 && vf[1] < 2)) {
 
  545            "Cannot read MED file older than V2.2");
 
  549  group_elem_map.clear();
 
  551  med_int num_families = MEDnFamily(fid, 
meshNames[index].c_str());
 
  552  if (num_families < 0) {
 
  554            "Could not read MED families");
 
  556  for (
int i = 0; 
i < num_families; 
i++) {
 
  559            ? MEDnFamily23Attribute(fid, 
meshNames[index].c_str(), 
i + 1)
 
  561    med_int num_groups = MEDnFamilyGroup(fid, 
meshNames[index].c_str(), 
i + 1);
 
  562    if (num_attrib < 0 || num_groups < 0) {
 
  564              "Could not read MED groups or attributes");
 
  567    std::vector<med_int> attribId(num_attrib + 1);
 
  568    std::vector<med_int> attrib_val(num_attrib + 1);
 
  569    std::vector<char> attrib_des(MED_COMMENT_SIZE * num_attrib + 1);
 
  570    std::vector<char> group_names(MED_LNAME_SIZE * num_groups + 1);
 
  571    char family_name[MED_NAME_SIZE + 1];
 
  575      if (MEDfamily23Info(fid, 
meshNames[index].c_str(), 
i + 1, family_name,
 
  576                          &attribId[0], &attrib_val[0], &attrib_des[0],
 
  577                          &family_num, &group_names[0]) < 0) {
 
  579                 "Could not read info for MED2 family %d", 
i + 1);
 
  582      if (MEDfamilyInfo(fid, 
meshNames[index].c_str(), 
i + 1, family_name,
 
  583                        &family_num, &group_names[0]) < 0) {
 
  585                 "Could not read info for MED3 family %d", 
i + 1);
 
  589    for (
int g = 0; 
g != num_groups; 
g++) {
 
  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()) {
 
  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",
 
  601        group_elem_map[name].merge(family_elem_map.at(family_num));
 
  606  if (MEDfileClose(fid) < 0) {
 
  608             "Unable to close file '%s'", file.c_str());
 
 
  700    boost::shared_ptr<std::vector<const CubitMeshSets *>> meshsets_ptr,
 
  701    boost::shared_ptr<Range> write_range_ptr, 
int verb) {
 
  704  moab::Interface &moab = m_field.
get_moab();
 
  708      << 
"WRITE_MED IS EXPERIMENTAL, MAY CONTAIN BUGS, ALWAYS CHECK THE OUTPUT " 
  712      MEDfileVersionOpen((
char *)file.c_str(), MED_ACC_CREAT, MED_MAJOR_NUM,
 
  713                         MED_MINOR_NUM, MED_RELEASE_NUM);
 
  717            "Cannot create MED file");
 
  722  CHKERR MEDfichDesEcr(fid, (
char *)
"MED file generated by MoFEM");
 
  725  char dtUnit[MED_SNAME_SIZE + 1] = 
"";       
 
  726  char axisName[3 * MED_SNAME_SIZE + 1] = 
""; 
 
  727  char axisUnit[3 * MED_SNAME_SIZE + 1] = 
""; 
 
  729  PetscBool is_cubit_meshset = PETSC_TRUE;
 
  730  int med_mesh_name_id = 0;
 
  733  char mesh_name_char[255];
 
  734  std::string mesh_name = 
"Mesh";
 
  737  PetscOptionsBegin(PETSC_COMM_WORLD, 
"", 
"MED mesh options", 
"");
 
  738  CHKERR PetscOptionsString(
"-med_mesh_name", 
"get med mesh name", 
"",
 
  739                            mesh_name.c_str(), mesh_name_char, 255, PETSC_NULLPTR);
 
  742  mesh_name = mesh_name_char;
 
  745  for (
auto &
m : *meshsets_ptr) {
 
  746    if (
m->getName() == mesh_name) {
 
  747      med_mesh_name_id = 
m->getMeshsetId();
 
  748      is_cubit_meshset = PETSC_FALSE;
 
  755  for (
auto &
m : *meshsets_ptr) {
 
  756    if (
m->getMeshsetId() == med_mesh_name_id)
 
  757      mesh_name = 
m->getName();
 
  758    max_id = (max_id < 
m->getMeshsetId()) ? 
m->getMeshsetId() : max_id;
 
  762  CHKERR MEDmeshCr(fid, mesh_name.c_str(), 3, 3, MED_UNSTRUCTURED_MESH,
 
  763                   "Mesh created", dtUnit, MED_SORT_DTIT, MED_CARTESIAN,
 
  767  med_int family_id = 0;
 
  768  std::map<std::vector<std::string>, std::tuple<med_int, std::vector<int>>>
 
  770  std::map<EntityHandle, med_int> entityHandle_family_map;
 
  773  shared_meshsets_map[std::vector<string>()] =
 
  774      std::make_tuple(family_id, std::vector<int>());
 
  778    std::string bc_type_name;
 
  779    if (iit->getBcTypeULong() & 
BLOCKSET) {
 
  781      bc_type_name = iit->getName();
 
  782      if (bc_type_name == 
"NoNameSet") {
 
  783        bc_type_name = 
"BLOCKSET_";
 
  784        bc_type_name += std::to_string(iit->getMeshsetId());
 
  787    } 
else if (iit->getBcTypeULong() & 
SIDESET ||
 
  788               iit->getBcTypeULong() & 
NODESET) {
 
  795        if ((iit->getBcType() & jj_bc_type).any()) {
 
  797          std::string temp_bc_type_name = string(
CubitBCNames[jj + 1]);
 
  798          if (temp_bc_type_name.length() > 10) {
 
  799            temp_bc_type_name = temp_bc_type_name.substr(0, 4);
 
  801          bc_type_name += temp_bc_type_name;
 
  806      bc_type_name += std::to_string(iit->getMeshsetId());
 
  808      bc_type_name = 
"UnknownSet";
 
  818  for (
auto &entity : *write_range_ptr) {
 
  820    std::vector<int> shared_meshsets;
 
  821    std::vector<string> shared_names;
 
  824    auto add_shared_meshset = [&](
auto &other_meshset) {
 
  825      std::string other_name = get_set_name(other_meshset);
 
  826      if (std::find(shared_names.begin(), shared_names.end(), other_name) ==
 
  827          shared_names.end()) {
 
  828        shared_meshsets.push_back(other_meshset->getMeshsetId());
 
  829        shared_names.push_back(other_name);
 
  834    for (
auto &other_meshset : *meshsets_ptr) {
 
  836      if (med_mesh_name_id == other_meshset->getMeshsetId())
 
  839      Range other_entities;
 
  841      moab.get_entities_by_handle(other_set, other_entities);
 
  842      if (other_entities.empty())
 
  845      bool is_in_meshset = moab.contains_entities(other_set, &entity, 1);
 
  848        add_shared_meshset(other_meshset);
 
  853    auto it = shared_meshsets_map.find(shared_names);
 
  854    if (it == shared_meshsets_map.end()) {
 
  858      it = shared_meshsets_map
 
  860                   {shared_names, std::make_tuple(family_id, shared_meshsets)})
 
  864    entityHandle_family_map[entity] = std::get<0>(it->second);
 
  868  for (
const auto &it : shared_meshsets_map) {
 
  870    std::string family_name = 
"F_";
 
  871    const auto &[family_id, shared_meshsets] = it.second;
 
  872    family_name += std::to_string(family_id);
 
  874    const auto &shared_meshset_names = it.first;
 
  875    std::string group_name;
 
  876    for (
const auto &name : shared_meshset_names) {
 
  878      std::string meshset_name = name;
 
  879      meshset_name.resize(MED_LNAME_SIZE, 
' ');
 
  880      group_name += meshset_name;
 
  884    CHKERR MEDfamilyCr(fid, mesh_name.c_str(), family_name.c_str(), family_id,
 
  885                       shared_meshset_names.size(), group_name.c_str());
 
  888        << 
"Creating family " << family_name << 
" with id " << family_id
 
  889        << 
" and " << shared_meshset_names.size() << 
" groups " << std::endl;
 
  894  verts = write_range_ptr->subset_by_type(MBVERTEX);
 
  897  std::vector<med_float> coord_med(3 * verts.size());
 
  898  std::vector<med_int> fam;
 
  899  std::vector<med_int> tags;
 
  900  std::map<EntityHandle, med_int> entityHandle_tag_map;
 
  901  med_int med_node_num = 1;
 
  904  for (Range::iterator it = verts.begin(); it != verts.end(); ++it) {
 
  906    moab.get_coords(&(*it), 1, coords);
 
  907    coord_med[3 * (it - verts.begin())] = coords[0];
 
  908    coord_med[3 * (it - verts.begin()) + 1] = coords[1];
 
  909    coord_med[3 * (it - verts.begin()) + 2] = coords[2];
 
  910    fam.push_back(entityHandle_family_map[*it]);
 
  913    entityHandle_tag_map[*it] = med_node_num;
 
  917  CHKERR MEDmeshNodeWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT,
 
  918                       MED_UNDEF_DT, MED_FULL_INTERLACE, verts.size(),
 
  919                       &coord_med[0], MED_FALSE, 
"", MED_TRUE, &tags[0],
 
  923  double last_tag = tags.back();
 
  925  for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
 
  928    moab.get_entities_by_type(0, ent_type, entities, 
true);
 
  930    ents_to_write = intersect(entities, *write_range_ptr);
 
  932    if (ents_to_write.empty())
 
  936    std::vector<med_int> tag_number;
 
  937    std::vector<med_int> family_number;
 
  938    std::vector<med_int> connectivity;
 
  941    for (
auto &entity : ents_to_write) {
 
  942      if (ent_type != MBVERTEX) {
 
  945        family_number.push_back(entityHandle_family_map[entity]);
 
  947        tag_number.push_back(last_tag);
 
  949        std::vector<EntityHandle> conn;
 
  950        moab.get_connectivity(&entity, 1, conn);
 
  952        for (
auto &
c : conn) {
 
  953          connectivity.push_back(entityHandle_tag_map[
c]);
 
  960    auto get_med_element_type = [](EntityType ent_type) {
 
  961      med_geometrie_element type;
 
  988    med_geometrie_element med_type = get_med_element_type(ent_type);
 
  989    if (ent_type == MBENTITYSET) {
 
  992    CHKERR MEDmeshElementWr(fid, mesh_name.c_str(), MED_NO_DT, MED_NO_IT, 0.,
 
  993                            MED_CELL, med_type, MED_NODAL, MED_FULL_INTERLACE,
 
  994                            family_number.size(), &connectivity[0], MED_FALSE,
 
  995                            nullptr, MED_TRUE, &tag_number[0], MED_TRUE,
 
  998    MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Writing %i elements of type %i (%s) ",
 
  999                family_number.size(), med_type,
 
 1000                moab::CN::EntityTypeName(ent_type));
 
 
 1010                                        const bool load_series,
 
 1011                                        const int only_step, 
int verb) {
 
 1015  med_idt fid = MEDfileOpen((
char *)file_name.c_str(), MED_LECTURE);
 
 1018             "Unable to open file '%s'", file_name.c_str());
 
 1021  med_int num_comp = MEDfieldnComponentByName(fid, 
field_name.c_str());
 
 1022  if (num_comp <= 0) {
 
 1024            "Could not get number of components for MED field");
 
 1027  char meshName[MED_NAME_SIZE + 1];
 
 1028  char dtUnit[MED_SNAME_SIZE + 1];
 
 1029  std::vector<char> compName(num_comp * MED_SNAME_SIZE + 1);
 
 1030  std::vector<char> compUnit(num_comp * MED_SNAME_SIZE + 1);
 
 1031  med_int numSteps = 0;
 
 1032  med_type_champ type;
 
 1034  if (MEDfieldInfoByName(fid, 
field_name.c_str(), meshName, &localMesh, &type,
 
 1035                         &compName[0], &compUnit[0], dtUnit, &numSteps) < 0) {
 
 1037            "Could not get MED field info");
 
 1050  int num_comp_msh = (num_comp <= 1)   ? 1
 
 1051                     : (num_comp <= 3) ? 3
 
 1052                     : (num_comp <= 9) ? 9
 
 1060    std::vector<double> def_val(num_comp_msh, 0);
 
 1062        tag_name.c_str(), num_comp_msh, MB_TYPE_DOUBLE, 
th,
 
 1063        MB_TAG_CREAT | MB_TAG_SPARSE, &def_val[0]);
 
 1070  const med_entity_type entType[] = {MED_NODE, MED_CELL, MED_NODE_ELEMENT};
 
 1071  const med_geometrie_element eleType[] = {
 
 1072      MED_NONE,   MED_SEG2,   MED_TRIA3, MED_QUAD4,  MED_TETRA4,  MED_HEXA8,
 
 1073      MED_PENTA6, MED_PYRA5,  MED_SEG3,  MED_TRIA6,  MED_QUAD9,   MED_TETRA10,
 
 1074      MED_HEXA27, MED_POINT1, MED_QUAD8, MED_HEXA20, MED_PENTA15, MED_PYRA13};
 
 1078  std::vector<std::pair<int, int>> pairs;
 
 1079  for (
unsigned int i = 0; 
i < 
sizeof(entType) / 
sizeof(entType[0]); 
i++) {
 
 1080    for (
unsigned int j = 0; 
j < 
sizeof(eleType) / 
sizeof(eleType[0]); 
j++) {
 
 1081      if ((!
i && !
j) || 
j) {
 
 1082        med_int 
n = numSteps;
 
 1084          pairs.push_back(std::pair<int, int>(
i, 
j));
 
 1085          numSteps = std::max(numSteps, 
n);
 
 1093  if (numSteps < 1 || pairs.empty()) {
 
 1095            "Nothing to import from MED file");
 
 1102  if (numSteps == 1 && only_step > 0) {
 
 1104            "Field has only one step. Cannot load step with -med_time_step, " 
 1105            "remove from command line.");
 
 1106  } 
else if (numSteps == only_step) {
 
 1108            "Cannot load final step %d. -med_time_step should be %d.",
 
 1109            only_step, only_step - 1);
 
 1112  for (
int step = (only_step == -1) ? 0 : only_step; step < numSteps; step++) {
 
 1113    MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Loading step %d of field %s", step + 1,
 
 1115    if (!load_series && only_step != step)
 
 1123    for (
unsigned int pair = 0; pair < pairs.size(); pair++) {
 
 1126      med_entite_maillage ent = entType[pairs[pair].first];
 
 1127      med_geometrie_element ele = eleType[pairs[pair].second];
 
 1128      med_int numdt, numit, ngauss;
 
 1130      if (MEDfieldComputingStepInfo(fid, 
field_name.c_str(), step + 1, &numdt,
 
 1133                "Could not read step info");
 
 1136      char locName[MED_NAME_SIZE + 1], profileName[MED_NAME_SIZE + 1];
 
 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);
 
 1151      std::vector<double> val(numVal * num_comp);
 
 1152      if (MEDfieldValueWithProfileRd(fid, 
field_name.c_str(), numdt, numit, ent,
 
 1153                                     ele, MED_COMPACT_STMODE, profileName,
 
 1154                                     MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
 
 1155                                     (
unsigned char *)&val[0]) < 0) {
 
 1157                "Could not read field values");
 
 1162        EntityType ent_type = MBMAXTYPE;
 
 1174                      "Not yet implemented for this cell %d", ele);
 
 1176        if (ent_type != MBMAXTYPE) {
 
 1179            CHKERR m_field.
get_moab().get_entities_by_type(meshset, ent_type,
 
 1181            double e_vals[num_comp_msh];
 
 1182            bzero(e_vals, 
sizeof(
double) * num_comp_msh);
 
 1183            std::vector<double>::iterator vit = val.begin();
 
 1184            for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
 
 1185              for (
int ii = 0; ii != num_comp; ii++, vit++) {
 
 1192            CHKERR m_field.
get_moab().get_entities_by_type(meshset, ent_type,
 
 1194            if (ents.size() * ngauss * num_comp != val.size()) {
 
 1196                      "data inconsistency");
 
 1199            double e_vals[num_comp_msh];
 
 1200            std::vector<double>::iterator vit = val.begin();
 
 1201            for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
 
 1202              bzero(e_vals, 
sizeof(
double) * num_comp_msh);
 
 1203              for (
int gg = 0; gg != ngauss; gg++) {
 
 1204                for (
int ii = 0; ii != num_comp; ii++, vit++) {
 
 1205                  e_vals[ii] += *vit / ngauss;
 
 1214      case MED_NODE_ELEMENT: {
 
 1215        EntityType ent_type = MBVERTEX;
 
 1217        CHKERR m_field.
get_moab().get_entities_by_type(meshset, ent_type, ents,
 
 1219        double e_vals[num_comp_msh];
 
 1220        bzero(e_vals, 
sizeof(
double) * num_comp_msh);
 
 1221        std::vector<double>::iterator vit = val.begin();
 
 1222        for (Range::iterator eit = ents.begin(); eit != ents.end(); eit++) {
 
 1223          for (
int ii = 0; ii != num_comp; ii++, vit++) {
 
 1230        MOFEM_LOG_C(
"MEDWORLD", Sev::inform, 
"Entity type %d not implemented",