1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208 | #include "meshkit/CAMALTetMesher.hpp"
#include "meshkit/MeshScheme.hpp"
#include "meshkit/ModelEnt.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
#include "meshkit/RegisterMeshOp.hpp"
#include "CAMALSurfEval.hpp"
#include "CAMALSizeEval.hpp"
#include "CMLTetMesher.hpp"
#include "RefEntity.hpp"
#include <vector>
#ifdef PARALLEL
#ifdef HAVE_PARALLEL_MOAB
#include "moab/ParallelComm.hpp"
#endif
#endif
const bool debug_camaltet = false;
namespace MeshKit
{
moab::EntityType CAMALTetMesher::meshTps[] = {moab::MBVERTEX, moab::MBTET, moab::MBMAXTYPE};
CAMALTetMesher::CAMALTetMesher(MKCore *mk_core, const MEntVector &me_vec)<--- Member variable 'CAMALTetMesher::cmlTetMesher' is not initialized in the constructor.
: MeshScheme(mk_core, me_vec)
{
}
CAMALTetMesher::~CAMALTetMesher()
{
}
MeshOp *CAMALTetMesher::get_tri_mesher()
{
MeshOpProxy * mproxy = mk_core()->meshop_proxy("CAMALTriAdvance") ;
return mk_core()->construct_meshop(mproxy);
}
void CAMALTetMesher::setup_this()
{
MeshOp *tri_mesher = NULL;
std::vector<MeshOp*> meshops;<--- Unused variable: meshops
MEntVector surfs;
for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {<--- Prefer prefix ++/-- operators for non-primitive types.
ModelEnt *me = (*sit).first;
if (me->dimension() != 3) throw Error(MK_BAD_INPUT, "Tet mesher assigned to an entity with dimension != 3.");
// get the bounding faces
surfs.clear();
me->get_adjacencies(2, surfs);
// check the mesh scheme on each one; if there is one, verify it can generate tris; if there
// isn't one, make one
bool inserted = false;
for (MEntVector::iterator vit = surfs.begin(); vit != surfs.end(); vit++) {<--- Prefer prefix ++/-- operators for non-primitive types.
if ((*vit)->is_meshops_list_empty()) {
// get a tri mesher if we haven't already
if (!tri_mesher) tri_mesher = get_tri_mesher();
// add this surface to it, and if first for the volume, make sure it's added upstream
tri_mesher->add_modelent(*vit);
if (!inserted) {
mk_core()->insert_node(tri_mesher, this);
inserted = true;
}
} // if no meshops
} // over surfs
} // over vols
}
void CAMALTetMesher::execute_this()
{
for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {<--- Prefer prefix ++/-- operators for non-primitive types.
// make a me, for convenience
ModelEnt *me = (*sit).first;
// assemble bounding mesh
std::vector<moab::EntityHandle> bdy;
std::vector<int> group_sizes, senses, bdy_ids;
me->boundary(2, bdy, &senses, &group_sizes);
// get connectivity
// convert the handles to integers for input to CMLTetMesher
moab::Range bdy_vrange;
std::vector<double> coords;
me->get_indexed_connect_coords(bdy, &senses, NULL, bdy_ids, coords, &bdy_vrange);
CMLTetMesher tet_mesher;
bool success = tet_mesher.set_boundary_mesh(bdy_vrange.size(), &coords[0], bdy.size(), &bdy_ids[0]);
if (!success)
ECERRCHK(MK_FAILURE, "Failed setting boundary mesh");
// generate the mesh
int num_tets, num_pts;
success = tet_mesher.generate_mesh(num_pts, num_tets);
if (!success) {
if (debug_camaltet) print_debug(me, coords, bdy_vrange, bdy, group_sizes, bdy_ids);
ECERRCHK(MK_FAILURE, "Trouble generating tet mesh.");
}
moab::Range &new_ents = (*sit).second;
moab::ErrorCode rval;
// resize the coords array, then get the coords of the new points
assert(num_pts >= (int)bdy_vrange.size());
if (num_pts > (int)bdy_vrange.size()) {
coords.resize(3*(num_pts-bdy_vrange.size()));
int pts_returned = tet_mesher.get_points_buf(coords.size(), &coords[0], bdy_vrange.size());
if (pts_returned != num_pts-(int)bdy_vrange.size())
ECERRCHK(MK_FAILURE, "Number of new points returned from TetMesher doesn't agree with previous value output.");
// create the new vertices' entities
rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents);
MBERRCHK(rval, mk_core()->moab_instance());
}
// for tets, pre-allocate connectivity
moab::ReadUtilIface *iface;
rval = mk_core()-> moab_instance() -> query_interface(iface);
MBERRCHK(rval, mk_core()->moab_instance());
//create the tris, get a direct ptr to connectivity
moab::EntityHandle starth, *connect;// *tmp_connect;
rval = iface->get_element_connect(num_tets, 4, moab::MBTET, 1, starth, connect);
MBERRCHK(rval, mk_core()->moab_instance());
// read connectivity directly into that array, as int's
int *connecti = (int*) connect;
int tets_returned = tet_mesher.get_tets_buf(4*num_tets, connecti);
if (tets_returned != num_tets)
ECERRCHK(MK_FAILURE, "Number of new tets returned from TetMesher doesn't agree with previous value output.");
// put vertex handles into an indexible array
std::vector<moab::EntityHandle> bdy_vvec;
std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
// now convert vertex indices into handles in-place, working from the back
for (int i = 4*num_tets-1; i >= 0; i--) {
assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size());
connect[i] = bdy_vvec[connecti[i]];
}
// put new tris into new entity range, then commit the mesh
new_ents.insert(starth, starth+num_tets-1);
me->commit_mesh(new_ents, COMPLETE_MESH);
}
}
void CAMALTetMesher::print_debug(ModelEnt *me, std::vector<double> &coords,
moab::Range &bdy_vrange, std::vector<moab::EntityHandle> &bdy,
std::vector<int> &group_sizes, std::vector<int> &bdy_ids)
{
std::cout << "Volume_boundary_mesh: mesh_size = "
<< me->mesh_interval_size() << std::endl;
for (int i = 0; i < (int)bdy_vrange.size(); i++) {
std::cout << coords[3 * i] << " " << coords[3 * i + 1]
<< " " << coords[3 * i + 2] << std::endl;
}
std::cout << "bdy_vertex_size:" << bdy_vrange.size()
<< ", group_size:" << group_sizes.size() << std::endl;
int index = 0;
for (size_t i = 0; i < group_sizes.size(); i++) {
int g_size = group_sizes[i];
std::cout << "boundary_order_group" << i + 1 << ", group_size="
<< g_size << std::endl;
for (int j = 0; j < g_size; j++) {
std::cout << bdy_ids[index + j] << " ";
}
std::cout << std::endl;
index += g_size;
}
moab::ErrorCode rval;
moab::EntityHandle outset;
std::string outfile;
std::stringstream ss;
RefEntity* entity = reinterpret_cast<RefEntity*> (me->geom_handle());
ss << "CAMALTet_boundary";
ss << entity->id();
ss << "_";
#ifdef PARALLEL
#ifdef HAVE_PARALLEL_MOAB
moab::ParallelComm* pcomm = moab::ParallelComm::get_pcomm(mk_core()->moab_instance(), 0);
ss << "proc";
ss << pcomm->proc_config().proc_rank();
#endif
#endif
ss >> outfile;
outfile += ".vtk";
rval = mk_core()->moab_instance()->create_meshset(0, outset);
MBERRCHK(rval, mk_core()->moab_instance());
rval = mk_core()->moab_instance()->add_entities(outset, &bdy[0], bdy.size());
MBERRCHK(rval, mk_core()->moab_instance());
rval = mk_core()->moab_instance()->write_mesh(outfile.c_str(), &outset, 1);
MBERRCHK(rval, mk_core()->moab_instance());
std::cout << outfile.c_str() << " is saved." << std::endl;
}
} // namespace MeshKit
|