19                                 {
   20 
   22 
   23  try {
   24 
   25    
   27    char mesh_out_file[255] = "out.h5m";
   28    int dim = 3;
   30    PetscBool shift = PETSC_TRUE;
   31    PetscBool flg_file = PETSC_FALSE;
   32    PetscBool 
debug = PETSC_FALSE;
 
   33 
   34    PetscOptionsBegin(PETSC_COMM_WORLD, "", "none", "none");
   35 
   36    CHKERR PetscOptionsString(
"-my_file", 
"mesh file name", 
"", 
"mesh.h5m",
 
   38    if (flg_file != PETSC_TRUE)
   39      CHKERR PetscOptionsString(
"-file_name", 
"mesh file name", 
"", 
"mesh.h5m",
 
   41    CHKERR PetscOptionsString(
"-output_file", 
"output mesh file name", 
"",
 
   42                              "mesh.h5m", mesh_out_file, 255, PETSC_NULLPTR);
   43    CHKERR PetscOptionsInt(
"-dim", 
"entities dim", 
"", dim, &dim, PETSC_NULLPTR);
 
   44    CHKERR PetscOptionsInt(
"-nb_levels", 
"number of refinement levels", 
"",
 
   46    CHKERR PetscOptionsBool(
"-shift",
 
   47                            "shift bits, squash entities of refined elements",
   48                            "", shift, &shift, PETSC_NULLPTR);
   49    CHKERR PetscOptionsBool(
"-debug",
 
   50                            "write additional files with bit ref levels", "",
   52 
   53    PetscOptionsEnd();
   54 
   55    moab::Core mb_instance;
   56    moab::Interface &moab = mb_instance;
   57    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, 
MYPCOMM_INDEX);
 
   58    if (pcomm == NULL)
   59      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
   60 
   61    const char *option;
   62    option = "";
   64 
   65    
   68 
   69    if (flg_file != PETSC_TRUE) {
   71              "*** ERROR -my_file (-file_name) (mesh file needed)");
   72    }
   73 
   76 
   78 
   80 
   85 
   87      CHKERR moab.get_adjacencies(ents, 1, 
true, edges, moab::Interface::UNION);
 
   88 
   90      CHKERR moab.create_meshset(MESHSET_SET, meshset_ref_edges);
 
   91      CHKERR moab.add_entities(meshset_ref_edges, edges);
 
   92 
   94 
   97      if (dim == 3) {
   99      } else if (dim == 2) {
  101      } else {
  103                "Refinement implemented only for three and two dimensions");
  104      }
  105 
  106      auto update_meshsets = [&]() {
  108        
  112              cubit_meshset, 
bit(
l + 1), cubit_meshset, MBMAXTYPE, 
true);
 
  113        }
  115      };
  116 
  117      auto update_partition_sets = [&]() {
  119 
  120        ParallelComm *pcomm = ParallelComm::get_pcomm(
  122        Tag part_tag = pcomm->part_tag();
 
  123 
  126            0, MBENTITYSET, &part_tag, NULL, 1, r_tagged_sets,
  127            moab::Interface::UNION);
  128 
  129        std::vector<EntityHandle> tagged_sets(r_tagged_sets.size());
  130        std::copy(r_tagged_sets.begin(), r_tagged_sets.end(),
  131                  tagged_sets.begin());
  132 
  133        auto order_tagged_sets = [&]() {
  135          std::vector<int> parts(tagged_sets.size());
  137              part_tag, &*tagged_sets.begin(), tagged_sets.size(),
  138              &*parts.begin());
  139          map<int, EntityHandle> m_tagged_sets;
  140          for (int p = 0; p != tagged_sets.size(); ++p) {
  141            m_tagged_sets[parts[p]] = tagged_sets[p];
  142          }
  143          for (int p = 0; p != tagged_sets.size(); ++p) {
  144            tagged_sets[p] = m_tagged_sets.at(p);
  145          }
  147        };
  148 
  149        auto add_children = [&]() {
  151          std::vector<Range> part_ents(tagged_sets.size());
  152 
  153          for (int p = 0; p != tagged_sets.size(); ++p) {
  155            CHKERR moab.get_entities_by_dimension(tagged_sets[p], dim, ents,
 
  156                                                  true);
  159 
  162            children = children.subset_by_dimension(dim);
  165 
  167            for (
auto d = 1; 
d != dim; ++
d) {
 
  168              CHKERR moab.get_adjacencies(children.subset_by_dimension(dim), d,
 
  169                                          false, adj, moab::Interface::UNION);
  170            }
  171 
  172            part_ents[p].merge(children);
  173            part_ents[p].merge(adj);
  174          }
  175 
  176          for (int p = 1; p != tagged_sets.size(); ++p) {
  177            for (int pp = 0; pp != p; pp++) {
  178              part_ents[p] = subtract(part_ents[p], part_ents[pp]);
  179            }
  180          }
  181 
  182          for (int p = 0; p != tagged_sets.size(); ++p) {
  183            CHKERR moab.add_entities(tagged_sets[p], part_ents[p]);
 
  184            CHKERR moab.tag_clear_data(part_tag, part_ents[p], &p);
 
  185          }
  186 
  188 
  194                                                   meshset_ptr->get_ptr(), 1);
  196            };
  197 
  198            for (int p = 0; p != tagged_sets.size(); ++p) {
  200                  << 
"Write part " << p << 
" level " << 
l;
 
  203                                                               ents, true);
  207                                    "_" + boost::lexical_cast<std::string>(
l) +
 
  208                                    ".vtk",
  209                                ents);
  210            }
  211          }
  212 
  214        };
  215 
  216        CHKERR order_tagged_sets();
 
  218 
  220      };
  221 
  222      CHKERR update_partition_sets();
 
  224 
  225      CHKERR moab.delete_entities(&meshset_ref_edges, 1);
 
  226    }
  227 
  230        MOFEM_LOG(
"WORLD", Sev::inform) << 
"Write level " << 
l;
 
  233            (
"level" + boost::lexical_cast<std::string>(
l) + 
".vtk").c_str(),
 
  234            "VTK", "");
  235      }
  236    }
  237 
  238    if (shift == PETSC_TRUE) {
  239      MOFEM_LOG(
"WORLD", Sev::inform) << 
"Shift bits";
 
  242    }
  243 
  244    CHKERR moab.write_file(mesh_out_file);
 
  245  }
  247 
  250 
  251  return 0;
  252}
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
MoFEMErrorCode updateMeshsetByEntitiesChildren(const EntityHandle parent, const BitRefLevel &parent_bit, const BitRefLevel &parent_mask, const BitRefLevel &child_bit, const BitRefLevel &child_mask, const EntityHandle child, EntityType child_type, const bool recursive=false, int verb=0)
Get child entities form meshset containing parent entities.
MoFEMErrorCode getEntitiesByDimAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const int dim, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode setBitRefLevelByDim(const EntityHandle meshset, const int dim, const BitRefLevel bit, int verb=QUIET) const
Set the Bit Ref Level By Dim object.
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childrens.
#define MOFEM_LOG(channel, severity)
Log.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
FTensor::Index< 'l', 3 > l
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
MoFEMErrorCode filterEntitiesByRefLevel(const BitRefLevel bit, const BitRefLevel mask, Range &ents, int verb=QUIET) const
filter entities by bit ref level
MoFEMErrorCode writeBitLevel(const BitRefLevel bit, const BitRefLevel mask, const char *file_name, const char *file_type, const char *options, const bool check_for_empty=true) const
Write bit ref level to file.
virtual moab::Interface & get_moab()=0
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
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.
Mesh refinement interface.
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
MoFEMErrorCode refineTris(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine triangles in the meshset
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.