12                                 {
   14 
   15  try {
   16 
   17    moab::Core mb_instance;
   18    moab::Interface &moab = mb_instance;
   19    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, 
MYPCOMM_INDEX);
 
   20    if (pcomm == NULL)
   21      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
   22 
   23    char mesh_out_file[255] = "out.h5m";
   25    char meshset_to_write[255] = "";
   26    PetscBool meshset_to_write_set = PETSC_FALSE;
   27    int time_step = 0;
   28    
   31    PetscBool bit_ref_set = PETSC_FALSE;
   32    int bit_ref_level = 0;
   33 
   34    PetscOptionsBegin(m_field.
get_comm(), 
"", 
"Read MED tool", 
"none");
 
   36                                 255, PETSC_NULLPTR);
   37    CHKERR PetscOptionsInt(
"-bit_ref_level", 
"bit ref level", 
"", bit_ref_level,
 
   38                           &bit_ref_level, &bit_ref_set);
   40                                 meshset_to_write, 255, &meshset_to_write_set);
   41    CHKERR PetscOptionsInt(
"-med_time_step", 
"time step", 
"", time_step,
 
   42                           &time_step, PETSC_NULLPTR);
   43    CHKERR PetscOptionsString(
"-output_file", 
"output mesh file name", 
"",
 
   44                              "out.h5m", mesh_out_file, 255, PETSC_NULLPTR);
   45    PetscOptionsEnd();
   46 
   47    const char *option = "";
   48 
   49    
   51 
   52    
   54 
   55    
   56    Range entities_to_write;
 
   57    if (bit_ref_set) {
   58      
   60      
   61      bit_level.set(bit_ref_level);
   62 
   63      
   66 
   67      
   68      entities_to_write = subtract(entities_to_write,
   69                                   entities_to_write.subset_by_type(MBPRISM));
   70      
   71      entities_to_write =
   72          subtract(entities_to_write, entities_to_write.subset_by_type(MBQUAD));
   74          << "Total number of entities in bit ref level " << bit_ref_level
   75          << ": " << entities_to_write.size() << std::endl;
   76 
   78          << "Removing interior edges and face entities"<< std::endl;
   79      
   80      Range entities_3d = entities_to_write.subset_by_dimension(3);
 
   81      Range entities_2d = entities_to_write.subset_by_dimension(2);
 
   82      Range entities_1d = entities_to_write.subset_by_dimension(1);
 
   83 
   84      
   87      Skinner skin(&moab);
   88      CHKERR skin.find_skin(0, entities_3d, 
false, skin_2d);
 
   89      CHKERR moab.get_adjacencies(skin_2d, 1, 
false, skin_1d,
 
   90                                  moab::Interface::UNION);
   91 
   92      
   93      entities_2d = subtract(entities_2d, skin_2d);
   94      entities_to_write = subtract(entities_to_write, entities_2d);
   95 
   96      entities_1d = subtract(entities_1d, skin_1d);
   97      entities_to_write = subtract(entities_to_write, entities_1d);
   98 
  100          << "Removing nodes with no connectivity"<< std::endl;
  101      
  103      
  104      nodes = entities_to_write.subset_by_type(MBVERTEX);
  105      
  107      for (auto node : nodes) {
  110        CHKERR moab.get_adjacencies(&node_handle, 1, 1, 
false, adj_edges,
 
  111                                    moab::Interface::UNION);
  112        if (adj_edges.empty()) {
  113          free_nodes.insert(node);
  114        }
  115      }
  116      
  117      for (auto node : free_nodes) {
  118        double coords[3];
  119        CHKERR moab.get_coords(&node, 1, coords);
 
  120      }
  121 
  122      
  123      entities_to_write = subtract(entities_to_write, free_nodes);
  124 
  125      
  126      
  129 
  131                                                                 meshset_10000);
  133                                                                 meshset_10001);
  134 
  135      Range entities_10000;
 
  136      Range entities_10001;
 
  137      moab.get_entities_by_handle(meshset_10000, entities_10000);
  138      moab.get_entities_by_handle(meshset_10001, entities_10001);
  139      
  140      auto check_face_orientation = [&moab](
const Range &faces) {
 
  141        for (auto face : faces) {
  142          
  143          
  145          CHKERR(moab.get_adjacencies(&face, 1, 3, 
false, adj_tets,
 
  146                                      moab::Interface::UNION));
  147          
  149          actual_tets = adj_tets.subset_by_type(MBTET);
  150 
  151          int side_number;
  152          int sense;
  153          int offset;
  154          for (auto tet : actual_tets) {
  155            CHKERR(moab.side_number(tet, face, side_number, sense, offset));
 
  156            if (sense == -1) {
  157              
  158              std::vector<EntityHandle> nodes_face;
  159              CHKERR(moab.get_connectivity(&face, 1, nodes_face));
 
  160              std::vector<EntityHandle> face_nodes_new;
  161              
  162              face_nodes_new = {nodes_face[1], nodes_face[0], nodes_face[2]};
  163              CHKERR(moab.set_connectivity(face, face_nodes_new.data(), 3));
 
  164              
  165              CHKERR(moab.side_number(tet, face, side_number, sense, offset));
 
  166              if (sense == -1) {
  168                    << "Face: " << face << " has wrong orientation";
  169              }
  170            }
  171          }
  172        }
  173      };
  174 
  175      check_face_orientation(entities_10000);
  176      check_face_orientation(entities_10001);
  177 
  178    } else if (meshset_to_write_set) {
  179      
  180      auto &meshsets_idx =
  182 
  183      for (auto &meshset : meshsets_idx) {
  184        auto meshset_name = meshset.getName();
  185        auto trim_name = [&](std::string &name) {
  186          name.erase(std::remove(name.begin(), name.end(), ' '), name.end());
  187        };
  188        trim_name(meshset_name);
  189 
  190        if (meshset_to_write == meshset_name) {
  191          CHKERR moab.get_entities_by_handle(meshset.getMeshset(),
 
  192                                             entities_to_write, true);
  193          break;
  194        }
  195      }
  197          << "Wrtiting all entitiies from meshset: " << meshset_to_write
  198          << std::endl;
  199    } else {
  200      
  201      CHKERR moab.get_entities_by_handle(0, entities_to_write, 
true);
 
  203          << "Wrtiting all entitiies" << std::endl;
  204    }
  206        << "Number of entities to write: " << entities_to_write.size()
  207        << std::endl;
  208    
  209    for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
  211      moab.get_entities_by_type(0, ent_type, entities, true);
  213      ents_to_write = intersect(entities, entities_to_write);
  214 
  215      if (ents_to_write.empty())
  216        continue;
  217 
  219          << "Number of entities to write: " << ents_to_write.size()
  220          << " of type: " << moab::CN::EntityTypeName(ent_type)
  221          << " from total of: " << entities.size() << std::endl;
  222    }
  223 
  224    
  225    boost::shared_ptr<Range> entities_to_write_ptr =
  226        boost::make_shared<Range>(entities_to_write);
  227 
  231  }
  233 
  235 
  236  return 0;
  237}
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
virtual MPI_Comm & get_comm() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Interface for load MED files.
MoFEMErrorCode writeMed(boost::shared_ptr< Range > range_ptr=nullptr, int verb=1)
write MED file
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.