Branch data Line data Source code
1 : : #include "meshkit/CAMALTriAdvance.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 "CMLTriAdvance.hpp"
11 : : #include "RefEntity.hpp"
12 : :
13 : : #ifdef PARALLEL
14 : : #ifdef HAVE_PARALLEL_MOAB
15 : : #include "moab/ParallelComm.hpp"
16 : : #endif
17 : : #endif
18 : :
19 : : #include <vector>
20 : :
21 : : const bool debug_camaltriadv = false;
22 : :
23 : : namespace MeshKit
24 : : {
25 : :
26 : : moab::EntityType CAMALTriAdvance::meshTps[] = {moab::MBVERTEX, moab::MBTRI, moab::MBMAXTYPE};
27 : :
28 : 7 : CAMALTriAdvance::CAMALTriAdvance(MKCore *mk_core, const MEntVector &me_vec)
29 : 7 : : MeshScheme(mk_core, me_vec)
30 : : {
31 : 7 : }
32 : :
33 : 21 : CAMALTriAdvance::~CAMALTriAdvance()
34 : : {
35 [ - + ]: 14 : }
36 : :
37 : 7 : void CAMALTriAdvance::setup_this()
38 : : {
39 : : // just call setup_boundary, since that's the only constraint we have
40 : 7 : setup_boundary();
41 : 7 : }
42 : :
43 : 7 : void CAMALTriAdvance::execute_this()
44 : : {
45 [ + - ][ + - ]: 34 : for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
[ + + ]
46 : : // make a me, for convenience
47 [ + - ]: 27 : ModelEnt *me = (*sit).first;
48 : :
49 : : // create a surface evaluator for this modelent, and a size evaluator
50 [ + - ]: 27 : CAMALSurfEval cse(me);
51 : :
52 [ + - ][ + - ]: 27 : SizingFunction * sz = mk_core()->sizing_function(me->sizing_function_index());
[ + - ]
53 [ + - ][ + + ]: 27 : if (!sz->variable())
54 : 26 : sz = NULL; // no function variable
55 : :
56 [ + - ][ + - ]: 54 : CAMALSizeEval mesize(me->mesh_interval_size(), sz);
57 : : // make sure the size isn't negative
58 [ + - ][ - + ]: 27 : if (mesize.get_size() == -1) mesize.set_size(1.0);
[ # # ]
59 : :
60 : : // assemble bounding mesh
61 [ + - ]: 54 : std::vector<moab::EntityHandle> bdy;
62 [ + - ][ + - ]: 54 : std::vector<int> group_sizes, bdy_ids;
63 [ + - ]: 27 : me->boundary(0, bdy, NULL, &group_sizes);
64 : :
65 : : // convert the handles to integers for input to TriAdv
66 [ + - ]: 54 : moab::Range bdy_vrange;
67 [ + - ]: 54 : std::vector<double> coords;
68 [ + - ]: 27 : me->get_indexed_connect_coords(bdy, NULL, NULL, bdy_ids, coords, &bdy_vrange);
69 : :
70 : : // now construct the CAMAL mesher, and pass it initial conditions
71 [ + - ]: 54 : CMLTriAdvance triadv(&cse, &mesize);
72 [ + - ][ + - ]: 27 : bool success = triadv.set_boundary_mesh(bdy_vrange.size(), &coords[0], group_sizes.size(), &group_sizes[0], &bdy_ids[0]);
[ + - ][ + - ]
[ + - ]
73 [ - + ]: 27 : if (!success)
74 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Trouble setting boundary mesh.");
75 : :
76 : : // ok, now generate the mesh
77 : : int num_pts, num_tris;
78 [ + - ]: 27 : success = triadv.generate_mesh(num_pts, num_tris);
79 [ - + ]: 27 : if (!success) {
80 : : if (debug_camaltriadv) print_debug(me, coords, bdy_vrange, group_sizes, bdy_ids);
81 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Trouble generating tri mesh.");
82 : : }
83 : :
84 [ + - ]: 27 : moab::Range &new_ents = (*sit).second;
85 : : moab::ErrorCode rval;
86 : :
87 : : // resize the coords array, then get the coords of the new points
88 [ + - ][ - + ]: 27 : assert(num_pts >= (int)bdy_vrange.size());
89 [ + - ][ + - ]: 27 : if (num_pts > (int)bdy_vrange.size()) {
90 [ + - ][ + - ]: 27 : coords.resize(3*(num_pts-bdy_vrange.size()));
91 [ + - ][ + - ]: 27 : unsigned int pts_returned = triadv.get_points_buf(coords.size(), &coords[0], bdy_vrange.size());
[ + - ]
92 [ + - ][ - + ]: 27 : if (pts_returned != num_pts-bdy_vrange.size())
93 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Number of new points returned from TriAdv doesn't agree with previous value output.");
94 : :
95 : : // create the new vertices' entities on the face
96 [ + - ][ + - ]: 27 : rval = mk_core()->moab_instance()->create_vertices(&coords[0], pts_returned, new_ents);
[ + - ][ + - ]
97 [ - + ][ # # ]: 27 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
98 : : }
99 : :
100 : : // for tris, pre-allocate connectivity
101 : : moab::ReadUtilIface *iface;
102 [ + - ][ + - ]: 27 : rval = mk_core()-> moab_instance() -> query_interface(iface);
[ + - ]
103 [ - + ][ # # ]: 27 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
104 : :
105 : : //create the tris, get a direct ptr to connectivity
106 : : moab::EntityHandle starth, *connect;
107 [ + - ]: 27 : rval = iface->get_element_connect(num_tris, 3, moab::MBTRI, 1, starth, connect);
108 [ - + ][ # # ]: 27 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
109 : :
110 : : // read connectivity directly into that array, as int's
111 : 27 : int *connecti = (int*) connect;
112 [ + - ]: 27 : int tris_returned = triadv.get_tris_buf(3*num_tris, connecti);
113 [ - + ]: 27 : if (tris_returned != num_tris)
114 [ # # ][ # # ]: 0 : ECERRCHK(MK_FAILURE, "Number of new tris returned from TriAdv doesn't agree with previous value output.");
115 : :
116 : : // put vertex handles into an indexible array
117 [ + - ]: 54 : std::vector<moab::EntityHandle> bdy_vvec;
118 [ + - ][ + - ]: 27 : std::copy(bdy_vrange.begin(), bdy_vrange.end(), std::back_inserter(bdy_vvec));
[ + - ][ + - ]
119 [ + - ][ + - ]: 27 : std::copy(new_ents.begin(), new_ents.end(), std::back_inserter(bdy_vvec));
[ + - ][ + - ]
120 : :
121 : : // now convert vertex indices into handles in-place, working from the back
122 [ + + ]: 22392 : for (int i = 3*num_tris-1; i >= 0; i--) {
123 [ + - ][ - + ]: 22365 : assert(connecti[i] >= 0 && connecti[i] < (int)bdy_vvec.size());
124 [ + - ]: 22365 : connect[i] = bdy_vvec[connecti[i]];
125 : : }
126 : :
127 : : // put new tris into new entity range, then commit the mesh
128 [ + - ]: 27 : new_ents.insert(starth, starth+num_tris-1);
129 [ + - ]: 27 : me->commit_mesh(new_ents, COMPLETE_MESH);
130 : 27 : }
131 : 7 : }
132 : :
133 : 0 : void CAMALTriAdvance::print_debug(ModelEnt *me, std::vector<double> &coords,
134 : : moab::Range &bdy_vrange, std::vector<int> &group_sizes,
135 : : std::vector<int> &bdy_ids)
136 : : {
137 [ # # ]: 0 : std::cout << "Surface_bounadry_mesh: mesh_size = "
138 [ # # ][ # # ]: 0 : << me->mesh_interval_size() << std::endl;
[ # # ]
139 : :
140 [ # # ][ # # ]: 0 : for (int i = 0; i < (int) bdy_vrange.size(); i++) {
141 [ # # ][ # # ]: 0 : std::cout << coords[3 * i] << " " << coords[3 * i + 1]
[ # # ][ # # ]
[ # # ]
142 [ # # ][ # # ]: 0 : << " " << coords[3 * i + 2] << std::endl;
[ # # ][ # # ]
143 : : }
144 : :
145 [ # # ][ # # ]: 0 : std::cout << "bdy_vertex_size:" << bdy_vrange.size()
[ # # ]
146 [ # # ][ # # ]: 0 : << ", group_size:" << group_sizes.size() << std::endl;
[ # # ]
147 : :
148 : 0 : int index = 0;
149 [ # # ]: 0 : for (int i = 0; i < (int) group_sizes.size(); i++) {
150 [ # # ]: 0 : int g_size = group_sizes[i];
151 [ # # ][ # # ]: 0 : std::cout << "boundary_order_group" << i + 1 << ", group_size="
[ # # ]
152 [ # # ][ # # ]: 0 : << g_size << std::endl;
153 [ # # ]: 0 : for (int j = 0; j < g_size; j++) {
154 [ # # ][ # # ]: 0 : std::cout << bdy_ids[index + j] << " ";
[ # # ]
155 : : }
156 [ # # ]: 0 : std::cout << std::endl;
157 : 0 : index += g_size;
158 : : }
159 : :
160 : : moab::ErrorCode rval;
161 : : moab::EntityHandle outset;
162 [ # # ]: 0 : std::string outfile;
163 [ # # ][ # # ]: 0 : std::stringstream ss;
164 : :
165 [ # # ]: 0 : RefEntity* entity = reinterpret_cast<RefEntity*> (me->geom_handle());
166 [ # # ]: 0 : ss << "CAMALTri_boundary_surf";
167 [ # # ][ # # ]: 0 : ss << entity->id();
168 [ # # ]: 0 : ss << "_";
169 : : #ifdef PARALLEL
170 : : #ifdef HAVE_PARALLEL_MOAB
171 : : moab::ParallelComm* pcomm = moab::ParallelComm::get_pcomm(mk_core()->moab_instance(), 0);
172 : : ss << "proc";
173 : : ss << pcomm->proc_config().proc_rank();
174 : : #endif
175 : : #endif
176 [ # # ]: 0 : ss >> outfile;
177 [ # # ]: 0 : outfile += ".vtk";
178 [ # # ][ # # ]: 0 : rval = mk_core()->moab_instance()->create_meshset(0, outset);
[ # # ]
179 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
180 [ # # ][ # # ]: 0 : rval = mk_core()->moab_instance()->add_entities(outset, bdy_vrange);
[ # # ]
181 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
182 [ # # ][ # # ]: 0 : rval = mk_core()->moab_instance()->write_mesh(outfile.c_str(), &outset, 1);
[ # # ]
183 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
184 [ # # ][ # # ]: 0 : std::cout << outfile.c_str() << " is saved." << std::endl;
[ # # ]
185 : 0 : }
186 [ + - ][ + - ]: 156 : } // namespace MeshKit
|
