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