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.