Branch data Line data Source code
1 : : //-----------------------------------C++-------------------------------------//
2 : : // File: algs/SCDMesh.cpp
3 : : // Author: Stuart R. Slattery
4 : : // Wednesday February 2 16:15:16 2011
5 : : // Brief: SCDMesh member definitions
6 : : //---------------------------------------------------------------------------//
7 : :
8 : : #include "meshkit/SCDMesh.hpp"
9 : : #include "meshkit/MKCore.hpp"
10 : : #include "meshkit/ModelEnt.hpp"
11 : : #include "meshkit/RegisterMeshOp.hpp"
12 : :
13 : : #include <vector>
14 : :
15 : : #include "moab/Core.hpp"
16 : : #include "moab/Range.hpp"
17 : : #include "moab/Interface.hpp"
18 : : #include "moab/ScdInterface.hpp"
19 : : #include "moab/OrientedBoxTreeTool.hpp"
20 : : #include "moab/CartVect.hpp"
21 : :
22 : : namespace MeshKit
23 : : {
24 : : //---------------------------------------------------------------------------//
25 : : // Initialize the entity types for SCDMesh
26 : : moab::EntityType SCDMesh_tps[] = {moab::MBVERTEX, moab::MBHEX, moab::MBMAXTYPE};
27 : 40 : const moab::EntityType* SCDMesh::output_types()
28 : 40 : { return SCDMesh_tps; }
29 : :
30 : : //---------------------------------------------------------------------------//
31 : : // Constructor for SCDMesh
32 : 4 : SCDMesh::SCDMesh(MKCore *mk_core, const MEntVector &me_vec)
33 [ + - ][ + - ]: 4 : : MeshScheme(mk_core, me_vec)
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
34 : : {
35 : 4 : boxIncrease = 0.0;
36 : 4 : useMeshGeom = false;
37 : :
38 : : // set bounding box size tag
39 : 4 : double bb_default[6] = { 0., 0., 0., 0., 0., 0. };
40 [ + - ]: 4 : rval = mk_core->moab_instance()->tag_get_handle("BOUNDING_BOX_SIZE",
41 : : 6, moab::MB_TYPE_DOUBLE,
42 [ + - ]: 4 : bb_tag, moab::MB_TAG_SPARSE|moab::MB_TAG_CREAT, bb_default);
43 [ - + ][ # # ]: 4 : if (moab::MB_SUCCESS != rval && moab::MB_ALREADY_ALLOCATED != rval)
44 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core->moab_instance());
[ # # ][ # # ]
[ # # ]
45 : 4 : }
46 : :
47 : : //---------------------------------------------------------------------------//
48 : : // setup function
49 : 4 : void SCDMesh::setup_this()
50 : : {
51 : : // if cfMesh chosen, check the grid entries to make sure they're correct
52 [ + - ]: 4 : if (gridType == 0) {
53 [ - + ]: 4 : if (coarse_i != fine_i.size()) {
54 : : throw Error(MK_INCOMPLETE_MESH_SPECIFICATION,
55 [ # # ]: 0 : "Number of coarse i divisions not equal to the number of fine i divisions.");
56 : : }
57 [ - + ]: 4 : if (coarse_j != fine_j.size()) {
58 : : throw Error(MK_INCOMPLETE_MESH_SPECIFICATION,
59 [ # # ]: 0 : "Number of coarse j divisions not equal to the number of fine j divisions.");
60 : : }
61 [ - + ]: 4 : if (coarse_k != fine_k.size()) {
62 : : throw Error(MK_INCOMPLETE_MESH_SPECIFICATION,
63 [ # # ]: 0 : "Number of coarse k divisions not equal to the number of fine k divisions.");
64 : : }
65 : : }
66 : :
67 : : // do setup for the model entities
68 [ + - ][ + - ]: 10 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) {
[ + + ]
69 [ + - ]: 6 : ModelEnt *me = mit->first;
70 : :
71 : : // check if the entity is already meshed
72 [ + - ][ + - ]: 6 : if (me->get_meshed_state() >= COMPLETE_MESH || me->mesh_intervals() > 0) continue;
[ + - ][ - + ]
[ - + ]
73 : : }
74 : 4 : }
75 : :
76 : : //---------------------------------------------------------------------------//
77 : : // Structured cartesian mesh generation execution
78 : 4 : void SCDMesh::execute_this()
79 : : {
80 [ + - ][ + - ]: 7 : for (MEntSelection::iterator mit = mentSelection.begin(); mit != mentSelection.end(); mit++) {
[ + + ]
81 [ + - ]: 6 : ModelEnt *me = mit->first;
82 : :
83 : : /* Generate the bounding box */
84 : :
85 : : // Cartesian case
86 [ + - ]: 6 : if (axisType == 0) {
87 [ + + ]: 6 : if (geomType == 0) {
88 [ + - ]: 3 : set_cart_box_all();
89 : : }
90 [ + - ]: 3 : else if (geomType == 1) {
91 [ + - ]: 3 : set_cart_box_individual(me);
92 : : }
93 [ + - ]: 6 : create_cart_edges();
94 : : }
95 : :
96 : : /* Generate the mesh */
97 : :
98 : : // create mesh vertex coordinates
99 [ + - ]: 6 : create_vertex_coords();
100 : :
101 : : // full representation case
102 [ + + ]: 6 : if (interfaceType == 0) {
103 [ + - ][ + - ]: 5 : create_full_mesh(mit->second);
104 : : }
105 : :
106 : : // lighweight ScdInterface case
107 [ + + ]: 6 : if (interfaceType == 1) {
108 [ + - ][ + - ]: 1 : create_light_mesh(mit->second);
109 : : }
110 : :
111 : : // stop if mesh is created for whole geometry
112 [ + + ]: 6 : if (geomType == 0) break;
113 : :
114 : : // commit the mesh when finished
115 [ + - ][ + - ]: 3 : me->commit_mesh(mit->second, COMPLETE_MESH);
116 : : }
117 : 4 : }
118 : :
119 : : //---------------------------------------------------------------------------//
120 : : // set the Cartesian bounding box dimension for the entire geometry
121 : 3 : void SCDMesh::set_cart_box_all()
122 : : {
123 [ + - ]: 3 : if (!useMeshGeom) {
124 : 3 : gerr = mk_core()->igeom_instance()->getBoundBox(minCoord[0], minCoord[1], minCoord[2],
125 : 3 : maxCoord[0], maxCoord[1], maxCoord[2]);
126 : 3 : IBERRCHK(gerr, "ERROR: couldn't get geometry bounding box");
127 : : }
128 : :
129 : : // increase box size
130 : : double box_size;
131 [ + + ]: 12 : for (int i = 0; i < 3; i++) {
132 : 9 : box_size = maxCoord[i] - minCoord[i];
133 : 9 : minCoord[i] -= box_size*boxIncrease;
134 : 9 : maxCoord[i] += box_size*boxIncrease;
135 : : }
136 : 3 : }
137 : :
138 : : //---------------------------------------------------------------------------//
139 : : // set the Cartesian bounding box dimension for an individual volume
140 : 3 : void SCDMesh::set_cart_box_individual(ModelEnt *this_me)
141 : : {
142 [ + - ]: 3 : gerr = this_me->igeom_instance()->getBoundBox(minCoord[0], minCoord[1], minCoord[2],
143 [ + - ]: 3 : maxCoord[0], maxCoord[1], maxCoord[2]);
144 [ + - ]: 3 : IBERRCHK(gerr, "ERROR: couldn't get geometry bounding box");
145 : :
146 : : // increase box size
147 : : double box_size;
148 [ + + ]: 12 : for (int i = 0; i < 3; i++) {
149 : 9 : box_size = maxCoord[i] - minCoord[i];
150 : 9 : minCoord[i] -= box_size*boxIncrease;
151 : 9 : maxCoord[i] += box_size*boxIncrease;
152 : : }
153 : :
154 : : // set bounding box size as tag
155 [ + - ]: 3 : moab::EntityHandle meshset = this_me->mesh_handle();
156 : 9 : double box_min_max[6] = { minCoord[0], minCoord[1], minCoord[2],
157 : 9 : maxCoord[0], maxCoord[1], maxCoord[2] };
158 [ + - ][ + - ]: 3 : rval = mk_core()->moab_instance()->tag_set_data(bb_tag, &meshset, 1,
159 [ + - ]: 3 : box_min_max);
160 : 3 : }
161 : :
162 : : //---------------------------------------------------------------------------//
163 : : // set the Cartesian bounding box dimension
164 : 0 : void SCDMesh::set_cart_box_min_max(double* min, double* max, double box_increase)
165 : : {
166 [ # # ]: 0 : for (int i = 0; i < 3; i++) {
167 : 0 : double box_size = max[i] - min[i];
168 : 0 : minCoord[i] = min[i] - box_size*box_increase;
169 : 0 : maxCoord[i] = max[i] + box_size*box_increase;
170 : : }
171 : 0 : }
172 : :
173 : : //---------------------------------------------------------------------------//
174 : : // get the axis cartesian bounding box edges if cfmesh has been chosen
175 : 6 : void SCDMesh::create_cart_edges()
176 : : {
177 : : // generate the vector arrays for cartesian grid.
178 : 6 : i_arr.push_back(minCoord[0]);
179 : 6 : double icrs_size = (maxCoord[0] - minCoord[0]) / coarse_i;
180 : : double ifn_size;
181 : : unsigned int crsi;
182 [ + + ]: 41 : for (crsi = 0; crsi != coarse_i; crsi++) {
183 : 35 : ifn_size = icrs_size / fine_i[crsi];
184 [ + + ]: 135 : for (int fni = 1; fni != fine_i[crsi] + 1; fni++) {
185 [ + - ]: 100 : i_arr.push_back(minCoord[0] + crsi*icrs_size + fni*ifn_size);
186 : : }
187 : : }
188 : :
189 : 6 : j_arr.push_back(minCoord[1]);
190 : 6 : double jcrs_size = (maxCoord[1] - minCoord[1]) / coarse_j;
191 : : double jfn_size;
192 : : unsigned int crsj;
193 [ + + ]: 41 : for (crsj = 0; crsj != coarse_j; crsj++) {
194 : 35 : jfn_size = jcrs_size / fine_j[crsj];
195 [ + + ]: 135 : for (int fnj = 1; fnj != fine_j[crsj] + 1; fnj++) {
196 [ + - ]: 100 : j_arr.push_back(minCoord[1] + crsj*jcrs_size + fnj*jfn_size);
197 : : }
198 : : }
199 : :
200 : 6 : k_arr.push_back(minCoord[2]);
201 : 6 : double kcrs_size = (maxCoord[2] - minCoord[2]) / coarse_k;
202 : : double kfn_size;
203 : : unsigned int crsk;
204 [ + + ]: 41 : for (crsk = 0; crsk != coarse_k; crsk++) {
205 : 35 : kfn_size = kcrs_size / fine_k[crsk];
206 [ + + ]: 135 : for (int fnk = 1; fnk != fine_k[crsk] + 1; fnk++) {
207 [ + - ]: 100 : k_arr.push_back(minCoord[2] + crsk*kcrs_size + fnk*kfn_size);
208 : : }
209 : : }
210 : 6 : }
211 : :
212 : : //---------------------------------------------------------------------------//
213 : : // create the mesh vertex coordinates
214 : 6 : void SCDMesh::create_vertex_coords()
215 : : {
216 : 6 : num_i = i_arr.size();
217 : 6 : num_j = j_arr.size();
218 : 6 : num_k = k_arr.size();
219 : 6 : num_verts = num_i*num_j*num_k;
220 : :
221 : 6 : full_coords.resize(3*num_verts);
222 : : unsigned int idx;
223 : :
224 : : // generate the vertex coordinate array in an interleaved fashion
225 [ + + ]: 169 : for (int kval = 0; kval != num_k; kval++) {
226 [ + + ]: 6060 : for (int jval = 0; jval != num_j; jval++) {
227 [ + + ]: 267870 : for (int ival = 0; ival != num_i; ival++) {
228 : 261973 : idx = ival + num_i*jval + num_i*num_j*kval;
229 : 261973 : full_coords[3*idx] = i_arr[ival];
230 : 261973 : full_coords[3*idx + 1] = j_arr[jval];
231 : 261973 : full_coords[3*idx + 2] = k_arr[kval];
232 : : }
233 : : }
234 : : }
235 : 6 : }
236 : :
237 : : //---------------------------------------------------------------------------//
238 : : // create the full mesh representation
239 : 5 : void SCDMesh::create_full_mesh(moab::Range& me_range)
240 : : {
241 : : // create the vertices from the coordinates array
242 [ + - ]: 5 : moab::Range vtx_range;
243 [ + - ][ + - ]: 5 : rval = mk_core()->moab_instance()->create_vertices(&full_coords[0],
[ + - ]
244 : : num_verts,
245 [ + - ]: 5 : vtx_range);
246 [ - + ][ # # ]: 5 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
247 : :
248 [ + + ][ + - ]: 5 : if (geomType == 0) me_range.merge(vtx_range);
249 : :
250 : : // generate hex entities from the vertices
251 : : // the entity handles here should be contiguous with the i blocks
252 : : // moving the fastest and the k blocks moving the slowest
253 : : moab::EntityHandle local_conn[8];
254 : : unsigned int idx;
255 [ + + ]: 144 : for (int kv = 0; kv != num_k - 1; kv++) {
256 [ + + ]: 5392 : for (int jv = 0; jv != num_j - 1; jv++) {
257 [ + + ]: 244186 : for (int iv = 0; iv != num_i - 1; iv++) {
258 : 238933 : idx = iv + jv*num_i + kv*num_i*num_j;
259 [ + - ]: 238933 : local_conn[0] = vtx_range[ idx];
260 [ + - ]: 238933 : local_conn[1] = vtx_range[ idx + 1];
261 [ + - ]: 238933 : local_conn[2] = vtx_range[ idx + 1 + num_i];
262 [ + - ]: 238933 : local_conn[3] = vtx_range[ idx + num_i];
263 [ + - ]: 238933 : local_conn[4] = vtx_range[ idx + num_i*num_j];
264 [ + - ]: 238933 : local_conn[5] = vtx_range[ idx + 1 + num_i*num_j];
265 [ + - ]: 238933 : local_conn[6] = vtx_range[ idx + 1 + num_i + num_i*num_j];
266 [ + - ]: 238933 : local_conn[7] = vtx_range[ idx + num_i + num_i*num_j];
267 : :
268 : : moab::EntityHandle this_elem;
269 [ + - ][ + - ]: 238933 : rval = mk_core()->moab_instance()->create_element( moab::MBHEX,
270 : : local_conn,
271 : : 8,
272 [ + - ]: 238933 : this_elem);
273 [ - + ][ # # ]: 238933 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
274 [ + + ][ + - ]: 238933 : if (geomType == 0) me_range.insert(this_elem);
275 : : }
276 : : }
277 : 5 : }
278 : 5 : }
279 : :
280 : : //---------------------------------------------------------------------------//
281 : : // create mesh using light-weight MOAB ScdInterface
282 : 1 : void SCDMesh::create_light_mesh(moab::Range& me_range)
283 : : {
284 : : // create an instance of the ScdInterface
285 [ + - ][ + - ]: 1 : rval = mk_core()->moab_instance()->query_interface(scdIface);
[ + - ]
286 [ - + ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
287 : :
288 : : // create an instance of the ScdBox class
289 : : // for now the coordinate indexing starts at 0
290 : : // this should be changed to account for multiple geometric instances being meshed
291 : : moab::ScdBox *scd_box;
292 : : rval = scdIface->construct_box(moab::HomCoord(0, 0, 0, 1),
293 : : moab::HomCoord(num_i-1, num_j-1, num_k-1, 1),
294 [ + - ]: 1 : &full_coords[0],
295 : : num_verts,
296 [ + - ][ + - ]: 1 : scd_box);
[ + - ]
297 [ - + ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
298 : :
299 : : // get rid of ScdInterface once we are done
300 [ + - ][ + - ]: 1 : rval = mk_core()->moab_instance()->release_interface(scdIface);
[ + - ]
301 [ - + ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
302 : :
303 : : // add created vertex and hexes
304 [ - + ]: 1 : if (geomType == 1) {
305 : : moab::EntityHandle start;
306 [ # # ]: 0 : start = scd_box->start_vertex();
307 [ # # ][ # # ]: 0 : moab::Range vert_range(start, start + scd_box->num_vertices()-1);
308 [ # # ]: 0 : me_range.merge(vert_range);
309 [ # # ]: 0 : start = scd_box->start_element();
310 [ # # ][ # # ]: 0 : moab::Range hex_range(start, start + scd_box->num_elements()-1);
311 [ # # ]: 0 : me_range.merge(hex_range);
312 : : }
313 : 1 : }
314 : :
315 : : //---------------------------------------------------------------------------//
316 : :
317 [ + - ][ + - ]: 156 : } // end namespace MeshKit
318 : :
319 : : //---------------------------------------------------------------------------//
320 : : // end algs/SCDMesh.cpp
321 : : //---------------------------------------------------------------------------//
|
