Branch data Line data Source code
1 : : #include "meshkit/EBMesher.hpp"
2 : : #include "meshkit/SCDMesh.hpp"
3 : : #include "meshkit/MKCore.hpp"
4 : : #include "meshkit/ModelEnt.hpp"
5 : : #include "meshkit/SizingFunction.hpp"
6 : : #include "moab/ReadUtilIface.hpp"
7 : : #include "meshkit/iGeom.hpp"
8 : : #include "meshkit/RegisterMeshOp.hpp"
9 : :
10 : : #include "moab/Core.hpp"
11 : : #include "moab/Range.hpp"
12 : : #include "moab/CartVect.hpp"
13 : : #include "moab/EntityType.hpp"
14 : : #include "moab/HomXform.hpp"
15 : : #include "moab/GeomTopoTool.hpp"
16 : :
17 : : #include <time.h>
18 : :
19 : : #ifdef HAVE_CGM
20 : : #include "CubitDefines.h"
21 : : #include "GeometryQueryTool.hpp"
22 : : #include "RefFace.hpp"
23 : : #endif
24 : :
25 : : #define PI 3.14159265
26 : : #define ERRORRF(a) {if (iBase_SUCCESS != err) {std::cerr << a << std::endl; return false;}}
27 : :
28 : : double rayDir[3][3] = {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}};
29 : :
30 : : const bool debug_ebmesh = false;
31 : 729 : inline bool less_intersect(const IntersectDist d1, const IntersectDist d2) {
32 : 729 : return d1.distance < d2.distance;
33 : : }
34 : :
35 : : inline bool equal_intersect(const IntersectDist d1, const IntersectDist d2) {
36 : : return std::abs(d1.distance - d2.distance) < 10e-7;
37 : : }
38 : :
39 : : namespace MeshKit
40 : : {
41 : : // static registration of this mesh scheme
42 : : moab::EntityType EBMesher_tps[] = {moab::MBVERTEX, moab::MBHEX, moab::MBMAXTYPE};
43 : 40 : const moab::EntityType* EBMesher::output_types()
44 : 40 : { return EBMesher_tps; }
45 : :
46 : 1 : EBMesher::EBMesher(MKCore *mkcore, const MEntVector &me_vec,
47 : : double size, bool use_geom, int add_layer)
48 [ + - ][ + - ]: 4 : : MeshScheme(mkcore, me_vec), m_nAddLayer(add_layer), m_dInputSize(size), m_bUseGeom(use_geom)
[ + - ][ + - ]
[ + - ][ + + ]
[ + - ][ + - ]
[ + - ][ + -
# # # # #
# ]
49 : : {
50 [ + - ][ + - ]: 1 : m_mesh = mkcore->imesh_instance()->instance();
51 [ + - ][ + - ]: 1 : m_hRootSet = mkcore->imesh_instance()->getRootSet();
52 : 1 : m_nStatus = OUTSIDE;
53 : 1 : m_iStartHex = 0;
54 : 1 : m_hObbTree = NULL;
55 : 1 : m_hTreeRoot = 0;
56 : :
57 : : // create tags with MOAB for dense tags
58 : 1 : int outside = 1;
59 : 1 : const void *out = &outside;
60 : : m_elemStatusTag = get_tag("ELEM_STATUS_TAG", 1,
61 [ + - ]: 1 : moab::MB_TAG_DENSE, moab::MB_TYPE_INTEGER, out);
62 : :
63 : 1 : int length = 1;
64 : 1 : const void *leng = &length;
65 : : m_edgeCutFracLengthTag = get_tag("EDGE_CUT_FRAC_LENGTH_TAG", // # of fractions
66 : : 3,
67 : : //moab::MB_TAG_SPARSE, // using dense for hdf5 exporting performance
68 : : moab::MB_TAG_DENSE,
69 [ + - ]: 1 : moab::MB_TYPE_INTEGER, leng);
70 : :
71 : : m_edgeCutFracTag = get_various_length_tag("EDGE_CUT_FRAC_TAG",
72 : : //moab::MB_TAG_SPARSE,
73 : : moab::MB_TAG_DENSE,
74 [ + - ]: 1 : moab::MB_TYPE_DOUBLE);
75 : :
76 : 1 : int m_id = 1;
77 : 1 : const void *id = &m_id;
78 : : m_volFracLengthTag = get_tag("VOL_FRAC_LENGTH_TAG", 1,
79 [ + - ]: 1 : moab::MB_TAG_SPARSE, moab::MB_TYPE_INTEGER, id);
80 : :
81 : : m_volFracHandleTag = get_various_length_tag("VOL_FRAC_HANDLE_TAG",
82 [ + - ]: 1 : moab::MB_TAG_SPARSE, moab::MB_TYPE_INTEGER);
83 : :
84 : : m_volFracTag = get_various_length_tag("VOL_FRAC_TAG",
85 [ + - ]: 1 : moab::MB_TAG_SPARSE, moab::MB_TYPE_DOUBLE);
86 : :
87 : : // set bounding box size tag
88 : 1 : double bb_default[6] = { 0., 0., 0., 0., 0., 0. };
89 : : m_bbTag = get_tag("BOUNDING_BOX_SIZE", 6,
90 [ + - ]: 1 : moab::MB_TAG_SPARSE, moab::MB_TYPE_DOUBLE, bb_default);
91 : :
92 [ + - ][ + - ]: 1 : m_GeomTopoTool = new moab::GeomTopoTool( moab_instance() );
[ + - ]
93 [ # # ]: 1 : }
94 : :
95 [ + - ][ + + ]: 6 : EBMesher::~EBMesher()
96 : : {
97 [ + - ]: 1 : delete m_GeomTopoTool;
98 [ - + ]: 2 : }
99 : :
100 : 0 : bool EBMesher::can_mesh(ModelEnt *me)
101 : : {
102 [ # # ]: 0 : if (me->dimension() == 3) return true;
103 : 0 : else return false;
104 : : }
105 : :
106 : 0 : bool EBMesher::add_modelent(ModelEnt *model_ent)
107 : : {
108 : : // make sure this represents geometric volumes
109 [ # # ]: 0 : if (3 != model_ent->dimension())
110 [ # # ]: 0 : throw Error(MK_WRONG_DIMENSION, "Wrong dimension entity added to EBMesher.");
111 : :
112 : 0 : return MeshOp::add_modelent(model_ent);
113 : : }
114 : :
115 : 1 : MeshOp *EBMesher::get_scd_mesher()
116 : : {
117 : 1 : MeshOpProxy* proxy = NULL;
118 : 1 : proxy = MKCore::meshop_proxy("SCDMesh");
119 [ - + ][ # # ]: 1 : if (proxy == NULL) throw Error(MK_FAILURE, "Couldn't find a MeshOp capable of producing the given mesh type.");
120 [ + - ][ + - ]: 1 : return mk_core()->construct_meshop(proxy);
121 : : }
122 : :
123 : 1 : void EBMesher::setup_this()
124 : : {
125 : 1 : MeshOp *scd_mesher = NULL;
126 [ + - ]: 1 : std::vector<MeshOp*> meshops;
127 : 1 : bool inserted = false;
128 : :
129 [ # # ][ + - ]: 1 : for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
[ + - ]
130 [ + - ]: 1 : ModelEnt *me = (*sit).first;
131 : 1 : bool found = false;
132 [ + - ][ + - ]: 1 : if (!scd_mesher) scd_mesher = get_scd_mesher();
133 : 1 : meshops.clear();
134 [ + - ]: 1 : me->get_meshops(meshops);
135 : 1 : int n_meshops = meshops.size();
136 : :
137 [ - + ]: 1 : if (n_meshops > 0) {
138 [ # # ]: 0 : for (int i = 0; i < n_meshops; i++) {
139 [ # # ][ # # ]: 0 : if (meshops[i] == scd_mesher) {
140 : 0 : found = true;
141 : 0 : break;
142 : : }
143 : : }
144 : : }
145 : :
146 [ + - ]: 1 : if (!found) { // if no scd meshop
147 : : // add this entity to it, and if first, make sure it's added upstream
148 [ + - ]: 1 : scd_mesher->add_modelent(me);
149 [ + - ]: 1 : if (!inserted) {
150 [ + - ][ + - ]: 1 : mk_core()->insert_node(scd_mesher, this);
151 : 1 : inserted = true;
152 : : }
153 : : }
154 : :
155 : : // no need to traverse all geometry for using whole geometry
156 [ + - ]: 1 : if (m_bUseWholeGeom) break;
157 : : }
158 : :
159 : 1 : SCDMesh *sm = (SCDMesh*) scd_mesher;
160 [ + - ]: 1 : sm->set_interface_scheme(SCDMesh::full);
161 [ + - ]: 1 : sm->set_grid_scheme(SCDMesh::cfMesh);
162 : :
163 [ + - ]: 1 : sm->set_axis_scheme(SCDMesh::cartesian);
164 [ + - ]: 1 : sm->set_box_increase_ratio(m_boxIncrease); // add some extra layer to box
165 : :
166 [ + - ][ + - ]: 1 : if (m_bUseWholeGeom) sm->set_geometry_scheme(SCDMesh::all);
167 [ # # ]: 0 : else sm->set_geometry_scheme(SCDMesh::individual);
168 : :
169 : : // use mesh based geometry in SCDMesh and make tree to get box dimension
170 [ - + ]: 1 : if (m_bUseMeshGeom) {
171 [ # # ]: 0 : sm->use_mesh_geometry(true);
172 [ # # ]: 0 : sm->set_cart_box_min_max(m_minCoord, m_maxCoord, m_boxIncrease);
173 : : }
174 : :
175 : : // set # of intervals for 3 directions
176 [ + - ]: 2 : std::vector<int> fine_i (m_nDiv[0], 1);
177 [ + - ]: 1 : sm->set_coarse_i_grid(m_nDiv[0]);
178 [ + - ][ + - ]: 1 : sm->set_fine_i_grid(fine_i);
179 [ + - ]: 2 : std::vector<int> fine_j (m_nDiv[1], 1);
180 [ + - ]: 1 : sm->set_coarse_j_grid(m_nDiv[1]);
181 [ + - ][ + - ]: 1 : sm->set_fine_j_grid(fine_j);
182 [ + - ]: 2 : std::vector<int> fine_k (m_nDiv[2], 1);
183 [ + - ]: 1 : sm->set_coarse_k_grid(m_nDiv[2]);
184 [ + - ][ + - ]: 2 : sm->set_fine_k_grid(fine_k);
185 : 1 : }
186 : :
187 : 1 : void EBMesher::execute_this()
188 : : {
189 : : iBase_ErrorType err;
190 [ + - ][ + - ]: 1 : GraphNode *scd_node = other_node(in_arcs());
191 : :
192 : 1 : double obb_time = .0;
193 : 1 : double intersection_time = .0;
194 : 1 : double set_info_time = .0;
195 : : time_t time1, time2, time3, time4;
196 : : unsigned long long mem1, mem3, mem4;
197 [ + - ][ + - ]: 1 : moab_instance()->estimated_memory_use(0, 0, 0, &mem1);
198 : : moab::ErrorCode rval;
199 : :
200 : : if (debug_ebmesh) {
201 : : rval = mk_core()->moab_instance()->write_mesh("input.vtk");
202 : : MBERRCHK(rval, mk_core()->moab_instance());
203 : : }
204 : :
205 : 1 : time(&time1);
206 [ + - ][ + - ]: 1 : if (m_bUseGeom && !m_bUseMeshGeom) { // construct obb tree for geometry surfaces and volumes by GeomTopoTool
207 [ + - ]: 1 : rval = m_GeomTopoTool->construct_obb_trees(m_bUseWholeGeom);
208 [ - + ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
209 : : }
210 : 1 : time(&time2);
211 : 1 : obb_time += difftime(time2, time1);
212 : :
213 [ # # ][ + - ]: 1 : for (MEntSelection::iterator sit = mentSelection.begin(); sit != mentSelection.end(); sit++) {
[ + - ]
214 : : // 1. set division
215 : : double box_min_max[6];
216 [ + - ]: 1 : if (m_bUseWholeGeom) {
217 [ + - ]: 1 : static_cast<SCDMesh*> (scd_node)->get_box_dimension(box_min_max, &box_min_max[3]);
218 : : }
219 : : else {
220 [ # # ][ # # ]: 0 : err = mk_core()->imesh_instance()->getData(reinterpret_cast< iBase_EntityHandle >
221 [ # # ][ # # ]: 0 : (sit ->first->mesh_handle()), m_bbTag, &box_min_max[0]);
[ # # ]
222 [ # # ]: 0 : IBERRCHK(err, "Failed to get hexes.\n");
223 : : }
224 [ + - ]: 1 : set_division(box_min_max, &box_min_max[3]);
225 : :
226 : : // 2. set or construct obb tree
227 : 1 : time(&time1);
228 [ + - ][ + - ]: 1 : set_tree_root(sit->first);
229 : 1 : time(&time2);
230 : 1 : obb_time += difftime(time2, time1);
231 : :
232 : : // 3. find intersected geometry surfaces by rays
233 [ + - ]: 1 : find_intersections();
234 : 1 : time(&time3);
235 : 1 : intersection_time += difftime(time3, time2);
236 [ + - ][ + - ]: 1 : moab_instance()->estimated_memory_use(0, 0, 0, &mem3);
237 : :
238 : : // 4. set hex status and boundary hex cut fraction info
239 [ + - ]: 1 : set_tag_info();
240 : 1 : time(&time4);
241 : 1 : set_info_time += difftime(time4, time3);
242 [ + - ][ + - ]: 1 : moab_instance()->estimated_memory_use(0, 0, 0, &mem4);
243 : :
244 [ + - ]: 1 : if (m_bUseWholeGeom) break;
245 : : }
246 : : if (debug_ebmesh) {
247 : : std::cout << "OBB_tree_construct_time: " << obb_time
248 : : << ", intersection_time: " << intersection_time
249 : : << ", set_info_time: " << set_info_time << std::endl;
250 : : std::cout << "start_memory: " << mem1
251 : : //<< ", OBB_tree_construct_moemory: " << mem2
252 : : << ", intersection_memory: " << mem3
253 : : << ", set_info_memory: " << mem4
254 : : << std::endl;
255 : : }
256 : 1 : }
257 : :
258 : 1 : void EBMesher::set_num_interval(int* n_interval)
259 : : {
260 [ + + ]: 4 : for (int i = 0; i < 3; i++) {
261 : 3 : m_nDiv[i] = n_interval[i];
262 : : }
263 : 1 : }
264 : :
265 : 1 : void EBMesher::set_tree_root(ModelEnt* me)
266 : : {
267 : : moab::ErrorCode rval;
268 [ + - ]: 1 : if (m_bUseGeom) {
269 [ + - ]: 1 : if (m_hObbTree == NULL) {
270 : 1 : m_hObbTree = m_GeomTopoTool->obb_tree();
271 : : }
272 [ + - ]: 1 : if (m_bUseWholeGeom) {
273 [ + - ]: 1 : if (!m_bUseMeshGeom) m_hTreeRoot = m_GeomTopoTool->get_one_vol_root();
274 : : }
275 : : else {
276 : 0 : rval = m_GeomTopoTool->get_root(me->mesh_handle(), m_hTreeRoot);
277 [ # # ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
278 : : }
279 : : }
280 : : else { // facet data input case
281 : : moab::EntityHandle meshset;
282 [ # # ]: 0 : if (m_bUseWholeGeom) meshset = 0;
283 [ # # ]: 0 : else meshset = me->mesh_handle();
284 : :
285 : : // get triangles
286 [ # # ]: 0 : moab::Range tris;
287 [ # # ][ # # ]: 0 : rval = moab_instance()->get_entities_by_dimension(meshset, 2, tris);
288 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
289 : :
290 : : // make tree
291 [ # # ]: 0 : if (m_hObbTree == NULL) {
292 [ # # ][ # # ]: 0 : m_hObbTree = new moab::OrientedBoxTreeTool( moab_instance() );
[ # # ]
293 : : }
294 [ # # ]: 0 : rval = m_hObbTree->build(tris, m_hTreeRoot);
295 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
296 : : }
297 : 1 : }
298 : :
299 : 1 : void EBMesher::find_intersections()
300 : : {
301 : 1 : std::cout << "starting to find intersections." << std::endl;
302 : : int i;
303 : : //m_vnHexStatus.resize(m_nHex, 1); // initialize all hex as outside
304 [ + - ]: 1 : m_vnHexStatus.resize(m_nHex, 0); // initialize all hex as inside
305 : :
306 : : // fire rays to 3 directions
307 [ + + ]: 4 : for (i = 0; i < 3; i++) {
308 : : //m_vnEdgeStatus[i].resize(m_nDiv[i]*m_nDiv[(i + 1)%3]*m_nDiv[(i + 2)%3], OUTSIDE);
309 [ + - ]: 3 : m_vnEdgeStatus[i].resize(m_nDiv[i]*m_nDiv[(i + 1)%3]*m_nDiv[(i + 2)%3], INSIDE);
310 : 3 : fire_rays(i);
311 : : }
312 : 1 : std::cout << "finished to find intersections." << std::endl;
313 : 1 : }
314 : :
315 : 1 : void EBMesher::export_mesh(const char* file_name, bool separate)
316 : : {
317 : : // get all hexes
318 : : time_t time1, time2, time3;
319 : 1 : time(&time1);
320 : : int i;
321 : : iBase_ErrorType err;
322 [ + - ]: 1 : std::vector<iBase_EntityHandle> hex_handles;
323 : : err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
324 [ + - ][ + - ]: 1 : iMesh_HEXAHEDRON, hex_handles);
[ + - ]
325 [ + - ]: 1 : IBERRCHK(err, "Failed to get hexes.\n");
326 : : /*
327 : : int hex_size = hex_handles.size();
328 : : int* hex_status = new int[hex_size];
329 : : int* hex_status = NULL;
330 : : std::vector<int> hex_status(hex_size);
331 : : err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
332 : : m_elemStatusTag, &hex_status[0]);*/
333 : :
334 : : int error;
335 : 1 : int hex_size = hex_handles.size();
336 [ + - ][ + - ]: 1 : int* hex_status = new int[hex_size];
337 : 1 : int hex_status_alloc = sizeof(int)*hex_size;
338 : 1 : int hex_status_size = 0;
339 [ + - ]: 1 : iMesh_getIntArrData(m_mesh, &hex_handles[0],
340 : : hex_size, m_elemStatusTag,
341 : : &hex_status, &hex_status_alloc,
342 [ + - ]: 1 : &hex_status_size, &error);
343 [ + - ]: 1 : IBERRCHK(error, "Failed to get hex status.\n");
344 : 1 : time(&time2);
345 : :
346 [ - + ]: 1 : if (separate) {
347 : 0 : int n_inside_hex = 0;
348 : 0 : int n_outside_hex = 0;
349 : 0 : int n_boundary_hex = 0;
350 : : int hex_stat;
351 : : (void) hex_stat;
352 [ # # ][ # # ]: 0 : std::vector<iBase_EntityHandle> insideHex, outsideHex, bndrHex;
[ # # ]
353 [ # # ]: 0 : for (i = 0; i < hex_size; i++) {
354 [ # # ]: 0 : if (hex_status[i] == 0) {
355 : 0 : n_inside_hex++;
356 : 0 : hex_stat = 0;
357 [ # # ][ # # ]: 0 : insideHex.push_back(hex_handles[i]);
358 : : }
359 [ # # ]: 0 : else if (hex_status[i] == 1) {
360 : 0 : n_outside_hex++;
361 : 0 : hex_stat = 1;
362 [ # # ][ # # ]: 0 : outsideHex.push_back(hex_handles[i]);
363 : : }
364 [ # # ]: 0 : else if (hex_status[i] == 2) {
365 : 0 : n_boundary_hex++;
366 : 0 : hex_stat = 2;
367 [ # # ][ # # ]: 0 : bndrHex.push_back(hex_handles[i]);
368 : : }
369 : : else {
370 [ # # ]: 0 : throw Error(MK_FAILURE, "Hex element status should be inside/outside/boundary.");
371 : : }
372 : : }
373 : :
374 [ # # ][ # # ]: 0 : std::cout << "# of exported inside hex:" << n_inside_hex
375 [ # # ][ # # ]: 0 : << ", # of exported outside hex:" << n_outside_hex
376 [ # # ][ # # ]: 0 : << ", # of exported boundary hex:" << n_boundary_hex
377 [ # # ]: 0 : << ", geom vol:"
378 [ # # ]: 0 : << n_inside_hex*m_dIntervalSize[0]*m_dIntervalSize[1]*m_dIntervalSize[2]
379 [ # # ][ # # ]: 0 : << ", vox vol:" << hex_size*m_dIntervalSize[0]*m_dIntervalSize[1]*m_dIntervalSize[2]
380 [ # # ]: 0 : << std::endl;
381 : 0 : time(&time3);
382 : :
383 : : // save inside/outside/boundary elements separately
384 [ # # ]: 0 : if (n_inside_hex > 0) {
385 [ # # ][ # # ]: 0 : write_mesh(file_name, 0, &insideHex[0], n_inside_hex);
386 : : }
387 [ # # ]: 0 : if (n_boundary_hex > 0) {
388 [ # # ][ # # ]: 0 : write_mesh(file_name, 2, &bndrHex[0], n_boundary_hex);
389 : : }
390 [ # # ]: 0 : if (n_outside_hex > 0) {
391 [ # # ][ # # ]: 0 : write_mesh(file_name, 1, &outsideHex[0], n_outside_hex);
392 : : }
393 : :
394 : : if (debug_ebmesh) {
395 : : std::cout << "hex_handle_get_time: "
396 : : << difftime(time2, time1)
397 : : << ", separate_write_time: "
398 : : << difftime(time3, time2)
399 : : << std::endl;
400 : 0 : }
401 : : }
402 : : else {
403 [ + - ][ + - ]: 1 : write_mesh(file_name, 3, &hex_handles[0], hex_size);
404 : 1 : time3 = clock();
405 : : if (debug_ebmesh) {
406 : : std::cout << "hex_handle_get_time: "
407 : : << difftime(time2, time1)
408 : : << ", actual_write_time: "
409 : : << difftime(time3, time2)
410 : : << std::endl;
411 : : }
412 : : }
413 [ + - ]: 1 : delete [] hex_status;
414 : 1 : }
415 : :
416 : 1 : void EBMesher::write_mesh(const char* file_name, int type,
417 : : iBase_EntityHandle* handles, int& n_elem)
418 : : {
419 : : time_t time1, time2, time3;
420 : 1 : time(&time1);
421 : 1 : int is_list = 1;
422 : : moab::ErrorCode rval;
423 : : iBase_EntitySetHandle set;
424 : : iBase_ErrorType err;
425 : :
426 [ + - ][ + - ]: 1 : err = mk_core()->imesh_instance()->createEntSet(is_list, set);
[ + - ]
427 [ + - ]: 1 : IBERRCHK(err, "Couldn't create set.");
428 : :
429 [ + - ][ + - ]: 1 : err = mk_core()->imesh_instance()->addEntArrToSet(handles, n_elem, set);
[ + - ]
430 [ + - ]: 1 : IBERRCHK(err, "Couldn't add hexes to set.");
431 : 1 : time(&time2);
432 : :
433 [ + - ]: 1 : std::string out_name;
434 [ + - ][ + - ]: 2 : std::stringstream ss;
435 [ - + ][ # # ]: 1 : if (type == 0) ss << "inside_";
436 [ - + ][ # # ]: 1 : else if (type == 1) ss << "outside_";
437 [ - + ][ # # ]: 1 : else if (type == 2) ss << "boundary_";
438 [ + - ]: 1 : ss << file_name;
439 [ + - ]: 1 : ss >> out_name;
440 : :
441 [ + - ]: 1 : rval = moab_instance()->write_mesh(out_name.c_str(),
442 [ + - ]: 1 : (const moab::EntityHandle*) &set, 1);
443 [ - + ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
444 : :
445 [ + - ][ + - ]: 1 : std::cout << "Elements are exported." << std::endl;
446 : 1 : time(&time3);
447 : :
448 : : if (debug_ebmesh) {
449 : : std::cout << "set_creation_time: "
450 : : << difftime(time2, time1)
451 : : << ", write_time: "
452 : : << difftime(time3, time2)
453 : : << std::endl;
454 : 1 : }
455 : 1 : }
456 : :
457 : 573 : EdgeStatus EBMesher::get_edge_status(const double dP, int& iSkip)
458 : : {
459 [ + + ]: 573 : if (m_nStatus == INSIDE) { // previous inside
460 [ - + ]: 207 : if (dP < m_dSecondP) {
461 : 0 : m_nMove = 0;
462 : 0 : iSkip = m_iSecondP;
463 : 0 : return INSIDE;
464 : : }
465 : : else {
466 [ - + ]: 207 : if (is_shared_overlapped_surf(m_iInter - 1)) m_nMove = 1;
467 : 207 : else m_nMove = 2;
468 : :
469 : : // skip shared and overlapped interesections
470 [ - + ]: 414 : while (m_vIntersection[m_iInter].distance -
471 [ + + - + ]: 300 : m_vIntersection[m_iInter - 1].distance < 1e-7 &&
472 : 93 : (unsigned int) (m_iInter + 1) < m_vIntersection.size()) m_iInter++;
473 : :
474 : 207 : return BOUNDARY;
475 : : }
476 : : }
477 [ + + ]: 366 : else if (m_nStatus == OUTSIDE) { // previous outside
478 [ - + ]: 159 : if (dP < m_dFirstP) {
479 : 0 : m_nMove = 0;
480 : 0 : return OUTSIDE;
481 : : }
482 : : else {
483 [ + - ]: 159 : if (dP < m_dSecondP) m_nMove = 0;
484 [ # # ]: 0 : else if (is_shared_overlapped_surf(m_iInter - 1)) m_nMove = 1;
485 : 0 : else m_nMove = 2;
486 : :
487 : : // skip shared and overlapped interesections
488 [ - + ]: 318 : while (m_vIntersection[m_iInter].distance -
489 [ + + - + ]: 218 : m_vIntersection[m_iInter - 1].distance < 1e-7 &&
490 : 59 : (unsigned int) (m_iInter + 1) < m_vIntersection.size()) m_iInter++;
491 : :
492 : 159 : return BOUNDARY;
493 : : }
494 : : }
495 [ + - ]: 207 : else if (m_nStatus == BOUNDARY) { // previous boundary
496 [ - + ]: 207 : if (dP < m_dFirstP) {
497 : 0 : m_nMove = 0;
498 : 0 : iSkip = m_iFirstP;
499 : 0 : return OUTSIDE;
500 : : }
501 : : else {
502 [ + - ]: 207 : if (dP < m_dSecondP) {
503 : 207 : m_nMove = 0;
504 [ - + ]: 207 : if (m_prevPnt < m_dFirstP) return BOUNDARY;
505 : : else {
506 : 207 : iSkip = m_iSecondP;
507 : 207 : return INSIDE;
508 : : }
509 : : }
510 : : else {
511 [ # # ]: 0 : if (is_shared_overlapped_surf(m_iInter - 1)) m_nMove = 1;
512 : 0 : else m_nMove = 2;
513 : :
514 : : // skip shared and overlapped interesections
515 [ # # ]: 0 : while (m_vIntersection[m_iInter].distance -
516 [ # # # # ]: 0 : m_vIntersection[m_iInter - 1].distance < 1e-7 &&
517 : 0 : (unsigned int) (m_iInter + 1) < m_vIntersection.size()) m_iInter++;
518 : :
519 [ # # ]: 0 : if (m_prevPnt < m_dSecondP) return BOUNDARY;
520 : 0 : else return OUTSIDE;
521 : : }
522 : : }
523 : : }
524 : : else {
525 : 0 : std::cerr << "Couldn't get edge status." << std::endl;
526 : 0 : return INSIDE;
527 : : }
528 : : }
529 : :
530 : 1413 : bool EBMesher::set_neighbor_hex_status(int dir, int i, int j, int k)
531 : : {
532 : : int iElem;
533 : 1413 : int otherDir1 = (dir + 1)%3;
534 : 1413 : int otherDir2 = (dir + 2)%3;
535 : :
536 [ + + ]: 1413 : if (dir == 0) { // x coordinate ray
537 : 471 : m_iElem = j*m_nDiv[0]*m_nDiv[1] + i*m_nDiv[0] + k;
538 [ - + ]: 471 : if (!set_hex_status(m_iElem, m_nStatus, dir)) return false;
539 : 471 : iElem = m_iElem - m_nDiv[dir];
540 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
541 : 471 : iElem -= m_nDiv[dir]*m_nDiv[otherDir1];
542 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
543 : 471 : iElem += m_nDiv[dir];
544 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
545 : : }
546 [ + + ]: 942 : else if (dir == 1) { // y coordinate ray
547 : 471 : m_iElem = i*m_nDiv[1]*m_nDiv[0] + k*m_nDiv[0] + j;
548 [ - + ]: 471 : if (!set_hex_status(m_iElem, m_nStatus, dir)) return false;
549 : 471 : iElem = m_iElem - 1;
550 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
551 : 471 : iElem -= m_nDiv[dir]*m_nDiv[otherDir2];
552 [ - + ]: 471 : if (!set_hex_status(iElem++, m_nStatus, dir)) return false;
553 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
554 : : }
555 [ + - ]: 471 : else if (dir == 2) { // z coordinate ray
556 : 471 : m_iElem = k*m_nDiv[0]*m_nDiv[1] + j*m_nDiv[0] + i;
557 [ - + ]: 471 : if (!set_hex_status(m_iElem, m_nStatus, dir)) return false;
558 : 471 : iElem = m_iElem - 1;
559 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
560 : 471 : iElem -= m_nDiv[otherDir1];
561 [ - + ]: 471 : if (!set_hex_status(iElem++, m_nStatus, dir)) return false;
562 [ - + ]: 471 : if (!set_hex_status(iElem, m_nStatus, dir)) return false;
563 : : }
564 : :
565 : 1413 : return true;
566 : : }
567 : :
568 : 5652 : bool EBMesher::set_hex_status(int index, EdgeStatus value, int dir)
569 : : {
570 [ + - ][ - + ]: 5652 : if (index < 0 || index > m_nHex - 1) {
571 : 0 : return false;
572 : : }
573 : :
574 [ + + ]: 5652 : if (m_vnHexStatus[index] != 2) {
575 [ + + ]: 3249 : if (value == INSIDE) {
576 : 457 : m_vnHexStatus[index] = 0;
577 : : }
578 [ + + ]: 2792 : else if (value == OUTSIDE) {
579 : 2400 : m_vnHexStatus[index] = 1;
580 : : }
581 [ + - ]: 392 : else if (value == BOUNDARY) {
582 : 392 : m_vnHexStatus[index] = 2;
583 : : }
584 : : }
585 : :
586 : 5652 : return true;
587 : : }
588 : :
589 : 366 : bool EBMesher::set_edge_status(int dir)
590 : : {
591 : : // set boundary cut information to edge
592 [ + - ]: 366 : std::vector<double> vdCutFraction;
593 [ + + ][ + - ]: 366 : if (m_nMove == 0) vdCutFraction.push_back((m_curPnt - m_dFirstP)/m_dIntervalSize[dir] - 1.);
594 [ + - ][ + - ]: 207 : else if (m_nMove == 1) vdCutFraction.push_back(1. - (m_curPnt - m_dSecondP)/m_dIntervalSize[dir]);
595 [ # # ]: 0 : else if (m_nMove == 2) {
596 [ # # ]: 0 : vdCutFraction.push_back(1. - (m_curPnt - m_dSecondP)/m_dIntervalSize[dir]);
597 [ # # ]: 0 : if (m_dFirstP < m_curPnt) {
598 [ # # ]: 0 : vdCutFraction.push_back((m_curPnt - m_dFirstP)/m_dIntervalSize[dir] - 1.);
599 : : }
600 : : }
601 : :
602 [ + - ]: 366 : std::map<int, CutFraction>::iterator iter = m_mdCutFraction.find(m_iElem);
603 [ + - ][ + + ]: 366 : if (iter == m_mdCutFraction.end()) { // not exist
604 [ + - ]: 268 : CutFraction cFraction(dir, vdCutFraction);
605 [ + - ][ + - ]: 268 : m_mdCutFraction[m_iElem] = cFraction;
606 : : }
607 : : else { // exist
608 [ + - ][ + - ]: 98 : iter->second.add(dir, vdCutFraction);
609 : : }
610 : :
611 : : //m_nFraction += vdCutFraction.size();
612 : :
613 : 366 : return true;
614 : : }
615 : :
616 : 0 : void EBMesher::set_division()
617 : : {
618 : : int i, j;
619 [ # # ][ # # ]: 0 : moab::CartVect box_center, box_axis1, box_axis2, box_axis3,
[ # # ][ # # ]
620 [ # # ]: 0 : min_cart_box(HUGE_VAL, HUGE_VAL, HUGE_VAL),
621 [ # # ]: 0 : max_cart_box(0., 0., 0.);
622 : :
623 : : moab::ErrorCode rval = m_hObbTree->box(m_hTreeRoot, box_center.array(),
624 : : box_axis1.array(), box_axis2.array(),
625 [ # # ][ # # ]: 0 : box_axis3.array());
[ # # ][ # # ]
[ # # ]
626 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
627 : :
628 : : // cartesian box corners
629 [ # # ][ # # ]: 0 : moab::CartVect corners[8] = {box_center + box_axis1 + box_axis2 + box_axis3,
630 [ # # ][ # # ]: 0 : box_center + box_axis1 + box_axis2 - box_axis3,
631 [ # # ][ # # ]: 0 : box_center + box_axis1 - box_axis2 + box_axis3,
632 [ # # ][ # # ]: 0 : box_center + box_axis1 - box_axis2 - box_axis3,
633 [ # # ][ # # ]: 0 : box_center - box_axis1 + box_axis2 + box_axis3,
634 [ # # ][ # # ]: 0 : box_center - box_axis1 + box_axis2 - box_axis3,
635 [ # # ][ # # ]: 0 : box_center - box_axis1 - box_axis2 + box_axis3,
636 [ # # ][ # # ]: 0 : box_center - box_axis1 - box_axis2 - box_axis3};
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
637 : :
638 : : // get the max, min cartesian box corners
639 [ # # ]: 0 : for (i = 0; i < 8; i++) {
640 [ # # ]: 0 : for (j = 0; j < 3; j++) {
641 [ # # ][ # # ]: 0 : if (corners[i][j] < min_cart_box[j]) min_cart_box[j] = corners[i][j];
[ # # ][ # # ]
[ # # ]
642 [ # # ][ # # ]: 0 : if (corners[i][j] > max_cart_box[j]) max_cart_box[j] = corners[i][j];
[ # # ][ # # ]
[ # # ]
643 : : }
644 : : }
645 [ # # ]: 0 : moab::CartVect length = max_cart_box - min_cart_box;
646 : :
647 : : // default value is adjusted to large geometry file
648 : : // interval_size_estimate : 2*L*sqrt(2*PI*sqrt(2)/# of tris)
649 [ # # ]: 0 : if (m_dInputSize < 0.) {
650 : : int n_tri;
651 [ # # ]: 0 : rval = moab_instance()->
652 : : get_number_entities_by_dimension(reinterpret_cast<moab::EntityHandle>(m_hRootSet),
653 [ # # ]: 0 : 2, n_tri);
654 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
655 : :
656 [ # # ][ # # ]: 0 : double box_length_ave = (length[0] + length[1] + length[2])/3.;
[ # # ]
657 : 0 : m_dInputSize = 2.*box_length_ave*sqrt(8.886/n_tri);
658 : : }
659 : :
660 [ # # ]: 0 : for (i = 0; i < 3; i++) {
661 [ # # ]: 0 : m_nDiv[i] = length[i]/m_dInputSize;
662 [ # # ]: 0 : if (m_nDiv[i] < m_nAddLayer) m_nDiv[i] = m_nAddLayer; // make "m_nAddLayer" elements larger than bounding box
663 : 0 : m_dIntervalSize[i] = m_dInputSize;
664 : : //if (m_nDiv[i]*.07 > m_nAddLayer) m_nDiv[i] += m_nDiv[i]*.07;
665 : : //else m_nDiv[i] += m_nAddLayer;
666 : 0 : m_nDiv[i] += m_nAddLayer;
667 : 0 : m_nNode[i] = m_nDiv[i] + 1;
668 [ # # ]: 0 : m_origin_coords[i] = box_center[i] - .5*m_nDiv[i]*m_dIntervalSize[i];
669 : : }
670 : :
671 : 0 : m_nHex = m_nDiv[0]*m_nDiv[1]*m_nDiv[2];
672 : :
673 [ # # ][ # # ]: 0 : std::cout << "# of hex: " << m_nHex << ", interval size: "
[ # # ]
674 [ # # ][ # # ]: 0 : << m_dInputSize << std::endl;
675 : :
676 [ # # ][ # # ]: 0 : std::cout << "# of division: " << m_nDiv[0] << ","
[ # # ]
677 [ # # ][ # # ]: 0 : << m_nDiv[1] << "," << m_nDiv[2] << std::endl;
[ # # ][ # # ]
678 : 0 : }
679 : :
680 : 1 : int EBMesher::set_division(double* min, double* max)
681 : : {
682 [ + + ]: 4 : for (int i = 0; i < 3; i++) {
683 : 3 : m_dIntervalSize[i] = (max[i] - min[i])/m_nDiv[i];
684 : 3 : m_nNode[i] = m_nDiv[i] + 1;
685 : 3 : m_origin_coords[i] = min[i];
686 : : }
687 : :
688 : 1 : m_nHex = m_nDiv[0]*m_nDiv[1]*m_nDiv[2];
689 : :
690 : 1 : std::cout << "# of hex: " << m_nHex << ", interval_size: "
691 : 1 : << m_dIntervalSize[0] << ", " << m_dIntervalSize[1] << ", "
692 : 1 : << m_dIntervalSize[2] << std::endl;
693 : :
694 : 1 : std::cout << "# of division: " << m_nDiv[0] << ","
695 : 1 : << m_nDiv[1] << "," << m_nDiv[2] << std::endl;
696 : :
697 : 1 : return iBase_SUCCESS;
698 : : }
699 : :
700 : 0 : int EBMesher::make_scd_hexes()
701 : : {
702 : : // create vertex and hex sequences
703 : : int i;
704 : : double max_coords[3];
705 : : (void) max_coords;
706 [ # # ]: 0 : for (i = 0; i < 3; i++) {
707 : 0 : max_coords[i] = m_origin_coords[i] + m_nDiv[i]*m_dIntervalSize[i];
708 : : }
709 : :
710 [ # # ]: 0 : moab::HomCoord coord_min(0, 0, 0);
711 [ # # ]: 0 : moab::HomCoord coord_max(m_nDiv[0], m_nDiv[1], m_nDiv[2]);
712 : 0 : moab::EntitySequence* vertex_seq = NULL;
713 : 0 : moab::EntitySequence* cell_seq = NULL;
714 : : moab::EntityHandle vs, cs;
715 [ # # ][ # # ]: 0 : moab::Core *mbcore = dynamic_cast<moab::Core*>(moab_instance());
716 : :
717 [ # # ]: 0 : moab::ErrorCode rval = mbcore->create_scd_sequence(coord_min, coord_max, moab::MBVERTEX, 1, vs, vertex_seq);
718 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
719 : :
720 [ # # ]: 0 : mbcore->create_scd_sequence(coord_min, coord_max, moab::MBHEX, 1, cs, cell_seq);
721 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
722 : :
723 [ # # ][ # # ]: 0 : moab::HomCoord p2(coord_max.i(), coord_min.j(), coord_min.k());
[ # # ][ # # ]
724 [ # # ][ # # ]: 0 : moab::HomCoord p3(coord_min.i(), coord_max.j(), coord_min.k());
[ # # ][ # # ]
725 : :
726 : : rval = mbcore->add_vsequence(vertex_seq, cell_seq, coord_min, coord_min,
727 [ # # ]: 0 : p2, p2, p3, p3);
728 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
729 : :
730 : 0 : int nTotNode = m_nNode[0]*m_nNode[1]*m_nNode[2];
731 : : int ii, jj, kk;
732 : : moab::EntityHandle handle;
733 : : double dumv[3];
734 [ # # ]: 0 : for (i = 0, handle = vs; i < nTotNode; i++, handle++) {
735 : 0 : ii = (i%(m_nNode[0]*m_nNode[1]))%m_nNode[0];
736 : 0 : jj = (i%(m_nNode[0]*m_nNode[1]))/m_nNode[0];
737 : 0 : kk = i/m_nNode[0]/m_nNode[1];
738 : 0 : dumv[0] = m_origin_coords[0] + ii*m_dIntervalSize[0];
739 : 0 : dumv[1] = m_origin_coords[1] + jj*m_dIntervalSize[1];
740 : 0 : dumv[2] = m_origin_coords[2] + kk*m_dIntervalSize[2];
741 [ # # ][ # # ]: 0 : rval = moab_instance()->set_coords(&handle, 1, dumv);
742 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
743 : : }
744 : :
745 [ # # ][ # # ]: 0 : m_iStartHex = moab_instance()->id_from_handle(cs);
746 : :
747 : 0 : return iBase_SUCCESS;
748 : : }
749 : :
750 : 4 : iBase_TagHandle EBMesher::get_tag(const char* name, int size,
751 : : unsigned store, moab::DataType type,
752 : : const void* def_value,
753 : : bool create_if_missing)
754 : : {
755 : 4 : moab::Tag retval = 0;
756 : : /*moab::ErrorCode result = moab_instance()->tag_create(name, size, store, type,
757 : : retval, def_value,
758 : : create_if_missing);*/
759 : 4 : store = store | moab::MB_TAG_CREAT;
760 [ - + ]: 4 : if(!create_if_missing)
761 : 0 : store = store | moab::MB_TAG_EXCL;
762 [ + - ]: 4 : moab::ErrorCode result = moab_instance()->tag_get_handle(name, size, type, retval, store,
763 [ + - ]: 4 : def_value);
764 [ + - ][ - + ]: 4 : if (create_if_missing && moab::MB_SUCCESS != result) {
765 [ # # ][ # # ]: 0 : std::cerr << "Couldn't find nor create tag named " << name << std::endl;
[ # # ]
766 : : }
767 : :
768 : 4 : return (iBase_TagHandle) retval;
769 : : }
770 : :
771 : 3 : iBase_TagHandle EBMesher::get_various_length_tag(const char* name,
772 : : unsigned store, moab::DataType type)
773 : : {
774 : 3 : moab::Tag retval = 0;
775 : 3 : store = store | moab::MB_TAG_VARLEN | moab::MB_TAG_CREAT;
776 [ + - ]: 3 : moab::ErrorCode result = moab_instance()->
777 [ + - ]: 3 : tag_get_handle( name, 1, type, retval, store);
778 [ - + ]: 3 : if (moab::MB_SUCCESS != result) {
779 [ # # ][ # # ]: 0 : std::cerr << "Couldn't find nor create tag named " << name << std::endl;
[ # # ]
780 : : }
781 : :
782 : 3 : return (iBase_TagHandle) retval;
783 : : }
784 : :
785 : 1 : void EBMesher::set_tag_info()
786 : : {
787 : : // get all hexes
788 : : int i, j, k;
789 : : iBase_ErrorType err;
790 [ + - ]: 1 : std::vector<iBase_EntityHandle> hex_handles;
791 : : err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
792 [ + - ][ + - ]: 1 : iMesh_HEXAHEDRON, hex_handles);
[ + - ]
793 [ + - ]: 1 : IBERRCHK(err, "Failed to get hexes.\n");
794 : :
795 [ + - ][ + - ]: 1 : err = mk_core()->imesh_instance()->setIntArrData(&hex_handles[0], m_nHex,
[ + - ]
796 : : m_elemStatusTag,
797 [ + - ][ + - ]: 2 : &m_vnHexStatus[0]);
798 [ + - ]: 1 : IBERRCHK(err, "Failed to set hex element status data.");
799 : :
800 : : // set cut fraction info to boundary hexes
801 : 1 : int nBndrHex = m_mdCutFraction.size();
802 [ + - ]: 2 : std::vector<iBase_EntityHandle> hvBndrHex(nBndrHex);
803 [ + - ][ + - ]: 1 : int* frac_size = new int[nBndrHex];
804 [ + - ][ + - ]: 1 : int* frac_leng = new int[3*nBndrHex];
805 : 1 : int nFracSize, nTempSize, iHex, nTolFrac = 0;
806 : 1 : int nDoubleSize = sizeof(double);
807 : 1 : std::map<int, CutFraction>::iterator iter = m_mdCutFraction.begin();
808 : 1 : std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
809 [ + - ][ + - ]: 269 : for (i = 0; iter != end_iter; iter++, i++) {
[ + + ]
810 [ + - ]: 268 : iHex = iter->first;
811 [ + - ][ + - ]: 268 : hvBndrHex[i] = hex_handles[iHex];
812 : :
813 : 268 : nFracSize = 0;
814 [ + + ]: 1072 : for (j = 0; j < 3; j++) {
815 [ + - ]: 804 : nTempSize = iter->second.fraction[j].size();
816 : 804 : nFracSize += nTempSize;
817 : 804 : frac_leng[3*i + j] = nTempSize;
818 : : }
819 : 268 : frac_size[i] = nDoubleSize*nFracSize; // sizeof(double)*
820 : 268 : nTolFrac += nFracSize;
821 : : }
822 : :
823 : 1 : int iFrac = 0;
824 [ + - ][ + - ]: 1 : double* fractions = new double[nTolFrac];
825 [ + - ]: 2 : std::vector<const void*> frac_data_pointer(nBndrHex);
826 : 1 : iter = m_mdCutFraction.begin();
827 [ + - ][ + - ]: 269 : for (i = 0; iter != end_iter; iter++, i++) {
[ + + ]
828 [ + + ]: 1072 : for (j = 0; j < 3; j++) {
829 [ + - ]: 804 : nFracSize = iter->second.fraction[j].size();
830 [ + - ]: 804 : frac_data_pointer[i] = &fractions[iFrac];
831 [ + + ]: 1170 : for (k = 0; k < nFracSize; k++) {
832 [ + - ][ + - ]: 366 : fractions[iFrac++] = iter->second.fraction[j][k];
833 : : }
834 : : }
835 : : }
836 : :
837 [ + - ][ + - ]: 1 : err = mk_core()->imesh_instance()->setIntArrData(&hvBndrHex[0], nBndrHex,
[ + - ]
838 : : m_edgeCutFracLengthTag,
839 [ + - ]: 1 : frac_leng);
840 [ + - ]: 1 : IBERRCHK(err, "Failed to set cut fraction sizes to hex.");
841 : :
842 [ + - ]: 1 : moab::ErrorCode rval = moab_instance()->tag_set_by_ptr(reinterpret_cast<moab::Tag> (m_edgeCutFracTag),
843 [ + - ]: 1 : reinterpret_cast<moab::EntityHandle*> (&hvBndrHex[0]),
844 [ + - ][ + - ]: 2 : nBndrHex, &frac_data_pointer[0], frac_size);
845 [ - + ][ # # ]: 1 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
846 : :
847 [ + - ]: 1 : delete [] frac_size;
848 [ + - ]: 1 : delete [] frac_leng;
849 [ + - ]: 2 : delete [] fractions;
850 : 1 : }
851 : :
852 : 3 : void EBMesher::fire_rays(int dir)
853 : : {
854 : : // ray fire
855 : : int i, j, k, l, index[3];
856 : 3 : double tolerance = 1e-12;
857 : 3 : double rayLength = m_nDiv[dir]*m_dIntervalSize[dir];
858 : : int iNodeStart, iNodeEnd, nIntersect, nNodeSlice;
859 : : double startPnt[3], endPnt[3];
860 : 3 : int otherDir1 = (dir + 1)%3;
861 : 3 : int otherDir2 = (dir + 2)%3;
862 : : // harmless as this expression (does nothing) so that a compiler sees it is used.
863 : : (void) iNodeEnd;
864 : : (void) nNodeSlice;
865 : :
866 [ + + ]: 30 : for (j = 1; j < m_nNode[otherDir2] - 1; j++) {
867 [ + + ]: 270 : for (i = 1; i < m_nNode[otherDir1] - 1; i++) {
868 : :
869 : : // get ray start and end points
870 [ + + ]: 243 : if (dir == 0) { // x coordinate ray
871 : 81 : iNodeStart = j*m_nNode[dir]*m_nNode[otherDir1] + i*m_nNode[dir];
872 : 81 : iNodeEnd = iNodeStart + m_nNode[dir] - 1;
873 : 81 : nNodeSlice = 1;
874 : 81 : index[0] = 0;
875 : 81 : index[1] = i;
876 : 81 : index[2] = j;
877 : : }
878 [ + + ]: 162 : else if (dir == 1) { // y coordinate ray
879 : 81 : iNodeStart = i*m_nNode[otherDir2]*m_nNode[dir] + j;
880 : 81 : iNodeEnd = iNodeStart + m_nNode[otherDir2]*(m_nNode[dir] - 1);
881 : 81 : nNodeSlice = m_nNode[otherDir2];
882 : 81 : index[0] = j;
883 : 81 : index[1] = 0;
884 : 81 : index[2] = i;
885 : : }
886 [ + - ]: 81 : else if (dir == 2) { // z coordinate ray
887 : 81 : iNodeStart = j*m_nNode[otherDir1] + i;
888 : 81 : iNodeEnd = iNodeStart + m_nNode[otherDir1]*m_nNode[otherDir2]*(m_nNode[dir] - 1);
889 : 81 : nNodeSlice = m_nNode[otherDir1]*m_nNode[otherDir2];
890 : 81 : index[0] = i;
891 : 81 : index[1] = j;
892 : 81 : index[2] = 0;
893 : : }
894 [ # # ]: 0 : else IBERRCHK(iBase_FAILURE, "Ray direction should be 0 to 2.");
895 : :
896 [ + + ]: 972 : for (l = 0; l < 3; l++) {
897 [ + + ]: 729 : if (l == dir) {
898 : 243 : startPnt[l] = m_origin_coords[l];
899 : 243 : endPnt[l] = m_origin_coords[l] + m_nDiv[dir]*m_dIntervalSize[l];
900 : : }
901 : : else {
902 : 486 : startPnt[l] = m_origin_coords[l] + index[l]*m_dIntervalSize[l];
903 : 486 : endPnt[l] = startPnt[l];
904 : : }
905 : : }
906 : :
907 : 243 : k = 0;
908 : 243 : m_iInter = 0;
909 [ + - ][ - + ]: 243 : if (!fire_ray(nIntersect, startPnt, endPnt, tolerance, dir, rayLength)) {
910 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't fire ray.");
911 : : }
912 : :
913 [ + + ]: 243 : if (nIntersect > 0) {
914 [ + - ]: 207 : m_iFirstP = m_vIntersection[m_iInter].distance/m_dIntervalSize[dir];
915 [ + - ]: 207 : m_dFirstP = startPnt[dir] + m_vIntersection[m_iInter++].distance;
916 [ + - ]: 207 : m_iSecondP = m_vIntersection[m_iInter].distance/m_dIntervalSize[dir];
917 [ + - ]: 207 : m_dSecondP = startPnt[dir] + m_vIntersection[m_iInter++].distance;
918 : :
919 : : // set outside before the first hit
920 [ + + ]: 447 : for (k = 0; k < m_iFirstP; k++) {
921 : 240 : m_nStatus = OUTSIDE;
922 : :
923 [ + - ][ - + ]: 240 : if (!set_neighbor_hex_status(dir, i, j, k)) {
924 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
925 : : }
926 : : }
927 : :
928 [ + - ]: 780 : for (; k < m_nNode[dir] - 1;) {
929 : 573 : int i_skip = 0;
930 : 573 : m_curPnt = startPnt[dir] + (k + 1)*m_dIntervalSize[dir];
931 : 573 : EdgeStatus preStat = m_nStatus;
932 [ + - ]: 573 : m_nStatus = get_edge_status(m_curPnt, i_skip);
933 [ + - ]: 573 : m_vnEdgeStatus[dir][i*m_nDiv[dir] + j*m_nDiv[dir]*m_nDiv[otherDir1] + k] = m_nStatus;
934 : :
935 : : // set status of all hexes sharing the edge
936 [ + - ][ - + ]: 573 : if (!set_neighbor_hex_status(dir, i, j, k)) {
937 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
938 : : }
939 : :
940 [ + + ]: 573 : if (m_nMove > 0) {
941 [ - + ]: 207 : if (m_iInter < nIntersect) {
942 [ # # ][ # # ]: 0 : if (!move_intersections(dir, nIntersect, startPnt)) {
943 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't move intersections.");
944 : : }
945 : : }
946 : : else {
947 : 207 : m_nMove = 1;
948 [ + - ][ + - ]: 207 : if (m_nStatus == BOUNDARY && !set_edge_status(dir)) {
[ - + ][ - + ]
949 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't set edge status.");
950 : : }
951 : 207 : k++;
952 : 207 : break; // rest is all outside
953 : : }
954 : : }
955 [ + + ][ + - ]: 366 : else if (m_nStatus == BOUNDARY && !set_edge_status(dir)) {
[ - + ][ - + ]
956 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't set edge status.");
957 : : }
958 [ - + ][ # # ]: 366 : else if (m_nStatus == OUTSIDE && preStat == BOUNDARY) { // set outside until next hit
959 : 0 : k++;
960 [ # # ]: 0 : for (; k < i_skip; k++) {
961 : 0 : m_nStatus = OUTSIDE;
962 [ # # ][ # # ]: 0 : if (!set_neighbor_hex_status(dir, i, j, k)) {
963 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
964 : : }
965 : : }
966 : : }
967 : :
968 : : // set cut-cell edge status
969 [ + + ]: 366 : if (i_skip > 0) {
970 : 207 : m_prevPnt = startPnt[dir] + i_skip*m_dIntervalSize[dir];
971 : 207 : k = i_skip;
972 : : }
973 : : else {
974 : 159 : m_prevPnt = m_curPnt;
975 : 366 : k++;
976 : : }
977 : : }
978 : : }
979 : :
980 : : // the rest are all outside
981 [ + + ]: 843 : for (; k < m_nNode[dir] - 1; k++) {
982 : 600 : m_nStatus = OUTSIDE;
983 [ + - ][ - + ]: 600 : if (!set_neighbor_hex_status(dir, i, j, k)) {
984 [ # # ]: 0 : IBERRCHK(iBase_FAILURE, "Couldn't set neighbor hex status.");
985 : : }
986 : : }
987 : : }
988 : : }
989 : 3 : }
990 : :
991 : 243 : bool EBMesher::fire_ray(int& nIntersect, double startPnt[3],
992 : : double endPnt[3], double tol, int dir,
993 : : double rayLength)
994 : : {
995 : 243 : m_vIntersection.clear();
996 : 243 : m_vhInterSurf.clear();
997 : 243 : m_vhInterFacet.clear();
998 : 243 : m_mhOverlappedSurf.clear();
999 [ + - ]: 243 : std::vector<double> temp_intersects;
1000 : : moab::ErrorCode rVal;
1001 [ + - ]: 486 : moab::OrientedBoxTreeTool::TrvStats stats;
1002 [ + - ]: 243 : if (m_bUseGeom) { // geometry input
1003 : : rVal = m_hObbTree->ray_intersect_sets(temp_intersects, m_vhInterSurf,
1004 : : m_vhInterFacet, m_hTreeRoot, tol,
1005 : : startPnt, rayDir[dir], &rayLength,
1006 [ + - ]: 243 : &stats);
1007 : : }
1008 : : else { // facet input
1009 [ # # ]: 0 : std::vector<moab::EntityHandle> dum_facets_out;
1010 : : rVal = m_hObbTree->ray_intersect_triangles(temp_intersects, dum_facets_out,
1011 : : m_hTreeRoot, tol,
1012 : : startPnt, rayDir[dir], &rayLength,
1013 [ # # ]: 0 : &stats);
1014 : : }
1015 : :
1016 : 243 : nIntersect = temp_intersects.size();
1017 [ - + ]: 243 : if (moab::MB_SUCCESS != rVal) {
1018 [ # # ][ # # ]: 0 : std::cerr << "Failed : ray-triangle intersection." << std::endl;
1019 : 0 : return false;
1020 : : }
1021 : :
1022 [ + + ]: 243 : if (nIntersect > 0) {
1023 [ + - ]: 207 : m_vIntersection.resize(nIntersect);
1024 : :
1025 [ + + ]: 730 : for (int l = 0; l < nIntersect; l++) {
1026 [ + - ][ + - ]: 523 : IntersectDist temp_inter_dist(temp_intersects[l], l);
1027 [ + - ]: 523 : m_vIntersection[l] = temp_inter_dist;
1028 : : }
1029 [ + - ]: 207 : std::sort(m_vIntersection.begin(), m_vIntersection.end(), less_intersect);
1030 : :
1031 [ + - ]: 207 : if (nIntersect > 1) {
1032 : : bool bMoveOnce;
1033 : 207 : m_nIteration = 0;
1034 : 207 : m_iOverlap = 0;
1035 : :
1036 [ + - ]: 207 : if (m_bUseGeom) { // if there are geometry info
1037 [ + - ]: 207 : remove_intersection_duplicates();
1038 : : }
1039 [ # # ][ # # ]: 0 : else if (is_ray_move_and_set_overlap_surf(bMoveOnce)) { // facet geom case
1040 [ # # ][ # # ]: 0 : if (!move_ray(nIntersect, startPnt, endPnt, tol, dir, bMoveOnce)) {
1041 [ # # ][ # # ]: 0 : std::cerr << "Number of Intersection between edges and ray should be even." << std::endl;
1042 : 0 : return false;
1043 : : }
1044 : : }
1045 : 207 : nIntersect = m_vIntersection.size();
1046 : : }
1047 : : }
1048 : :
1049 : 486 : return true;
1050 : : }
1051 : :
1052 : 0 : bool EBMesher::move_intersections(int n_dir, int n_inter, double start_pnt[3])
1053 : : {
1054 [ # # ]: 0 : if (m_nMove > 0) {
1055 [ # # ]: 0 : if (m_iInter < n_inter) {
1056 [ # # ]: 0 : if (m_nMove == 1) {
1057 [ # # ]: 0 : do {
1058 [ # # ][ # # ]: 0 : if (m_nStatus == BOUNDARY && !set_edge_status(n_dir)) return false;
[ # # ]
1059 : 0 : m_iFirstP = m_iSecondP;
1060 : 0 : m_dFirstP = m_dSecondP;
1061 : 0 : m_iSecondP = m_vIntersection[m_iInter].distance/m_dIntervalSize[n_dir];
1062 : 0 : m_dSecondP = start_pnt[n_dir] + m_vIntersection[m_iInter++].distance;
1063 : : }
1064 [ # # ]: 0 : while (m_dSecondP < m_curPnt && m_iInter < n_inter);
1065 : : }
1066 [ # # ]: 0 : else if (m_nMove == 2) {
1067 [ # # ]: 0 : do {
1068 : 0 : m_iFirstP = m_vIntersection[m_iInter].distance/m_dIntervalSize[n_dir];
1069 : 0 : m_dFirstP = start_pnt[n_dir] + m_vIntersection[m_iInter++].distance;
1070 [ # # ][ # # ]: 0 : if (m_nStatus == BOUNDARY && !set_edge_status(n_dir)) return false;
[ # # ]
1071 [ # # ]: 0 : if (m_iInter < n_inter) {
1072 : 0 : m_iSecondP = m_vIntersection[m_iInter].distance/m_dIntervalSize[n_dir];
1073 : 0 : m_dSecondP = start_pnt[n_dir] + m_vIntersection[m_iInter++].distance;
1074 : : }
1075 : : else {
1076 : 0 : m_iSecondP = m_iFirstP;
1077 : 0 : m_dSecondP = m_dFirstP;
1078 : : }
1079 : : }
1080 [ # # ]: 0 : while (m_dSecondP < m_curPnt && m_iInter < n_inter);
1081 : : }
1082 : : }
1083 : : }
1084 : :
1085 : 0 : return true;
1086 : : }
1087 : :
1088 : 207 : bool EBMesher::is_shared_overlapped_surf(int index)
1089 : : {
1090 : : int nParent;
1091 : : iBase_ErrorType err;
1092 : : moab::EntityHandle hSurf;
1093 [ + - ]: 207 : if (m_bUseGeom) {
1094 [ + - ][ + - ]: 207 : hSurf = m_vhInterSurf[m_vIntersection[index].index];
1095 [ + - ][ + - ]: 207 : err = mk_core()->imesh_instance()->getNumPrnt(reinterpret_cast<iBase_EntitySetHandle> (hSurf),
1096 [ + - ]: 207 : 1, nParent);
1097 [ + - ]: 207 : IBERRCHK(err, "Failed to get number of surface parents.\n");
1098 : :
1099 [ - + ]: 207 : if (nParent > 1) return true;
1100 : : }
1101 [ # # ]: 0 : else hSurf = m_vIntersection[index].index;
1102 : :
1103 [ + - ]: 207 : return m_mhOverlappedSurf.count(hSurf) > 0;
1104 : : }
1105 : :
1106 : 1 : void EBMesher::get_grid_and_edges_techX(double* boxMin, double* boxMax, int* nDiv,
1107 : : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellSurfEdge,
1108 : : std::vector<int>& rvnInsideCell, bool isCornerExterior)
1109 : : {
1110 : : int i, j, err, ii, jj, kk, iHex;
1111 [ + + ]: 4 : for (i = 0; i < 3; i++) {
1112 : 3 : boxMin[i] = m_origin_coords[i];
1113 : 3 : boxMax[i] = m_origin_coords[i] + m_dIntervalSize[i]*m_nDiv[i];
1114 : 3 : nDiv[i] = m_nDiv[i];
1115 : : }
1116 : :
1117 : : // get all hexes
1118 [ + - ]: 1 : std::vector<iBase_EntityHandle> hex_handles;
1119 : : err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
1120 [ + - ][ + - ]: 1 : iMesh_HEXAHEDRON, hex_handles);
[ + - ]
1121 [ + - ]: 1 : IBERRCHK(err, "Failed to get hexes.\n");
1122 : :
1123 : : // get hex status
1124 : 1 : int hex_size = hex_handles.size();
1125 : : /*
1126 : : //int* hex_status = new int[hex_size];
1127 : : //int* hex_status;
1128 : : //std::vector<int> hex_status(hex_size);
1129 : : int temp = 0;
1130 : : int* hex_status = &temp;
1131 : : err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
1132 : : m_elemStatusTag, hex_status);
1133 : : */
1134 : : int error;
1135 [ + - ][ + - ]: 1 : int* hex_status = new int[hex_size];
1136 : 1 : int hex_status_alloc = sizeof(int)*hex_size;
1137 : 1 : int hex_status_size = 0;
1138 [ + - ]: 1 : iMesh_getIntArrData(m_mesh, &hex_handles[0],
1139 : : hex_size, m_elemStatusTag,
1140 : : &hex_status, &hex_status_alloc,
1141 [ + - ]: 1 : &hex_status_size, &error);
1142 [ + - ]: 1 : IBERRCHK(error, "Failed to get hex status.");
1143 : :
1144 : : // get inside and boundary hexes
1145 : 1 : int nInside = 0;
1146 : 1 : int nOutside = 0;
1147 : 1 : int nBoundary = 0;
1148 : 1 : rvnInsideCell.clear();
1149 [ + + ]: 1001 : for (i = 0; i < hex_size; i++) {
1150 [ + + ]: 1000 : if (hex_status[i] == 0) { // if inside
1151 [ + - ]: 304 : iHex = moab_instance()->id_from_handle(reinterpret_cast<moab::EntityHandle>
1152 [ + - ][ + - ]: 304 : (hex_handles[i])) - m_iStartHex;
1153 [ + - ]: 304 : rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0]);
1154 [ + - ]: 304 : rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0]);
1155 [ + - ]: 304 : rvnInsideCell.push_back(iHex/m_nDiv[0]/m_nDiv[1]);
1156 : 304 : nInside++;
1157 : : }
1158 [ + + ]: 696 : else if (hex_status[i] == 1) nOutside++;
1159 [ + - ]: 392 : else if (hex_status[i] == 2) nBoundary++;
1160 [ # # ]: 0 : else IBERRCHK(err, "Element status should be one of inside/outside/boundary.");
1161 : : }
1162 [ + - ][ + - ]: 1 : std::cout << "# of inside, outside, boundary elem: " << nInside
1163 [ + - ][ + - ]: 1 : << ", " << nOutside << ", " << nBoundary << std::endl;
[ + - ][ + - ]
[ + - ]
1164 [ + - ]: 1 : delete [] hex_status;
1165 : :
1166 : : // get cut-cell fractions
1167 : : double dFrac[4];
1168 : 1 : std::map<int, CutFraction>::iterator iter = m_mdCutFraction.begin();
1169 : 1 : std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
1170 : :
1171 [ + - ][ + - ]: 269 : for (; iter != end_iter; iter++) { // for each cut-cell
[ + + ]
1172 : : // get i, j, k index from handle
1173 [ + - ]: 268 : iHex = iter->first;
1174 : 268 : ii = (iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0];
1175 : 268 : jj = (iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0];
1176 : 268 : kk = iHex/m_nDiv[0]/m_nDiv[1];
1177 : :
1178 : : // surface 1
1179 [ + - ][ + - ]: 268 : CutFraction cutFrac = iter->second;
1180 [ + + ]: 268 : if (cutFrac.fraction[1].size() > 0) {
1181 [ + - ]: 122 : dFrac[0] = cutFrac.fraction[1][0];
1182 [ + + ]: 122 : if (dFrac[0] < 0.) dFrac[0] *= -1.;
1183 : : }
1184 : 146 : else dFrac[0] = -1.;
1185 [ + + ]: 268 : if (cutFrac.fraction[2].size() > 0) {
1186 [ + - ]: 122 : dFrac[3] = cutFrac.fraction[2][0];
1187 [ + + ]: 122 : if (dFrac[3] < 0.) dFrac[3] *= -1.;
1188 : : }
1189 : 146 : else dFrac[3] = -1.;
1190 [ + - ]: 268 : dFrac[1] = get_edge_fraction(iHex + m_nDiv[0], 2);
1191 [ + - ]: 268 : dFrac[2] = get_edge_fraction(iHex + m_nDiv[0]*m_nDiv[1], 1);
1192 : :
1193 [ + + ][ + + ]: 268 : if (dFrac[0] > 0. || dFrac[1] > 0. || dFrac[2] > 0. || dFrac[3] > 0.) { // if surface is cut
[ + + ][ + + ]
1194 [ + - ]: 211 : CutCellSurfEdgeKey surfKey(ii, jj, kk, 0);
1195 [ + - ]: 211 : std::vector<double> edges;
1196 [ + + ][ + - ]: 211 : if (dFrac[0] < 0.) dFrac[0] = get_uncut_edge_fraction(ii, jj, kk, 1);
1197 [ + + ][ + - ]: 211 : if (dFrac[1] < 0.) dFrac[1] = get_uncut_edge_fraction(ii, jj + 1, kk, 2);
1198 [ + + ][ + - ]: 211 : if (dFrac[2] < 0.) dFrac[2] = get_uncut_edge_fraction(ii, jj, kk + 1, 1);
1199 [ + + ][ + - ]: 211 : if (dFrac[3] < 0.) dFrac[3] = get_uncut_edge_fraction(ii, jj, kk, 2);
1200 [ + - ][ + + ]: 1055 : for (j = 0; j < 4; j++) edges.push_back(dFrac[j]);
1201 [ + - ][ + - ]: 211 : rmdCutCellSurfEdge[surfKey] = edges;
1202 : : }
1203 : :
1204 : : // surface 2
1205 [ + - ][ + - ]: 268 : cutFrac = iter->second;
1206 [ + + ]: 268 : if (cutFrac.fraction[0].size() > 0) {
1207 [ + - ]: 122 : dFrac[0] = cutFrac.fraction[0][0];
1208 [ + + ]: 122 : if (dFrac[0] < 0.) dFrac[0] *= -1.;
1209 : : }
1210 : 146 : else dFrac[0] = -1.;
1211 [ + + ]: 268 : if (cutFrac.fraction[1].size() > 0) {
1212 [ + - ]: 122 : dFrac[3] = cutFrac.fraction[1][0];
1213 [ + + ]: 122 : if (dFrac[3] < 0.) dFrac[3] *= -1.;
1214 : : }
1215 : 146 : else dFrac[3] = -1.;
1216 [ + - ]: 268 : dFrac[1] = get_edge_fraction(iHex + 1, 2);
1217 [ + - ]: 268 : dFrac[2] = get_edge_fraction(iHex + m_nDiv[0]*m_nDiv[1], 0);
1218 : :
1219 [ + + ][ + + ]: 268 : if (dFrac[0] > 0. || dFrac[1] > 0. || dFrac[2] > 0. || dFrac[3] > 0.) { // if surface is cut
[ + + ][ + + ]
1220 [ + - ]: 263 : CutCellSurfEdgeKey surfKey(ii, jj, kk, 1);
1221 [ + - ]: 263 : std::vector<double> edges;
1222 [ + + ][ + - ]: 263 : if (dFrac[0] < 0.) dFrac[0] = get_uncut_edge_fraction(ii, jj, kk, 0);
1223 [ + + ][ + - ]: 263 : if (dFrac[1] < 0.) dFrac[1] = get_uncut_edge_fraction(ii + 1, jj, kk, 2);
1224 [ + + ][ + - ]: 263 : if (dFrac[2] < 0.) dFrac[2] = get_uncut_edge_fraction(ii, jj, kk + 1, 0);
1225 [ + + ][ + - ]: 263 : if (dFrac[3] < 0.) dFrac[3] = get_uncut_edge_fraction(ii, jj, kk, 2);
1226 [ + - ][ + + ]: 1315 : for (j = 0; j < 4; j++) edges.push_back(dFrac[j]);
1227 [ + - ][ + - ]: 263 : rmdCutCellSurfEdge[surfKey] = edges;
1228 : : }
1229 : :
1230 : : // surface 3
1231 [ + - ][ + - ]: 268 : cutFrac = iter->second;
1232 [ + + ]: 268 : if (cutFrac.fraction[0].size() > 0) {
1233 [ + - ]: 122 : dFrac[0] = cutFrac.fraction[0][0];
1234 [ + + ]: 122 : if (dFrac[0] < 0.) dFrac[0] *= -1.;
1235 : : }
1236 : 146 : else dFrac[0] = -1.;
1237 [ + + ]: 268 : if (cutFrac.fraction[1].size() > 0) {
1238 [ + - ]: 122 : dFrac[3] = cutFrac.fraction[1][0];
1239 [ + + ]: 122 : if (dFrac[3] < 0.) dFrac[3] *= -1.;
1240 : : }
1241 : 146 : else dFrac[3] = -1.;
1242 [ + - ]: 268 : dFrac[1] = get_edge_fraction(iHex + 1, 1);
1243 [ + - ]: 268 : dFrac[2] = get_edge_fraction(iHex + m_nDiv[0], 0);
1244 : :
1245 [ + + ][ + + ]: 268 : if (dFrac[0] > 0. || dFrac[1] > 0. || dFrac[2] > 0. || dFrac[3] > 0.) { // if surface is cut
[ + + ][ + + ]
1246 [ + - ]: 211 : CutCellSurfEdgeKey surfKey(ii, jj, kk, 2);
1247 [ + - ]: 211 : std::vector<double> edges;
1248 [ + + ][ + - ]: 211 : if (dFrac[0] < 0.) dFrac[0] = get_uncut_edge_fraction(ii, jj, kk, 0);
1249 [ + + ][ + - ]: 211 : if (dFrac[1] < 0.) dFrac[1] = get_uncut_edge_fraction(ii + 1, jj, kk, 1);
1250 [ + + ][ + - ]: 211 : if (dFrac[2] < 0.) dFrac[2] = get_uncut_edge_fraction(ii, jj + 1, kk, 0);
1251 [ + + ][ + - ]: 211 : if (dFrac[3] < 0.) dFrac[3] = get_uncut_edge_fraction(ii, jj, kk, 1);
1252 [ + - ][ + + ]: 1055 : for (j = 0; j < 4; j++) edges.push_back(dFrac[j]);
1253 [ + - ][ + - ]: 211 : rmdCutCellSurfEdge[surfKey] = edges;
1254 : : }
1255 : 268 : }
1256 : :
1257 : 1 : if (debug_ebmesh) export_fraction_edges(rmdCutCellSurfEdge);
1258 : 1 : }
1259 : :
1260 : 1 : bool EBMesher::get_grid_and_edges(double* boxMin, double* boxMax, int* nDiv,
1261 : : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& rmdCutCellEdge,
1262 : : std::vector<int>& rvnInsideCell, bool isCornerExterior)
1263 : : {
1264 : : int i, ii, jj, kk, iHex;
1265 [ + + ]: 4 : for (i = 0; i < 3; i++) {
1266 : 3 : boxMin[i] = m_origin_coords[i];
1267 : 3 : boxMax[i] = m_origin_coords[i] + m_dIntervalSize[i]*m_nDiv[i];
1268 : 3 : nDiv[i] = m_nDiv[i];
1269 : : }
1270 : :
1271 [ + - ][ - + ]: 1 : if (!get_inside_hex(rvnInsideCell)) return false;
1272 : :
1273 : : // get cut-cell fractions
1274 : 1 : std::map<int, CutFraction>::iterator iter = m_mdCutFraction.begin();
1275 : 1 : std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
1276 : :
1277 [ + - ][ + - ]: 269 : for (; iter != end_iter; iter++) { // for each cut-cell
[ + + ]
1278 : : // get i, j, k index from handle
1279 [ + - ]: 268 : iHex = iter->first;
1280 : 268 : ii = (iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0];
1281 : 268 : jj = (iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0];
1282 : 268 : kk = iHex/m_nDiv[0]/m_nDiv[1];
1283 : :
1284 [ + + ]: 1072 : for (i = 0; i < 3; i++) {
1285 [ + - ][ + - ]: 804 : std::vector<double> fractions(iter->second.fraction[i]);
1286 [ + - ]: 804 : CutCellSurfEdgeKey edgeKey(ii, jj, kk, i);
1287 [ + - ][ + - ]: 804 : rmdCutCellEdge[edgeKey] = fractions;
1288 : 804 : }
1289 : : }
1290 : :
1291 : : if (debug_ebmesh) return export_fraction_points(rmdCutCellEdge);
1292 : :
1293 : 1 : return true;
1294 : : }
1295 : :
1296 : 1 : bool EBMesher::get_inside_hex(std::vector<int>& rvnInsideCell)
1297 : : {
1298 : : int i, err, iHex;
1299 : :
1300 : : // get all hexes
1301 [ + - ]: 1 : std::vector<iBase_EntityHandle> hex_handles;
1302 : : err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
1303 [ + - ][ + - ]: 1 : iMesh_HEXAHEDRON, hex_handles);
[ + - ]
1304 [ + - ]: 1 : IBERRCHK(err, "Failed to get hexes.\n");
1305 : :
1306 : : // get hex status
1307 : 1 : int hex_size = hex_handles.size();
1308 : : /*
1309 : : //int* hex_status = new int[hex_size];
1310 : : //int* hex_status;
1311 : : std::vector<int> hex_status(hex_size);
1312 : : err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
1313 : : m_elemStatusTag, &hex_status[0]);
1314 : : */
1315 : : int error;
1316 [ + - ][ + - ]: 1 : int* hex_status = new int[hex_size];
1317 : 1 : int hex_status_alloc = sizeof(int)*hex_size;
1318 : 1 : int hex_status_size = 0;
1319 [ + - ]: 1 : iMesh_getIntArrData(m_mesh, &hex_handles[0],
1320 : : hex_size, m_elemStatusTag,
1321 : : &hex_status, &hex_status_alloc,
1322 [ + - ]: 1 : &hex_status_size, &error);
1323 [ - + ][ # # ]: 1 : ERRORRF("Failed to get hex status.\n");
[ # # ]
1324 : :
1325 : : // get inside and boundary hexes
1326 : 1 : int nInside = 0;
1327 : 1 : int nOutside = 0;
1328 : 1 : int nBoundary = 0;
1329 : 1 : rvnInsideCell.clear();
1330 [ + + ]: 1001 : for (i = 0; i < hex_size; i++) {
1331 [ + + ]: 1000 : if (hex_status[i] == 0) { // if inside
1332 [ + - ]: 304 : iHex = moab_instance()->id_from_handle(reinterpret_cast<moab::EntityHandle>
1333 [ + - ][ + - ]: 304 : (hex_handles[i])) - m_iStartHex;
1334 [ + - ]: 304 : rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0]);
1335 [ + - ]: 304 : rvnInsideCell.push_back((iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0]);
1336 [ + - ]: 304 : rvnInsideCell.push_back(iHex/m_nDiv[0]/m_nDiv[1]);
1337 : 304 : nInside++;
1338 : : }
1339 [ + + ]: 696 : else if (hex_status[i] == 1) nOutside++;
1340 [ + - ]: 392 : else if (hex_status[i] == 2) nBoundary++;
1341 [ # # ][ # # ]: 0 : else ERRORRF("Element status should be one of inside/outside/boundary.\n");
[ # # ]
1342 : : }
1343 [ + - ][ + - ]: 1 : std::cout << "# of inside, outside, boundary elem: " << nInside
1344 [ + - ][ + - ]: 1 : << ", " << nOutside << ", " << nBoundary << std::endl;
[ + - ][ + - ]
[ + - ]
1345 : :
1346 [ + - ]: 1 : delete [] hex_status;
1347 : 1 : return true;
1348 : : }
1349 : :
1350 : 0 : bool EBMesher::export_fraction_edges(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellSurfEdge)
1351 : : {
1352 : : // export fractions as edge
1353 : : double curPnt[3], ePnt[6];
1354 [ # # ]: 0 : std::vector<iBase_EntityHandle> edge_handles;
1355 : 0 : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator itr = mdCutCellSurfEdge.begin();
1356 : 0 : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator e_itr = mdCutCellSurfEdge.end();
1357 [ # # ][ # # ]: 0 : for (; itr != e_itr; itr++) {
[ # # ]
1358 [ # # ]: 0 : curPnt[0] = m_origin_coords[0] + itr->first.i*m_dIntervalSize[0];
1359 [ # # ]: 0 : curPnt[1] = m_origin_coords[1] + itr->first.j*m_dIntervalSize[1];
1360 [ # # ]: 0 : curPnt[2] = m_origin_coords[2] + itr->first.k*m_dIntervalSize[2];
1361 [ # # ][ # # ]: 0 : std::vector<double> edges = itr->second;
1362 : :
1363 [ # # ][ # # ]: 0 : if (itr->first.l == 0) {
1364 [ # # ][ # # ]: 0 : if (edges[0] > 0. && edges[0] < 1.) {
[ # # ][ # # ]
[ # # ]
1365 : 0 : ePnt[0] = curPnt[0];
1366 : 0 : ePnt[1] = curPnt[1];
1367 : 0 : ePnt[2] = curPnt[2];
1368 : 0 : ePnt[3] = ePnt[0];
1369 [ # # ]: 0 : ePnt[4] = ePnt[1] + edges[0]*m_dIntervalSize[1];
1370 : 0 : ePnt[5] = ePnt[2];
1371 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1372 : : }
1373 [ # # ][ # # ]: 0 : if (edges[1] > 0. && edges[1] < 1.) {
[ # # ][ # # ]
[ # # ]
1374 : 0 : ePnt[0] = curPnt[0];
1375 : 0 : ePnt[1] = curPnt[1] + m_dIntervalSize[1];
1376 : 0 : ePnt[2] = curPnt[2];
1377 : 0 : ePnt[3] = ePnt[0];
1378 : 0 : ePnt[4] = ePnt[1];
1379 [ # # ]: 0 : ePnt[5] = ePnt[2] + edges[1]*m_dIntervalSize[2];
1380 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1381 : : }
1382 [ # # ][ # # ]: 0 : if (edges[2] > 0. && edges[2] < 1.) {
[ # # ][ # # ]
[ # # ]
1383 : 0 : ePnt[0] = curPnt[0];
1384 : 0 : ePnt[1] = curPnt[1];
1385 : 0 : ePnt[2] = curPnt[2] + m_dIntervalSize[2];
1386 : 0 : ePnt[3] = ePnt[0];
1387 [ # # ]: 0 : ePnt[4] = ePnt[1] + edges[2]*m_dIntervalSize[1];
1388 : 0 : ePnt[5] = ePnt[2];
1389 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1390 : : }
1391 [ # # ][ # # ]: 0 : if (edges[3] > 0. && edges[3] < 1.) {
[ # # ][ # # ]
[ # # ]
1392 : 0 : ePnt[0] = curPnt[0];
1393 : 0 : ePnt[1] = curPnt[1];
1394 : 0 : ePnt[2] = curPnt[2];
1395 : 0 : ePnt[3] = ePnt[0];
1396 : 0 : ePnt[4] = ePnt[1];
1397 [ # # ]: 0 : ePnt[5] = ePnt[2] + edges[3]*m_dIntervalSize[2];
1398 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1399 : : }
1400 : : }
1401 [ # # ][ # # ]: 0 : if (itr->first.l == 1) {
1402 [ # # ][ # # ]: 0 : if (edges[0] > 0. && edges[0] < 1.) {
[ # # ][ # # ]
[ # # ]
1403 : 0 : ePnt[0] = curPnt[0];
1404 : 0 : ePnt[1] = curPnt[1];
1405 : 0 : ePnt[2] = curPnt[2];
1406 [ # # ]: 0 : ePnt[3] = ePnt[0] + edges[0]*m_dIntervalSize[0];
1407 : 0 : ePnt[4] = ePnt[1];
1408 : 0 : ePnt[5] = ePnt[2];
1409 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1410 : : }
1411 [ # # ][ # # ]: 0 : if (edges[1] > 0. && edges[1] < 1.) {
[ # # ][ # # ]
[ # # ]
1412 : 0 : ePnt[0] = curPnt[0] + m_dIntervalSize[0];
1413 : 0 : ePnt[1] = curPnt[1];
1414 : 0 : ePnt[2] = curPnt[2];
1415 : 0 : ePnt[3] = ePnt[0];
1416 : 0 : ePnt[4] = ePnt[1];
1417 [ # # ]: 0 : ePnt[5] = ePnt[2] + edges[1]*m_dIntervalSize[2];
1418 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1419 : : }
1420 [ # # ][ # # ]: 0 : if (edges[2] > 0. && edges[2] < 1.) {
[ # # ][ # # ]
[ # # ]
1421 : 0 : ePnt[0] = curPnt[0];
1422 : 0 : ePnt[1] = curPnt[1];
1423 : 0 : ePnt[2] = curPnt[2] + m_dIntervalSize[2];
1424 [ # # ]: 0 : ePnt[3] = ePnt[0] + edges[2]*m_dIntervalSize[0];
1425 : 0 : ePnt[4] = ePnt[1];
1426 : 0 : ePnt[5] = ePnt[2];
1427 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1428 : : }
1429 [ # # ][ # # ]: 0 : if (edges[3] > 0. && edges[3] < 1.) {
[ # # ][ # # ]
[ # # ]
1430 : 0 : ePnt[0] = curPnt[0];
1431 : 0 : ePnt[1] = curPnt[1];
1432 : 0 : ePnt[2] = curPnt[2];
1433 : 0 : ePnt[3] = ePnt[0];
1434 : 0 : ePnt[4] = ePnt[1];
1435 [ # # ]: 0 : ePnt[5] = ePnt[2] + edges[3]*m_dIntervalSize[2];
1436 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1437 : : }
1438 : : }
1439 [ # # ][ # # ]: 0 : if (itr->first.l == 2) {
1440 [ # # ][ # # ]: 0 : if (edges[0] > 0. && edges[0] < 1.) {
[ # # ][ # # ]
[ # # ]
1441 : 0 : ePnt[0] = curPnt[0];
1442 : 0 : ePnt[1] = curPnt[1];
1443 : 0 : ePnt[2] = curPnt[2];
1444 [ # # ]: 0 : ePnt[3] = ePnt[0] + edges[0]*m_dIntervalSize[0];
1445 : 0 : ePnt[4] = ePnt[1];
1446 : 0 : ePnt[5] = ePnt[2];
1447 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1448 : : }
1449 [ # # ][ # # ]: 0 : if (edges[1] > 0. && edges[1] < 1.) {
[ # # ][ # # ]
[ # # ]
1450 : 0 : ePnt[0] = curPnt[0] + m_dIntervalSize[0];
1451 : 0 : ePnt[1] = curPnt[1];
1452 : 0 : ePnt[2] = curPnt[2];
1453 : 0 : ePnt[3] = ePnt[0];
1454 [ # # ]: 0 : ePnt[4] = ePnt[1] + edges[1]*m_dIntervalSize[1];
1455 : 0 : ePnt[5] = ePnt[2];
1456 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1457 : : }
1458 [ # # ][ # # ]: 0 : if (edges[2] > 0. && edges[2] < 1.) {
[ # # ][ # # ]
[ # # ]
1459 : 0 : ePnt[0] = curPnt[0];
1460 : 0 : ePnt[1] = curPnt[1] + m_dIntervalSize[2];
1461 : 0 : ePnt[2] = curPnt[2];
1462 [ # # ]: 0 : ePnt[3] = ePnt[0] + edges[2]*m_dIntervalSize[0];
1463 : 0 : ePnt[4] = ePnt[1];
1464 : 0 : ePnt[5] = ePnt[2];
1465 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
1466 : : }
1467 [ # # ][ # # ]: 0 : if (edges[3] > 0. && edges[3] < 1.) {
[ # # ][ # # ]
[ # # ]
1468 : 0 : ePnt[0] = curPnt[0];
1469 : 0 : ePnt[1] = curPnt[1];
1470 : 0 : ePnt[2] = curPnt[2];
1471 : 0 : ePnt[3] = ePnt[0];
1472 [ # # ]: 0 : ePnt[4] = ePnt[1] + edges[3]*m_dIntervalSize[1];
1473 : 0 : ePnt[5] = ePnt[2];
1474 [ # # ][ # # ]: 0 : if (!make_edge(ePnt, edge_handles)) return false;
[ # # ]
1475 : : }
1476 : : }
1477 : 0 : }
1478 : :
1479 : 0 : int is_list = 1;
1480 : : iBase_EntitySetHandle set;
1481 : : iBase_ErrorType err;
1482 [ # # ][ # # ]: 0 : err = mk_core()->imesh_instance()->createEntSet(is_list, set);
[ # # ]
1483 [ # # ]: 0 : IBERRCHK(err, "Couldn't create set.");
1484 : :
1485 [ # # ][ # # ]: 0 : err = mk_core()->imesh_instance()->addEntArrToSet(&edge_handles[0],
[ # # ]
1486 [ # # ]: 0 : edge_handles.size(), set);
1487 [ # # ]: 0 : IBERRCHK(err, "Couldn't add edges to set.");
1488 : :
1489 [ # # ]: 0 : moab::ErrorCode rval = moab_instance()->write_mesh("edges.vtk",
1490 [ # # ]: 0 : (const moab::EntityHandle*) &set, 1);
1491 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
1492 : :
1493 : 0 : return true;
1494 : : }
1495 : :
1496 : 0 : bool EBMesher::export_fraction_points(std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >& mdCutCellEdge)
1497 : : {
1498 : : // export fractions as edge
1499 : : int i, j, dir, nFrc;
1500 : : iBase_ErrorType err;
1501 : : double curPnt[3], fracPnt[3], frac;
1502 : : iBase_EntityHandle v_handle;
1503 [ # # ]: 0 : std::vector<iBase_EntityHandle> vertex_handles;
1504 : 0 : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator itr = mdCutCellEdge.begin();
1505 : 0 : std::map< CutCellSurfEdgeKey, std::vector<double>, LessThan >::iterator e_itr = mdCutCellEdge.end();
1506 [ # # ][ # # ]: 0 : for (; itr != e_itr; itr++) {
[ # # ]
1507 [ # # ]: 0 : curPnt[0] = m_origin_coords[0] + itr->first.i*m_dIntervalSize[0];
1508 [ # # ]: 0 : curPnt[1] = m_origin_coords[1] + itr->first.j*m_dIntervalSize[1];
1509 [ # # ]: 0 : curPnt[2] = m_origin_coords[2] + itr->first.k*m_dIntervalSize[2];
1510 [ # # ]: 0 : dir = itr->first.l;
1511 [ # # ]: 0 : nFrc = itr->second.size();
1512 : :
1513 [ # # ]: 0 : for (i = 0; i < nFrc; i++) {
1514 [ # # ][ # # ]: 0 : frac = itr->second[i];
1515 [ # # ][ # # ]: 0 : if (frac > 0. && frac < 1.) {
1516 [ # # ]: 0 : for (j = 0; j < 3; j++) {
1517 [ # # ]: 0 : if (j == dir) fracPnt[j] = curPnt[j] + frac*m_dIntervalSize[j];
1518 : 0 : else fracPnt[j] = curPnt[j];
1519 : : }
1520 : : err = mk_core()->imesh_instance()->createVtx(fracPnt[0], fracPnt[1],
1521 [ # # ][ # # ]: 0 : fracPnt[2], v_handle);
[ # # ]
1522 [ # # ]: 0 : IBERRCHK(err, "Couldn't create vertex.");
1523 [ # # ]: 0 : vertex_handles.push_back(v_handle);
1524 : : }
1525 : : }
1526 : : }
1527 : :
1528 : 0 : int is_list = 1;
1529 : : moab::ErrorCode result;
1530 : : iBase_EntitySetHandle set;
1531 [ # # ][ # # ]: 0 : err = mk_core()->imesh_instance()->createEntSet(is_list, set);
[ # # ]
1532 [ # # ]: 0 : IBERRCHK(err, "Couldn't create set.");
1533 : :
1534 [ # # ][ # # ]: 0 : err = mk_core()->imesh_instance()->addEntArrToSet(&vertex_handles[0],
[ # # ]
1535 [ # # ]: 0 : vertex_handles.size(), set);
1536 [ # # ]: 0 : IBERRCHK(err, "Couldn't add vertices to set.");
1537 : :
1538 [ # # ]: 0 : result = moab_instance()->write_mesh("frac_vertices.vtk",
1539 [ # # ]: 0 : (const moab::EntityHandle*) &set, 1);
1540 [ # # ]: 0 : if (moab::MB_SUCCESS != result) {
1541 [ # # ][ # # ]: 0 : std::cerr << "Failed to write fraction vertices." << std::endl;
1542 : 0 : return false;
1543 : : }
1544 : :
1545 : 0 : return true;
1546 : : }
1547 : :
1548 : 0 : bool EBMesher::make_edge(double ePnt[6], std::vector<iBase_EntityHandle>& edge_handles)
1549 : : {
1550 : : iBase_ErrorType err;
1551 : : iBase_EntityHandle vertex_handle[2], edge_handle;
1552 : 0 : iBase_EntityHandle* pVertexHandle = &vertex_handle[0];
1553 : : err = mk_core()->imesh_instance()->createVtxArr(2, iBase_INTERLEAVED,
1554 [ # # ][ # # ]: 0 : ePnt, pVertexHandle);
[ # # ]
1555 [ # # ]: 0 : IBERRCHK(err, "Failed to create vertices.\n");
1556 : :
1557 : : err = mk_core()->imesh_instance()->createEnt(iMesh_LINE_SEGMENT,
1558 : : &vertex_handle[0], 2,
1559 [ # # ][ # # ]: 0 : edge_handle);
[ # # ]
1560 [ # # ]: 0 : IBERRCHK(err, "Failed to create edge.\n");
1561 : :
1562 [ # # ]: 0 : edge_handles.push_back(edge_handle);
1563 : :
1564 : 0 : return true;
1565 : : }
1566 : :
1567 : 1608 : double EBMesher::get_edge_fraction(int idHex, int dir)
1568 : : {
1569 : 1608 : std::map<int, CutFraction>::iterator end_iter = m_mdCutFraction.end();
1570 [ + - ]: 1608 : std::map<int, CutFraction>::iterator iter = m_mdCutFraction.find(idHex);
1571 [ + - ][ + + ]: 1608 : if (iter != end_iter && iter->second.fraction[dir].size() > 0) {
[ + - ][ + + ]
[ + + ]
1572 [ + - ][ + - ]: 516 : double frac = iter->second.fraction[dir][0];
1573 [ + + ]: 516 : if (frac < 0.) frac *= -1.;
1574 : 516 : return frac;
1575 : : }
1576 : 1608 : else return -1.;
1577 : : }
1578 : :
1579 : 1492 : double EBMesher::get_uncut_edge_fraction(int i, int j, int k, int dir)
1580 : : {
1581 : : int iEdge;
1582 : :
1583 [ + + ]: 1492 : if (dir == 0) {
1584 [ + + ][ + + ]: 532 : if (j == m_nDiv[1] || k == m_nDiv[2]) return 0.; // return outside
1585 : 476 : iEdge = k*m_nDiv[0]*m_nDiv[1] + j*m_nDiv[0] + i;
1586 : : }
1587 [ + + ]: 960 : else if (dir == 1) {
1588 [ + + ][ + + ]: 428 : if (i == m_nDiv[0] || k == m_nDiv[2]) return 0.; // return outside
1589 : 376 : iEdge = i*m_nDiv[1]*m_nDiv[2] + k*m_nDiv[1] + j;
1590 : : }
1591 [ + - ]: 532 : else if (dir == 2) {
1592 [ + + ][ + + ]: 532 : if (i == m_nDiv[0] || j == m_nDiv[1]) return 0.; // return outside
1593 : 476 : iEdge = j*m_nDiv[0]*m_nDiv[2] + i*m_nDiv[2] + k;
1594 : : }
1595 : 0 : else return -1.;
1596 : :
1597 : 1328 : EdgeStatus edgeStatus = m_vnEdgeStatus[dir][iEdge];
1598 [ + + ]: 1328 : if (edgeStatus == INSIDE) return 1.;
1599 [ - + ]: 80 : else if (edgeStatus == OUTSIDE) return 0.;
1600 : 80 : else return -1.;
1601 : : }
1602 : :
1603 : 0 : bool EBMesher::move_ray(int& nIntersect, double* startPnt, double* endPnt,
1604 : : double tol, int dir, bool bSpecialCase)
1605 : : {
1606 : : //int i, nIteration;
1607 : : int i;
1608 : 0 : bool bMove = true;
1609 : :
1610 [ # # ]: 0 : if (bSpecialCase) m_iOverlap = 0;
1611 : :
1612 [ # # ]: 0 : while (bMove) {
1613 [ # # ]: 0 : if (m_nIteration > 10) {
1614 [ # # ]: 0 : if (bSpecialCase) return true; // special case
1615 [ # # ]: 0 : else if (m_bUseGeom) { // shared/overlapped, faceting case
1616 [ # # ]: 0 : m_mhOverlappedSurf[m_iOverlap] = 0;
1617 [ # # ]: 0 : m_mhOverlappedSurf[m_iOverlap + 1] = 0;
1618 : 0 : return true;
1619 : : }
1620 : : else {
1621 : 0 : return false;
1622 : : }
1623 : : }
1624 : :
1625 [ # # ]: 0 : for (i = 0; i < 3; i++) {
1626 [ # # ]: 0 : if (i != dir) {
1627 : 0 : startPnt[i] += tol;
1628 : 0 : endPnt[i] += tol;
1629 : : }
1630 : : }
1631 : :
1632 [ # # ]: 0 : moab::CartVect ray(endPnt[0] - startPnt[0], endPnt[1] - startPnt[1], endPnt[2] - startPnt[2]);
1633 [ # # ]: 0 : double rayLength = ray.length();
1634 : : moab::ErrorCode rVal;
1635 [ # # ]: 0 : ray.normalize();
1636 : 0 : m_vIntersection.clear();
1637 : 0 : m_vhInterSurf.clear();
1638 : 0 : m_vhInterFacet.clear();
1639 : :
1640 [ # # ]: 0 : std::vector<double> temp_intersects;
1641 [ # # ][ # # ]: 0 : moab::OrientedBoxTreeTool::TrvStats stats;
1642 [ # # ]: 0 : if (m_bUseGeom) {
1643 : : rVal = m_hObbTree->ray_intersect_sets(temp_intersects, m_vhInterSurf,
1644 : : m_vhInterFacet, m_hTreeRoot, tol,
1645 [ # # ]: 0 : startPnt, ray.array(), &rayLength,
1646 [ # # ]: 0 : &stats);
1647 : : }
1648 : : else { // facet input
1649 [ # # ]: 0 : std::vector<moab::EntityHandle> dum_facets_out;
1650 : : rVal = m_hObbTree->ray_intersect_triangles(temp_intersects, dum_facets_out,
1651 : : m_hTreeRoot, tol,
1652 [ # # ]: 0 : startPnt, ray.array(), &rayLength,
1653 [ # # ]: 0 : &stats);
1654 [ # # ]: 0 : m_vhInterSurf.resize(temp_intersects.size());
1655 : : }
1656 : :
1657 [ # # ]: 0 : if (moab::MB_SUCCESS != rVal) {
1658 [ # # ][ # # ]: 0 : std::cerr << "Failed : ray-triangle intersection." << std::endl;
1659 : 0 : return false;
1660 : : }
1661 : :
1662 : 0 : nIntersect = temp_intersects.size();
1663 [ # # ]: 0 : m_vIntersection.resize(nIntersect);
1664 [ # # ]: 0 : for (i = 0; i < nIntersect; i++) {
1665 [ # # ][ # # ]: 0 : IntersectDist temp_inter_dist(temp_intersects[i], i);
1666 [ # # ]: 0 : m_vIntersection[i] = temp_inter_dist;
1667 : : }
1668 [ # # ]: 0 : std::sort(m_vIntersection.begin(), m_vIntersection.end(), less_intersect);
1669 [ # # ]: 0 : bMove = is_ray_move_and_set_overlap_surf(bSpecialCase);
1670 [ # # ]: 0 : m_nIteration++;
1671 : 0 : }
1672 : :
1673 : 0 : std::cout << "ray is moved successfully." << std::endl;
1674 : :
1675 : 0 : return true;
1676 : : }
1677 : :
1678 : 0 : bool EBMesher::is_ray_move_and_set_overlap_surf(bool& bSpecialCase)
1679 : : {
1680 : 0 : int nInter = m_vIntersection.size() - 1;
1681 : :
1682 [ # # ]: 0 : while (m_iOverlap < nInter) {
1683 [ # # ]: 0 : if (m_vIntersection[m_iOverlap + 1].distance - m_vIntersection[m_iOverlap].distance < 1e-7) {
1684 [ # # ]: 0 : if (m_bUseGeom) {
1685 [ # # ][ # # ]: 0 : moab::EntityHandle h1 = m_vhInterSurf[m_vIntersection[m_iOverlap].index];
1686 [ # # ][ # # ]: 0 : moab::EntityHandle h2 = m_vhInterSurf[m_vIntersection[m_iOverlap + 1].index];
1687 : :
1688 [ # # ]: 0 : if (h1 == h2) { // remove too close case
1689 : 0 : bSpecialCase = false;
1690 : 0 : return true;
1691 : : }
1692 [ # # ]: 0 : else if (m_nIteration < 10) { // when ray intesect shared edge by 2 triangles
1693 : 0 : bSpecialCase = true;
1694 : 0 : return true;
1695 : : }
1696 : : else { // overlapped surfaces
1697 [ # # ]: 0 : m_mhOverlappedSurf[h1] = 0;
1698 [ # # ]: 0 : m_mhOverlappedSurf[h2] = 0;
1699 : 0 : m_nIteration = 0;
1700 : 0 : m_iOverlap++;
1701 : : }
1702 : : }
1703 : : else {
1704 [ # # ]: 0 : if (m_nIteration < 10) { // when ray intesect shared edge by 2 triangles
1705 : 0 : bSpecialCase = true;
1706 : 0 : return true;
1707 : : }
1708 : : else { // overlapped surfaces
1709 : : //bSpecialCase = false;
1710 : : //m_nIteration = 0;
1711 : : //m_iOverlap++;
1712 [ # # ][ # # ]: 0 : m_vIntersection.erase(m_vIntersection.begin() + m_iOverlap + 1);
[ # # ]
1713 : 0 : nInter = m_vIntersection.size();
1714 : 0 : m_vhInterSurf.resize(nInter);
1715 : : //m_nIteration = 0;
1716 : : }
1717 : : }
1718 : : }
1719 : 0 : else m_iOverlap++;
1720 : : }
1721 : :
1722 : 0 : return false;
1723 : : }
1724 : :
1725 : 207 : void EBMesher::remove_intersection_duplicates()
1726 : : {
1727 : 207 : int ind = 0;
1728 : 207 : int nInter = m_vIntersection.size() - 1;
1729 : :
1730 [ + + ]: 523 : while (ind < nInter) {
1731 [ + + ]: 316 : if (m_vIntersection[ind + 1].distance - m_vIntersection[ind].distance < 1e-7) {
1732 [ + - ][ + - ]: 109 : moab::EntityHandle h1 = m_vhInterSurf[m_vIntersection[ind].index];
1733 [ + - ][ + - ]: 109 : moab::EntityHandle h2 = m_vhInterSurf[m_vIntersection[ind + 1].index];
1734 : :
1735 [ + + ]: 109 : if (h1 != h2) { // overlapped/shared surfaces
1736 [ + - ][ + - ]: 142 : std::vector<iBase_EntitySetHandle> parents_out1, parents_out2;
1737 [ + - ][ + - ]: 71 : int err = mk_core()->imesh_instance()->getPrnts(reinterpret_cast<iBase_EntitySetHandle> (h1), 1, parents_out1);
[ + - ]
1738 [ + - ]: 71 : IBERRCHK(err, "Failed to get number of surface parents.\n");
1739 [ + - ][ + - ]: 71 : err = mk_core()->imesh_instance()->getPrnts(reinterpret_cast<iBase_EntitySetHandle> (h2), 1, parents_out2);
[ + - ]
1740 [ + - ]: 71 : IBERRCHK(err, "Failed to get number of surface parents.\n");
1741 [ + - ][ + - ]: 71 : if (parents_out1.size() == 1 && parents_out2.size() == 1) {
[ + - ]
1742 [ + - ][ + - ]: 71 : if (parents_out1[0] != parents_out2[0]) { // parent volues are also different
[ - + ]
1743 [ # # ]: 0 : m_mhOverlappedSurf[h1] = 0;
1744 [ # # ]: 0 : m_mhOverlappedSurf[h2] = 0;
1745 : : }
1746 : 71 : }
1747 : : }
1748 [ + - ][ + - ]: 109 : m_vIntersection.erase(m_vIntersection.begin() + ind + 1);
[ + - ]
1749 : 109 : nInter--;
1750 : : }
1751 : 207 : else ind++;
1752 : : }
1753 : 207 : }
1754 : :
1755 : 0 : bool EBMesher::get_volume_fraction(int vol_frac_div)
1756 : : {
1757 [ # # ][ # # ]: 0 : std::cout << "starting get_volume_fraction." << std::endl;
1758 : : int i, j, k, l, n, iHex, dir, nIntersect, rayIndex[3], index[3],
1759 : : otherDir1, otherDir2, err, nParent;
1760 : : (void) otherDir1;
1761 : : (void) otherDir2;
1762 : 0 : double tolerance = 1e-12, dDistance, elem_origin[3],
1763 : : elem_interval_size[3], startPnt[3], endPnt[3], rayMaxEnd[3];
1764 : : moab::ErrorCode rval;
1765 : 0 : bool bOverlapSecond, bClosed = true;
1766 [ # # ]: 0 : std::vector<iBase_EntityHandle> edge_handles;
1767 : :
1768 : 0 : double dTotEdgeElem = 0.;
1769 [ # # ]: 0 : for (i = 0; i < 3; i++) {
1770 : 0 : rayMaxEnd[i] = m_origin_coords[i] + m_nDiv[i]*m_dIntervalSize[i];
1771 : 0 : dTotEdgeElem += m_dIntervalSize[i]*(vol_frac_div + 1)*(vol_frac_div + 1);
1772 : : }
1773 : :
1774 : : // get all hexes
1775 [ # # ]: 0 : std::vector<iBase_EntityHandle> hex_handles;
1776 : : err = mk_core()->imesh_instance()->getEntities(m_hRootSet, iBase_REGION,
1777 [ # # ][ # # ]: 0 : iMesh_HEXAHEDRON, hex_handles);
[ # # ]
1778 [ # # ]: 0 : IBERRCHK(err, "Failed to get hexes.\n");
1779 : :
1780 : 0 : int hex_size = hex_handles.size();
1781 : : /*
1782 : : std::vector<int> hex_status(hex_size);
1783 : : //int* hex_status = NULL;
1784 : : err = mk_core()->imesh_instance()->getIntArrData(&hex_handles[0], hex_size,
1785 : : m_elemStatusTag, &hex_status[0]);
1786 : : */
1787 : : int error;
1788 [ # # ]: 0 : std::vector<int> frac_length;
1789 [ # # ][ # # ]: 0 : int* hex_status = new int[hex_size];
1790 : 0 : int hex_status_alloc = sizeof(int)*hex_size;
1791 : 0 : int hex_status_size = 0;
1792 [ # # ]: 0 : iMesh_getIntArrData(m_mesh, &hex_handles[0],
1793 : : hex_size, m_elemStatusTag,
1794 : : &hex_status, &hex_status_alloc,
1795 [ # # ]: 0 : &hex_status_size, &error);
1796 [ # # ][ # # ]: 0 : ERRORRF("Failed to get hex status.\n");
[ # # ]
1797 : :
1798 [ # # ]: 0 : for (n = 0; n < hex_size; n++) {
1799 [ # # ]: 0 : if (hex_status[n] == 2) { // just boundary
1800 [ # # ]: 0 : std::map<moab::EntityHandle, VolFrac> vol_fraction;
1801 [ # # ]: 0 : std::map<moab::EntityHandle, VolFrac>::iterator vf_iter;
1802 [ # # ]: 0 : std::map<moab::EntityHandle, VolFrac>::iterator vf_end_iter;
1803 [ # # ]: 0 : iHex = moab_instance()->id_from_handle(reinterpret_cast<moab::EntityHandle>
1804 [ # # ][ # # ]: 0 : (hex_handles[n])) - m_iStartHex;
1805 : 0 : index[0] = (iHex%(m_nDiv[0]*m_nDiv[1]))%m_nDiv[0];
1806 : 0 : index[1] = (iHex%(m_nDiv[0]*m_nDiv[1]))/m_nDiv[0];
1807 : 0 : index[2] = iHex/m_nDiv[0]/m_nDiv[1];
1808 : :
1809 [ # # ]: 0 : for (i = 0; i < 3; i++) {
1810 : 0 : elem_origin[i] = m_origin_coords[i] + index[i]*m_dIntervalSize[i];
1811 : 0 : elem_interval_size[i] = m_dIntervalSize[i]/vol_frac_div;
1812 : : }
1813 : :
1814 [ # # ]: 0 : for (dir = 0; dir < 3; dir++) { // x, y, z directions
1815 : 0 : otherDir1 = (dir + 1)%3;
1816 : 0 : otherDir2 = (dir + 2)%3;
1817 : :
1818 [ # # ]: 0 : for (j = 0; j < vol_frac_div + 1; j++) {
1819 [ # # ]: 0 : for (i = 0; i < vol_frac_div + 1; i++) {
1820 : : // get ray start and end points
1821 [ # # ]: 0 : if (dir == 0) { // x coordinate ray
1822 : 0 : rayIndex[0] = 0;
1823 : 0 : rayIndex[1] = i;
1824 : 0 : rayIndex[2] = j;
1825 : : }
1826 [ # # ]: 0 : else if (dir == 1) { // y coordinate ray
1827 : 0 : rayIndex[0] = j;
1828 : 0 : rayIndex[1] = 0;
1829 : 0 : rayIndex[2] = i;
1830 : : }
1831 [ # # ]: 0 : else if (dir == 2) { // z coordinate ray
1832 : 0 : rayIndex[0] = i;
1833 : 0 : rayIndex[1] = j;
1834 : 0 : rayIndex[2] = 0;
1835 : : }
1836 [ # # ]: 0 : else IBERRCHK(iBase_FAILURE, "Ray direction should be from 0 to 2.");
1837 : :
1838 [ # # ]: 0 : for (k = 0; k < 3; k++) {
1839 [ # # ]: 0 : if (k == dir) {
1840 : 0 : startPnt[k] = elem_origin[k];
1841 : 0 : endPnt[k] = startPnt[k] + m_dIntervalSize[k];
1842 : : }
1843 : : else {
1844 : 0 : startPnt[k] = elem_origin[k] + rayIndex[k]*elem_interval_size[k];
1845 : 0 : endPnt[k] = startPnt[k];
1846 : : }
1847 : : }
1848 : :
1849 : : // ray-tracing
1850 [ # # ]: 0 : if (!fire_ray(nIntersect, startPnt, endPnt,
1851 [ # # ]: 0 : tolerance, dir, m_dIntervalSize[dir])) return iBase_FAILURE;
1852 : :
1853 [ # # ]: 0 : if (nIntersect > 0) { // ray is intersected
1854 : 0 : bOverlapSecond = false;
1855 : 0 : bClosed = true;
1856 [ # # ]: 0 : for (k = 0; k < nIntersect; k++) {
1857 [ # # ]: 0 : std::vector<moab::EntityHandle> parents;
1858 : : //moab::Range parents;
1859 [ # # ][ # # ]: 0 : rval = moab_instance()->get_parent_meshsets(m_vhInterSurf[m_vIntersection[k].index],
[ # # ]
1860 [ # # ]: 0 : parents);
1861 [ # # ]: 0 : if (rval != moab::MB_SUCCESS) {
1862 [ # # ][ # # ]: 0 : std::cerr << "Couldn't get parents." << std::endl;
1863 : 0 : return false;
1864 : : }
1865 : :
1866 : 0 : nParent = parents.size();
1867 [ # # ]: 0 : dDistance = m_vIntersection[k].distance;
1868 : :
1869 [ # # ][ # # ]: 0 : if (is_shared_overlapped_surf(k)) {
1870 [ # # ]: 0 : if (nParent == 2) { // shared surface
1871 [ # # ]: 0 : for (l = 0; l < 2; l++) {
1872 [ # # ]: 0 : if (l == 1) {
1873 : 0 : dDistance *= -1.;
1874 : 0 : bClosed = false;
1875 : : }
1876 : 0 : else bClosed = true;
1877 : :
1878 [ # # ][ # # ]: 0 : vf_iter = vol_fraction.find(parents[l]);
1879 [ # # ][ # # ]: 0 : if (vf_iter == vol_fraction.end()) {
1880 [ # # ]: 0 : VolFrac temp_vf(dDistance, bClosed);
1881 [ # # ][ # # ]: 0 : vol_fraction[parents[l]] = temp_vf;
1882 : : //std::cout << "iHex=" << iHex << ", vh="
1883 : : // << parents[l] << ", dir=" << dir << ", dDistance1="
1884 : : // << dDistance << ", vol="
1885 : : // << temp_vf.vol << std::endl;
1886 : : }
1887 : : else {
1888 [ # # ]: 0 : vf_iter->second.vol += dDistance;
1889 [ # # ]: 0 : vf_iter->second.closed = bClosed;
1890 : : //std::cout << "iHex=" << iHex << ", vh="
1891 : : // << vf_iter->first << ", dir=" << dir << ", dDistance1="
1892 : : // << dDistance << ", vol="
1893 : : // << vf_iter->second.vol << std::endl;
1894 : : }
1895 : : }
1896 : : }
1897 [ # # ]: 0 : else if (nParent == 1) { // overlapped surface
1898 [ # # ]: 0 : if (bOverlapSecond) { // second intersection
1899 : : /*
1900 : : for (int m = 0; m < 3; m++) {
1901 : : ePnt[m] = startPnt[m];
1902 : : }
1903 : : ePnt[dir] += dDistance;
1904 : : */
1905 : 0 : dDistance *= -1.;
1906 : 0 : bClosed = false;
1907 : 0 : bOverlapSecond = false;
1908 : : }
1909 : : else {
1910 : : /*
1911 : : // make edge
1912 : : for (int m = 0; m < 3; m++) {
1913 : : if (bClosed) ePnt[m] = startPnt[m];
1914 : : ePnt[m + 3] = ePnt[m];
1915 : : }
1916 : : ePnt[dir + 3] += dDistance;
1917 : : if (!make_edge(ePnt, edge_handles)) return iBase_FAILURE;
1918 : : */
1919 : 0 : bOverlapSecond = true;// first intersection
1920 : 0 : bClosed = true;
1921 : : }
1922 : :
1923 [ # # ][ # # ]: 0 : vf_iter = vol_fraction.find(parents[0]);
1924 [ # # ][ # # ]: 0 : if (vf_iter == vol_fraction.end()) {
1925 [ # # ]: 0 : VolFrac temp_vf(dDistance, bClosed);
1926 [ # # ][ # # ]: 0 : vol_fraction[parents[0]] = temp_vf;
1927 : : //std::cout << "iHex=" << iHex << ", vh="
1928 : : // << parents[0] << ", dDistance2="
1929 : : // << dDistance << ", vol="
1930 : : // << temp_vf.vol << std::endl;
1931 : : }
1932 : : else {
1933 [ # # ]: 0 : vf_iter->second.vol += dDistance;
1934 [ # # ]: 0 : vf_iter->second.closed = bClosed;
1935 : : //std::cout << "iHex=" << iHex << ", vh="
1936 : : // << vf_iter->first << ", dir=" << dir << ", dDistance2="
1937 : : // << dDistance << ", vol="
1938 : : // << vf_iter->second.vol << std::endl;
1939 : : }
1940 : : }
1941 : 0 : else return iBase_FAILURE;
1942 : : }
1943 : : else { // normal case
1944 [ # # ]: 0 : if (nParent != 1) return iBase_FAILURE;
1945 : :
1946 [ # # ][ # # ]: 0 : if (!is_same_direct_to_ray(k, dir)) {
1947 : 0 : dDistance *= -1.; // outside
1948 : 0 : bClosed = false;
1949 : : }
1950 : 0 : else bClosed = true;
1951 : :
1952 [ # # ][ # # ]: 0 : vf_iter = vol_fraction.find(parents[0]);
1953 [ # # ][ # # ]: 0 : if (vf_iter == vol_fraction.end()) {
1954 [ # # ]: 0 : VolFrac temp_vf(dDistance, bClosed);
1955 [ # # ][ # # ]: 0 : vol_fraction[parents[0]] = temp_vf;
1956 : : //std::cout << "iHex=" << iHex << ", vh="
1957 : : // << parents[0] << ", dir=" << dir << ", dDistance3="
1958 : : // << dDistance << ", vol="
1959 : : // << temp_vf.vol << std::endl;
1960 : : }
1961 : : else {
1962 [ # # ]: 0 : vf_iter->second.vol += dDistance;
1963 [ # # ][ # # ]: 0 : vf_iter->second.closed = bClosed;
1964 : : //std::cout << "iHex=" << iHex << ", vh="
1965 : : // << vf_iter->first << ", dir=" << dir << ", dDistance3="
1966 : : // << dDistance << ", vol="
1967 : : // << vf_iter->second.vol << std::endl;
1968 : : }
1969 : : }
1970 : 0 : }
1971 : :
1972 : : // if fraction is not closed, add interval size to close
1973 : 0 : vf_iter = vol_fraction.begin();
1974 : 0 : vf_end_iter = vol_fraction.end();
1975 [ # # ][ # # ]: 0 : for (; vf_iter != vf_end_iter; vf_iter++) {
[ # # ]
1976 [ # # ][ # # ]: 0 : if (!vf_iter->second.closed) {
1977 [ # # ]: 0 : vf_iter->second.vol += m_dIntervalSize[dir];
1978 [ # # ]: 0 : vf_iter->second.closed = true;
1979 : : /*
1980 : : for (int m = 0; m < 3; m++) {
1981 : : ePnt[m + 3] = startPnt[m];
1982 : : }
1983 : : ePnt[dir + 3] += m_dIntervalSize[dir];
1984 : : if (!make_edge(ePnt, edge_handles)) return iBase_FAILURE;
1985 : : */
1986 : : //std::cout << "iHex=" << iHex << ", vh="
1987 : : // << vf_iter->first << ", dir=" << dir << ", dDistance4="
1988 : : // << m_dIntervalSize[dir] << ", vol="
1989 : : // << vf_iter->second.vol << std::endl;
1990 : : }
1991 : : }
1992 : : }
1993 : : else { // decide if it is fully inside
1994 [ # # ]: 0 : if (!fire_ray(nIntersect, startPnt, rayMaxEnd,
1995 [ # # ]: 0 : tolerance, dir, m_nDiv[dir]*m_dIntervalSize[dir])) return iBase_FAILURE;
1996 : :
1997 [ # # ]: 0 : if (nIntersect > 0) { // fully inside
1998 [ # # ][ # # ]: 0 : if (is_shared_overlapped_surf(0) || // shared or overlapped
[ # # ][ # # ]
1999 [ # # ]: 0 : is_same_direct_to_ray(0, dir)) { // other inside case
2000 [ # # ]: 0 : moab::Range parents;
2001 [ # # ][ # # ]: 0 : rval = moab_instance()->get_parent_meshsets(m_vhInterSurf[m_vIntersection[0].index],
[ # # ]
2002 [ # # ]: 0 : parents);
2003 [ # # ]: 0 : if (rval != moab::MB_SUCCESS) {
2004 [ # # ][ # # ]: 0 : std::cerr << "Couldn't get parents." << std::endl;
2005 : 0 : return false;
2006 : : }
2007 : :
2008 [ # # ]: 0 : moab::EntityHandle parent_entity = parents.pop_front();
2009 [ # # ]: 0 : vf_iter = vol_fraction.find(parent_entity);
2010 [ # # ][ # # ]: 0 : if (vf_iter == vf_end_iter) {
2011 [ # # ]: 0 : VolFrac temp_vf(m_dIntervalSize[dir], true);
2012 [ # # ]: 0 : vol_fraction[parent_entity] = temp_vf;
2013 : : //std::cout << "iHex=" << iHex << ", vh="
2014 : : // << parents[0] << ", dir=" << dir << ", dDistance5="
2015 : : // << dDistance << ", vol="
2016 : : // << temp_vf.vol << std::endl;
2017 : : }
2018 : : else {
2019 [ # # ]: 0 : vf_iter->second.vol += m_dIntervalSize[dir];
2020 [ # # ][ # # ]: 0 : vf_iter->second.closed = bClosed;
2021 : : //std::cout << "iHex=" << iHex << ", vh="
2022 : : // << vf_iter->first << ", dir=" << dir << ", dDistance5="
2023 : : // << dDistance << ", vol="
2024 : : // << vf_iter->second.vol << std::endl;
2025 : 0 : }
2026 : : }
2027 : : }
2028 : : }
2029 : : }
2030 : : }
2031 : : }
2032 : :
2033 : 0 : int vol_frac_length = vol_fraction.size();
2034 [ # # ][ # # ]: 0 : std::vector<int> vol_frac_id(vol_frac_length);
2035 [ # # ]: 0 : std::vector<double> vol_frac(vol_frac_length);
2036 : 0 : int vol_frac_id_size = vol_frac_length*sizeof(int);
2037 : 0 : int vol_frac_size = vol_frac_length*sizeof(double);
2038 : 0 : vf_iter = vol_fraction.begin();
2039 : 0 : vf_end_iter = vol_fraction.end();
2040 [ # # ][ # # ]: 0 : for (int m = 0; vf_iter != vf_end_iter; vf_iter++, m++) {
[ # # ]
2041 [ # # ][ # # ]: 0 : vol_frac_id[m] = vf_iter->first;
2042 [ # # ][ # # ]: 0 : vol_frac[m] = vf_iter->second.vol/dTotEdgeElem;
2043 : : //std::cout << "iHex=" << iHex << ", i=" << index[0]
2044 : : // << ", j=" << index[1] << ", k=" << index[2]
2045 : : // << ", vol_handle=" << vf_iter->first
2046 : : // << ", vol=" << vf_iter->second.vol
2047 : : // << ", vol_fraction=" << vf_iter->second.vol/dTotEdgeElem
2048 : : // << std::endl;
2049 : : }
2050 [ # # ]: 0 : const void* vol_frac_ids = &vol_frac_id[0];
2051 [ # # ]: 0 : const void* vol_fracs = &vol_frac[0];
2052 : :
2053 : : // set volume fraction information as tag
2054 [ # # ]: 0 : rval = moab_instance()->tag_set_data(reinterpret_cast<moab::Tag> (m_volFracLengthTag),
2055 [ # # ]: 0 : reinterpret_cast<moab::EntityHandle*> (&hex_handles[n]),
2056 [ # # ]: 0 : 1, &vol_frac_length);
2057 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
2058 : :
2059 [ # # ]: 0 : rval = moab_instance()->tag_set_by_ptr(reinterpret_cast<moab::Tag> (m_volFracHandleTag),
2060 [ # # ]: 0 : reinterpret_cast<moab::EntityHandle*> (&hex_handles[n]),
2061 [ # # ]: 0 : 1, &vol_frac_ids, &vol_frac_id_size);
2062 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
2063 : :
2064 [ # # ]: 0 : rval = moab_instance()->tag_set_by_ptr(reinterpret_cast<moab::Tag> (m_volFracTag),
2065 [ # # ]: 0 : reinterpret_cast<moab::EntityHandle*> (&hex_handles[n]),
2066 [ # # ]: 0 : 1, &vol_fracs, &vol_frac_size);
2067 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
2068 : :
2069 : : if (debug_ebmesh) {
2070 : : for (int v = 0; v < vol_frac_length; v++) {
2071 : : std::cout << "#_boundary_hex_handles=" << hex_handles[n]
2072 : : << ",vol_frac_handle[" << v << "]=" << vol_frac_id[v]
2073 : : << ",vol_fracs[" << v << "]=" << vol_frac[v]
2074 : : << std::endl;
2075 : : }
2076 : 0 : }
2077 : : }
2078 : : }
2079 : :
2080 [ # # ][ # # ]: 0 : std::cout << "end get_volume_fraction." << std::endl;
2081 [ # # ]: 0 : delete [] hex_status;
2082 : :
2083 : 0 : return true;
2084 : : }
2085 : :
2086 : 0 : bool EBMesher::is_same_direct_to_ray(int i, int dir)
2087 : : {
2088 [ # # ][ # # ]: 0 : moab::CartVect coords[3], normal(0.0), ray_dir(rayDir[dir]);
[ # # ][ # # ]
2089 : : const moab::EntityHandle *conn;
2090 : : int len;
2091 : :
2092 : : // get triangle normal vector
2093 [ # # ][ # # ]: 0 : moab::ErrorCode rval = moab_instance()->get_connectivity(m_vhInterFacet[m_vIntersection[i].index],
[ # # ]
2094 [ # # ]: 0 : conn, len);
2095 [ # # ][ # # ]: 0 : assert(moab::MB_SUCCESS == rval && 3 == len);
2096 : :
2097 [ # # ][ # # ]: 0 : rval = moab_instance()->get_coords(conn, 3, coords[0].array());
[ # # ]
2098 : : (void) rval;
2099 [ # # ]: 0 : assert(moab::MB_SUCCESS == rval);
2100 : :
2101 [ # # ]: 0 : coords[1] -= coords[0];
2102 [ # # ]: 0 : coords[2] -= coords[0];
2103 [ # # ][ # # ]: 0 : normal += coords[1] * coords[2];
2104 [ # # ]: 0 : normal.normalize();
2105 : :
2106 [ # # ]: 0 : return angle(normal, ray_dir) < .5*PI;
2107 : : }
2108 : :
2109 : 0 : void EBMesher::set_obb_tree_box_dimension()
2110 : : {
2111 [ # # ][ # # ]: 0 : std::cout << "starting to construct obb tree." << std::endl;
2112 [ # # ]: 0 : moab::ErrorCode rval = m_GeomTopoTool->construct_obb_trees(m_bUseWholeGeom);
2113 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
2114 : :
2115 [ # # ]: 0 : m_hObbTree = m_GeomTopoTool->obb_tree();
2116 [ # # ]: 0 : m_hTreeRoot = m_GeomTopoTool->get_one_vol_root();
2117 : :
2118 [ # # ][ # # ]: 0 : moab::CartVect box_center, box_axis1, box_axis2, box_axis3;
[ # # ][ # # ]
2119 [ # # ]: 0 : for (int i = 0; i < 3; i++) {
2120 : 0 : m_minCoord[i] = HUGE_VAL;
2121 : 0 : m_maxCoord[i] = 0.;
2122 : : }
2123 : : rval = m_hObbTree->box(m_hTreeRoot, box_center.array(),
2124 : : box_axis1.array(), box_axis2.array(),
2125 [ # # ][ # # ]: 0 : box_axis3.array());
[ # # ][ # # ]
[ # # ]
2126 [ # # ][ # # ]: 0 : MBERRCHK(rval, mk_core()->moab_instance());
[ # # ][ # # ]
[ # # ][ # # ]
2127 : :
2128 : : // cartesian box corners
2129 [ # # ][ # # ]: 0 : moab::CartVect corners[8] = {box_center + box_axis1 + box_axis2 + box_axis3,
2130 [ # # ][ # # ]: 0 : box_center + box_axis1 + box_axis2 - box_axis3,
2131 [ # # ][ # # ]: 0 : box_center + box_axis1 - box_axis2 + box_axis3,
2132 [ # # ][ # # ]: 0 : box_center + box_axis1 - box_axis2 - box_axis3,
2133 [ # # ][ # # ]: 0 : box_center - box_axis1 + box_axis2 + box_axis3,
2134 [ # # ][ # # ]: 0 : box_center - box_axis1 + box_axis2 - box_axis3,
2135 [ # # ][ # # ]: 0 : box_center - box_axis1 - box_axis2 + box_axis3,
2136 [ # # ][ # # ]: 0 : box_center - box_axis1 - box_axis2 - box_axis3};
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
2137 : :
2138 : : // get the max, min cartesian box corners
2139 [ # # ]: 0 : for (int i = 0; i < 8; i++) {
2140 [ # # ]: 0 : for (int j = 0; j < 3; j++) {
2141 [ # # ][ # # ]: 0 : if (corners[i][j] < m_minCoord[j]) m_minCoord[j] = corners[i][j];
[ # # ]
2142 [ # # ][ # # ]: 0 : if (corners[i][j] > m_maxCoord[j]) m_maxCoord[j] = corners[i][j];
[ # # ]
2143 : : }
2144 : : }
2145 [ # # ][ # # ]: 0 : std::cout << "finished to construct obb tree." << std::endl;
2146 : 0 : }
2147 [ + - ][ + - ]: 156 : } // namespace MeshKit
|
