v0.15.0
Loading...
Searching...
No Matches
mofem_part.cpp
Go to the documentation of this file.
1/** \file reading_med.cpp
2
3 \brief Partition mesh and configuring blocksets
4
5*/
6
7
8
9#include <MoFEM.hpp>
10using namespace MoFEM;
11
12static char help[] = "-file_name <filename> : mesh file name\n"
13 "-output_file <filename> : output mesh file name\n"
14 "-nparts <int> : number of parts\n"
15 "-dim <int> : entities dim\n"
16 "-adj_dim <int> : adjacency dim\n"
17 "-my_create_lower_dim_ents <bool> : "
18 "if true create lower dimension entireties\n"
19 "-block_tags <bool> : only block and meshsests tags \n\n";
20
21int main(int argc, char *argv[]) {
22
23 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
24
25 try {
26
27 moab::Core mb_instance;
28 moab::Interface &moab = mb_instance;
29 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
30 if (pcomm == NULL)
31 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
32
33 // global variables
34 char mesh_file_name[255];
35 char mesh_out_file[255] = "out.h5m";
36 PetscBool flg_file = PETSC_FALSE;
37 PetscBool flg_n_part = PETSC_FALSE;
38 PetscBool flg_part = PETSC_FALSE;
39 PetscBool only_tags = PETSC_FALSE;
40 PetscInt n_partas = 1;
41 PetscBool create_lower_dim_ents = PETSC_TRUE;
42 PetscInt dim = 3;
43 PetscInt adj_dim = 0;
44
45 PetscOptionsBegin(PETSC_COMM_WORLD, "", "none", "none");
46
47 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
48 mesh_file_name, 255, &flg_file);
49 if (flg_file != PETSC_TRUE)
50 CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
51 mesh_file_name, 255, &flg_file);
52 if (flg_file != PETSC_TRUE)
53 SETERRQ(PETSC_COMM_SELF, 1,
54 "*** ERROR -my_file (-file_name) (mesh file needed)");
55
56 const char *option;
57 option = "";
58 CHKERR moab.load_file(mesh_file_name, 0, option);
59
60 CHKERR PetscOptionsString("-output_file", "output mesh file name", "",
61 "out.h5m", mesh_out_file, 255, PETSC_NULLPTR);
62 CHKERR PetscOptionsInt("-my_nparts", "number of parts", "", n_partas,
63 &n_partas, &flg_n_part);
64 CHKERR PetscOptionsInt("-nparts", "number of parts", "", n_partas,
65 &n_partas, &flg_part);
66
67 if (!flg_n_part && !flg_part)
68 SETERRQ(PETSC_COMM_SELF, 1,
69 "*** ERROR partitioning number not given (-nparts)");
70
71 auto get_nb_ents_by_dim = [&](const int dim) {
72 int nb;
73 CHKERR moab.get_number_entities_by_dimension(0, dim, nb);
74 return nb;
75 };
76 for (; dim >= 0; dim--) {
77 if (get_nb_ents_by_dim(dim))
78 break;
79 }
80
81 if (!dim)
82 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
83 "Dimension of entities to partition not found");
84
85 CHKERR PetscOptionsInt("-dim", "entities dim", "", dim, &dim, PETSC_NULLPTR);
86 // adj_dim = dim - 1;
87 CHKERR PetscOptionsInt("-adj_dim", "adjacency dim", "", adj_dim, &adj_dim,
88 PETSC_NULLPTR);
89 CHKERR PetscOptionsBool(
90 "-my_create_lower_dim_ents", "if true create lower dimension entireties",
91 "", create_lower_dim_ents, &create_lower_dim_ents, PETSC_NULLPTR);
92 CHKERR PetscOptionsBool("-block_tags", "only block and meshsests tags", "",
93 only_tags, &only_tags, PETSC_NULLPTR);
94
95 PetscOptionsEnd();
96
97 auto get_dim = [](const Range &ents) -> int {
98 for (int d : {3, 2, 1})
99 if (ents.num_of_dimension(d))
100 return d;
101 return 0;
102 };
103
104 // Create MoFEM database
105 MoFEM::Core core(moab);
106 MoFEM::Interface &m_field = core;
107
108 auto meshsets_interface_ptr = m_field.getInterface<MeshsetsManager>();
109 CHKERR meshsets_interface_ptr->setMeshsetFromFile();
110
111 MOFEM_LOG_CHANNEL("WORLD");
112 MOFEM_LOG_TAG("WORLD", "mofem_part");
113 int min_dim = 3;
114 for (auto cit = meshsets_interface_ptr->getBegin();
115 cit != meshsets_interface_ptr->getEnd(); cit++) {
116 Range ents;
117 CHKERR m_field.get_moab().get_entities_by_handle(cit->getMeshset(), ents,
118 true);
119 min_dim = std::min(min_dim, get_dim(ents));
120 MOFEM_LOG("WORLD", Sev::inform) << *cit;
121 }
122
123 int g_dim; // global dimension
124 MPI_Allreduce(&min_dim, &g_dim, 1, MPI_INT, MPI_MIN, m_field.get_comm());
125 if (g_dim < adj_dim) {
126 MOFEM_LOG("WORLD", Sev::warning)
127 << "The minimum meshsets dimension is = " << min_dim;
128 }
129 if (adj_dim >= dim) {
130 MOFEM_LOG("WORLD", Sev::warning)
131 << "The -adj_dim >= dim, adj_dim = " << adj_dim << " dim = " << dim;
132 }
133 MOFEM_LOG_CHANNEL("WORLD");
134
135 {
136 Range ents_dim;
137 CHKERR moab.get_entities_by_dimension(0, dim, ents_dim, false);
138 if (create_lower_dim_ents) {
139 if (dim == 3) {
140 Range faces;
141 CHKERR moab.get_adjacencies(ents_dim, 2, true, faces,
142 moab::Interface::UNION);
143 }
144 if (dim > 2) {
145 Range edges;
146 CHKERR moab.get_adjacencies(ents_dim, 1, true, edges,
147 moab::Interface::UNION);
148 }
149 }
150 CommInterface *comm_interafce_ptr = m_field.getInterface<CommInterface>();
151 CHKERR comm_interafce_ptr->partitionMesh(
152 ents_dim, dim, adj_dim, n_partas, nullptr, nullptr, nullptr, VERBOSE);
153 }
154
155 auto get_tag_list = [&]() {
156 std::vector<Tag> tags_list;
157 auto meshsets_mng = m_field.getInterface<MeshsetsManager>();
158 auto &list = meshsets_mng->getMeshsetsMultindex();
159 for (auto &m : list) {
160 auto meshset = m.getMeshset();
161 std::vector<Tag> tmp_tags_list;
162 CHKERR m_field.get_moab().tag_get_tags_on_entity(meshset,
163 tmp_tags_list);
164 for (auto t : tmp_tags_list) {
165 tags_list.push_back(t);
166 }
167 }
168
169 return tags_list;
170 };
171
172 EntityHandle root_mesh = 0;
173 if (m_field.get_comm_rank() == 0) {
174 if (only_tags) {
175 auto tags_list = get_tag_list();
176 std::sort(tags_list.begin(), tags_list.end());
177 auto new_end = std::unique(tags_list.begin(), tags_list.end());
178 tags_list.resize(std::distance(tags_list.begin(), new_end));
179 tags_list.push_back(pcomm->part_tag());
180 // tags_list.push_back(m_field.get_moab().globalId_tag());
181 CHKERR moab.write_file(mesh_out_file, "MOAB", "", &root_mesh, 1,
182 &*tags_list.begin(), tags_list.size());
183 } else {
184 CHKERR moab.write_file(mesh_out_file, "MOAB", "");
185 }
186 MOFEM_LOG("WORLD", Sev::inform) << "Wrote file " << mesh_out_file;
187 }
188 }
190
192 CHKERRQ(ierr);
193
194 return 0;
195}
int main()
@ VERBOSE
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
MoFEMErrorCode partitionMesh(const Range &ents, const int dim, const int adj_dim, const int n_parts, Tag *th_vertex_weights=nullptr, Tag *th_edge_weights=nullptr, Tag *th_part_weights=nullptr, int verb=VERBOSE, const bool debug=false)
Set partition tag to each finite element in the problem.
char mesh_file_name[255]
static char help[]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
constexpr double t
plate stiffness
Definition plate.cpp:58
FTensor::Index< 'm', 3 > m
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
Interface for managing meshsets containing materials and boundary conditions.
CubitMeshSet_multiIndex & getMeshsetsMultindex()
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.