v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
write_med.cpp File Reference

write med files More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Detailed Description

write med files

Definition in file write_med.cpp.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 12 of file write_med.cpp.

12 {
13 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
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";
24 char mesh_file_name[255] = "mesh.h5m";
25 char meshset_to_write[255] = "";
26 PetscBool meshset_to_write_set = PETSC_FALSE;
27 int time_step = 0;
28 // Create MoFEM database
29 MoFEM::Core core(moab, PETSC_COMM_WORLD, VERBOSE);
30 MoFEM::Interface &m_field = core;
31 PetscBool bit_ref_set = PETSC_FALSE;
32 int bit_ref_level = 0;
33
34 CHKERR PetscOptionsBegin(m_field.get_comm(), "", "Read MED tool", "none");
35 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-file_name", mesh_file_name,
36 255, PETSC_NULL);
37 CHKERR PetscOptionsInt("-bit_ref_level", "bit ref level", "", bit_ref_level,
38 &bit_ref_level, &bit_ref_set);
39 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-meshset_to_write",
40 meshset_to_write, 255, &meshset_to_write_set);
41 CHKERR PetscOptionsInt("-med_time_step", "time step", "", time_step,
42 &time_step, PETSC_NULL);
43 CHKERR PetscOptionsString("-output_file", "output mesh file name", "",
44 "out.h5m", mesh_out_file, 255, PETSC_NULL);
45 ierr = PetscOptionsEnd();
46 CHKERRQ(ierr);
47
48 const char *option = "";
49
50 // load mesh
51 CHKERR moab.load_file(mesh_file_name, 0, option);
52
53 // read meshsets
54 CHKERR m_field.getInterface<MeshsetsManager>()->readMeshsets();
55
56 // define range to write
57 Range entities_to_write;
58 if (bit_ref_set) {
59 // for writing mesh with crack
60 BitRefLevel bit_level;
61 // write mesh with crack
62 bit_level.set(bit_ref_level);
63
64 // get all entities in bit ref level
65 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
66 bit_level, BitRefLevel().set(), entities_to_write);
67
68 // remove prism entities
69 entities_to_write = subtract(entities_to_write,
70 entities_to_write.subset_by_type(MBPRISM));
71 // remove quad entities
72 entities_to_write =
73 subtract(entities_to_write, entities_to_write.subset_by_type(MBQUAD));
74 MOFEM_LOG("MEDWORLD", Sev::inform)
75 << "Total number of entities in bit ref level " << bit_ref_level
76 << ": " << entities_to_write.size() << std::endl;
77
78 MOFEM_LOG("MEDWORLD", Sev::inform)
79 << "Removing interior edges and face entities"<< std::endl;
80 // get all subset ranges
81 Range entities_3d = entities_to_write.subset_by_dimension(3);
82 Range entities_2d = entities_to_write.subset_by_dimension(2);
83 Range entities_1d = entities_to_write.subset_by_dimension(1);
84
85 // find skin
86 Range skin_2d;
87 Range skin_1d;
88 Skinner skin(&moab);
89 CHKERR skin.find_skin(0, entities_3d, false, skin_2d);
90 CHKERR moab.get_adjacencies(skin_2d, 1, false, skin_1d,
91 moab::Interface::UNION);
92
93 // remove interior entities from entities to write
94 entities_2d = subtract(entities_2d, skin_2d);
95 entities_to_write = subtract(entities_to_write, entities_2d);
96
97 entities_1d = subtract(entities_1d, skin_1d);
98 entities_to_write = subtract(entities_to_write, entities_1d);
99
100 MOFEM_LOG("MEDWORLD", Sev::inform)
101 << "Removing nodes with no connectivity"<< std::endl;
102 // find nodes with no connectivity
103 Range nodes;
104 //moab.get_entities_by_type(0, MBVERTEX, nodes, true);
105 nodes = entities_to_write.subset_by_type(MBVERTEX);
106 // loop to find nodes with no adjacencies
107 Range free_nodes;
108 for (auto node : nodes) {
109 Range adj_edges;
110 EntityHandle node_handle = node;
111 CHKERR moab.get_adjacencies(&node_handle, 1, 1, false, adj_edges,
112 moab::Interface::UNION);
113 if (adj_edges.empty()) {
114 free_nodes.insert(node);
115 }
116 }
117 // get coordinates of nodes with no connectivity
118 for (auto node : free_nodes) {
119 double coords[3];
120 CHKERR moab.get_coords(&node, 1, coords);
121 }
122
123 // remove nodes with no connectivity
124 entities_to_write = subtract(entities_to_write, free_nodes);
125
126 // Check orientation of crack surfaces
127 // Get Meshest 10000 and 10001
128 EntityHandle meshset_10000;
129 EntityHandle meshset_10001;
130
131 CHKERR m_field.getInterface<MeshsetsManager>()->getMeshset(10000, SIDESET,
132 meshset_10000);
133 CHKERR m_field.getInterface<MeshsetsManager>()->getMeshset(10001, SIDESET,
134 meshset_10001);
135
136 Range entities_10000;
137 Range entities_10001;
138 moab.get_entities_by_handle(meshset_10000, entities_10000);
139 moab.get_entities_by_handle(meshset_10001, entities_10001);
140
141 auto check_face_orientation = [&moab](const Range &faces) {
142 for (auto face : faces) {
143 // check orientation of face
144 // get adjacent tets
145 Range adj_tets;
146 CHKERR(moab.get_adjacencies(&face, 1, 3, false, adj_tets,
147 moab::Interface::UNION));
148 // remove other entities from adj_tets
149 Range actual_tets;
150 actual_tets = adj_tets.subset_by_type(MBTET);
151
152 int side_number;
153 int sense;
154 int offset;
155 for (auto tet : actual_tets) {
156 CHKERR(moab.side_number(tet, face, side_number, sense, offset));
157 if (sense == -1) {
158 // get nodes of face
159 std::vector<EntityHandle> nodes_face;
160 CHKERR(moab.get_connectivity(&face, 1, nodes_face));
161 std::vector<EntityHandle> face_nodes_new;
162 // change orientation
163 face_nodes_new = {nodes_face[1], nodes_face[0], nodes_face[2]};
164 CHKERR(moab.set_connectivity(face, face_nodes_new.data(), 3));
165 // recheck orientation
166 CHKERR(moab.side_number(tet, face, side_number, sense, offset));
167 if (sense == -1) {
168 MOFEM_LOG("MEDWORLD", Sev::warning)
169 << "Face: " << face << " has wrong orientation";
170 }
171 }
172 }
173 }
174 };
175
176 check_face_orientation(entities_10000);
177 check_face_orientation(entities_10001);
178
179 } else if (meshset_to_write_set) {
180 // loop meshsets to find entities to write
181 auto &meshsets_idx =
182 m_field.getInterface<MeshsetsManager>()->getMeshsetsMultindex();
183
184 for (auto &meshset : meshsets_idx) {
185 auto meshset_name = meshset.getName();
186 auto trim_name = [&](std::string &name) {
187 name.erase(std::remove(name.begin(), name.end(), ' '), name.end());
188 };
189 trim_name(meshset_name);
190
191 if (meshset_to_write == meshset_name) {
192 CHKERR moab.get_entities_by_handle(meshset.getMeshset(),
193 entities_to_write, true);
194 break;
195 }
196 }
197 MOFEM_LOG("MEDWORLD", Sev::inform)
198 << "Wrtiting all entitiies from meshset: " << meshset_to_write
199 << std::endl;
200 } else {
201 // get all entities
202 CHKERR moab.get_entities_by_handle(0, entities_to_write, true);
203 MOFEM_LOG("MEDWORLD", Sev::inform)
204 << "Wrtiting all entitiies" << std::endl;
205 }
206 MOFEM_LOG("MEDWORLD", Sev::inform)
207 << "Number of entities to write: " << entities_to_write.size()
208 << std::endl;
209 // loop to print size by entity type
210 for (EntityType ent_type = MBVERTEX; ent_type < MBMAXTYPE; ent_type++) {
211 Range entities;
212 moab.get_entities_by_type(0, ent_type, entities, true);
213 Range ents_to_write;
214 ents_to_write = intersect(entities, entities_to_write);
215
216 if (ents_to_write.empty())
217 continue;
218
219 MOFEM_LOG("MEDWORLD", Sev::inform)
220 << "Number of entities to write: " << ents_to_write.size()
221 << " of type: " << moab::CN::EntityTypeName(ent_type)
222 << " from total of: " << entities.size() << std::endl;
223 }
224
225 // make shared pointer to range
226 boost::shared_ptr<Range> entities_to_write_ptr =
227 boost::make_shared<Range>(entities_to_write);
228
229 MedInterface *med_interface_ptr;
230 CHKERR m_field.getInterface(med_interface_ptr);
231 CHKERR med_interface_ptr->writeMed(entities_to_write_ptr);
232 }
233 CATCH_ERRORS;
234
235 MoFEM::Core::Finalize();
236
237 return 0;
238}
VERBOSE
@ VERBOSE
Definition definitions.h:222
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
SIDESET
@ SIDESET
Definition definitions.h:160
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:548
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:308
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::MedInterface
Interface for load MED files.
Definition MedInterface.hpp:23
MoFEM::MedInterface::writeMed
MoFEMErrorCode writeMed(boost::shared_ptr< Range > range_ptr=nullptr, int verb=1)
write MED file
Definition MedInterface.cpp:686
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:104
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
help
static char help[]
Definition write_med.cpp:10

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 10 of file write_med.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow