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.