v0.14.0
Loading...
Searching...
No Matches
mesh_refine_atom.cpp
Go to the documentation of this file.
1
2
3#include <MoFEM.hpp>
4
5using namespace MoFEM;
6
7static char help[] = "testing mesh refinement algorithm\n\n";
8
9int main(int argc, char *argv[]) {
10
11 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
12
13 try {
14
15 PetscBool flg = PETSC_TRUE;
16 char mesh_file_name[255];
17 ierr = PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
18 255, &flg);
20 if (flg != PETSC_TRUE) {
21 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
22 "*** ERROR -my_file (MESH FILE NEEDED)");
23 }
24
25 moab::Core mb_instance;
26 moab::Interface &moab = mb_instance;
27
28 const char *option;
29 option = "";
30 CHKERR moab.load_file(mesh_file_name, 0, option);
31
32 MoFEM::Core core(moab);
33 MoFEM::Interface &m_field = core;
34
35 auto get_dim = [&]() {
36 int dim;
37 int nb_ents_3d;
38 CHKERR m_field.get_moab().get_number_entities_by_dimension(
39 0, 3, nb_ents_3d, true);
40 if (nb_ents_3d > 0) {
41 dim = 3;
42 } else {
43 int nb_ents_2d;
44 CHKERR m_field.get_moab().get_number_entities_by_dimension(
45 0, 2, nb_ents_2d, true);
46 if (nb_ents_2d > 0) {
47 dim = 2;
48 } else {
49 dim = 1;
50 }
51 }
52 return dim;
53 };
54
55 auto dim = get_dim();
56
57 auto refine = m_field.getInterface<MeshRefinement>();
58
59 BitRefLevel bit_level0;
60 bit_level0.set(0);
61 if (dim == 3 || dim == 2) {
62 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
63 0, dim, bit_level0);
64 } else {
65 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
66 "Dimension not handled by test");
67 }
68
69 BitRefLevel bit_level1;
70 bit_level1.set(1);
71
72 auto refine_edges = [&](auto bit0, auto bit) {
74 auto meshset_ptr = get_temp_meshset_ptr(moab);
75 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
76 bit0, BitRefLevel().set(), *meshset_ptr);
77
78 // random mesh refinement
79 auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
80 Range edges_to_refine;
81 CHKERR moab.get_entities_by_type(*meshset_ptr, MBEDGE, edges_to_refine);
82 int ii = 0;
83 for (Range::iterator eit = edges_to_refine.begin();
84 eit != edges_to_refine.end(); eit++, ii++) {
85 int numb = ii % 2;
86 if (numb == 0) {
87 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*eit, 1);
88 }
89 }
90 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit,
91 false, QUIET, 10000);
92 if (dim == 3) {
93 CHKERR refine->refineTets(*meshset_ptr, bit, QUIET, true);
94 } else if (dim == 2) {
95 CHKERR refine->refineTris(*meshset_ptr, bit, QUIET, true);
96 } else {
97 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
98 "Dimension not handled by test");
99 }
101 };
102
103 auto refine_ents_hanging_nodes = [&](auto bit0, auto bit) {
105 auto meshset_ptr = get_temp_meshset_ptr(moab);
106 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
107 bit0, BitRefLevel().set(), *meshset_ptr);
108
109 // random mesh refinement
110 auto meshset_ref_edges_ptr = get_temp_meshset_ptr(moab);
111 Range ents_dim;
112 CHKERR moab.get_entities_by_dimension(*meshset_ptr, dim, ents_dim);
113 int ii = 0;
114 for (Range::iterator eit = ents_dim.begin(); eit != ents_dim.end();
115 eit++, ii++) {
116 int numb = ii % 2;
117 if (numb == 0) {
118 std::vector<EntityHandle> adj_ents;
119 CHKERR moab.get_adjacencies(&*eit, 1, 1, false, adj_ents);
120 CHKERR moab.add_entities(*meshset_ref_edges_ptr, &*adj_ents.begin(),
121 adj_ents.size());
122 }
123 }
124
125 CHKERR refine->addVerticesInTheMiddleOfEdges(*meshset_ref_edges_ptr, bit,
126 false, QUIET, 10000);
127 if (dim == 3) {
128 CHKERR refine->refineTetsHangingNodes(*meshset_ptr, bit, QUIET, true);
129 } else if (dim == 2) {
130 CHKERR refine->refineTrisHangingNodes(*meshset_ptr, bit, QUIET, true);
131 } else {
132 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
133 "Dimension not handled by test");
134 }
136 };
137
138 auto save_blessed_field = [&](auto bit) {
140
141 std::ofstream myfile;
142 myfile.open("mesh_refine.txt");
143
144 auto out_meshset_tet_ptr = get_temp_meshset_ptr(moab);
145 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
146 bit, BitRefLevel().set(), dim, *out_meshset_tet_ptr);
147 Range tets;
148 CHKERR moab.get_entities_by_handle(*out_meshset_tet_ptr, tets);
149 {
150 int ii = 0;
151 for (Range::iterator tit = tets.begin(); tit != tets.end(); tit++) {
152 int num_nodes;
153 const EntityHandle *conn;
154 CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
155
156 for (int nn = 0; nn < num_nodes; nn++) {
157 // cout << conn[nn] << " ";
158 myfile << conn[nn] << " ";
159 }
160 // cout << std::endl;
161 myfile << std::endl;
162 if (ii > 25)
163 break;
164 }
165 }
166
167 myfile.close();
168
170 };
171
172 auto save_vtk = [&](auto bit) {
174 auto out_meshset_tet_ptr = get_temp_meshset_ptr(moab);
175 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByDimAndRefLevel(
176 bit, BitRefLevel().set(), dim, *out_meshset_tet_ptr);
177 CHKERR moab.write_file("out_mesh_refine.vtk", "VTK", "",
178 out_meshset_tet_ptr->get_ptr(), 1);
180 };
181
182 CHKERR refine_edges(BitRefLevel().set(0), BitRefLevel().set(1));
183 CHKERR refine_ents_hanging_nodes(BitRefLevel().set(1),
184 BitRefLevel().set(2));
185 CHKERR save_blessed_field(BitRefLevel().set(2));
186 CHKERR save_vtk(BitRefLevel().set(2));
187 }
189
191
192 return 0;
193}
int main()
@ QUIET
#define CATCH_ERRORS
Catch errors.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_INVALID_DATA
Definition definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
const int dim
auto bit
set bit
static char help[]
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=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:112
Deprecated interface functions.
Mesh refinement interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.