Branch data Line data Source code
1 : : #include "moab/Interface.hpp"
2 : : #include "moab/BoundBox.hpp"
3 : :
4 : : // this is used for spherical elements, bounding box needs to be updated due to curvature
5 : : #include "moab/IntxMesh/IntxUtils.hpp"
6 : :
7 : : namespace moab
8 : : {
9 : 22 : ErrorCode BoundBox::update( Interface& iface, const Range& elems, bool spherical, double radius )
10 : : {
11 : : ErrorCode rval;
12 [ + - ]: 22 : bMin = CartVect( HUGE_VAL );
13 [ + - ]: 22 : bMax = CartVect( -HUGE_VAL );
14 : :
15 [ + - ]: 22 : CartVect coords;
16 : 22 : EntityHandle const *conn = NULL, *conn2 = NULL;
17 : 22 : int len = 0, len2 = 0;
18 [ + - ]: 22 : Range::const_iterator i;
19 [ + - ][ + + ]: 616 : CartVect coordverts[27]; // maximum nodes per element supported by MOAB
20 : :
21 : : // vertices
22 [ + - ]: 22 : const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
23 [ + - ][ + - ]: 887 : for( i = elems.begin(); i != elem_begin; ++i )
[ + - ][ + + ]
24 : : {
25 [ + - ][ + - ]: 865 : rval = iface.get_coords( &*i, 1, coords.array() );
[ + - ]
26 [ - + ]: 865 : if( MB_SUCCESS != rval ) return rval;
27 [ + - ][ + - ]: 865 : update_min( coords.array() );
28 [ + - ][ + - ]: 865 : update_max( coords.array() );
29 : : }
30 : :
31 : : // elements with vertex-handle connectivity list
32 [ + - ]: 22 : const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
33 [ + - ]: 22 : std::vector< EntityHandle > dum_vector;
34 [ + - ][ + - ]: 5595 : for( i = elem_begin; i != poly_begin; ++i )
[ + + ]
35 : : {
36 [ + - ][ + - ]: 5573 : rval = iface.get_connectivity( *i, conn, len, true, &dum_vector );
37 [ - + ]: 5573 : if( MB_SUCCESS != rval ) return rval;
38 : :
39 [ + - ][ + - ]: 5573 : rval = iface.get_coords( conn, len, &( coordverts[0][0] ) );
40 [ - + ]: 5573 : if( MB_SUCCESS != rval ) return rval;
41 [ + + ]: 35761 : for( int j = 0; j < len; ++j )
42 : : {
43 [ + - ][ + - ]: 30188 : update_min( coordverts[j].array() );
44 [ + - ][ + - ]: 30188 : update_max( coordverts[j].array() );
45 : : }
46 : : // if spherical, use gnomonic projection to improve the computed box
47 : : // it is used now for coupling only
48 [ - + ][ # # ]: 5573 : if( spherical ) update_box_spherical_elem( coordverts, len, radius );
49 : : }
50 : :
51 : : // polyhedra
52 [ + - ]: 22 : const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
53 [ # # ][ + - ]: 22 : for( i = poly_begin; i != set_begin; ++i )
[ - + ]
54 : : {
55 [ # # ][ # # ]: 0 : rval = iface.get_connectivity( *i, conn, len, true );
56 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
57 : :
58 [ # # ]: 0 : for( int j = 0; j < len; ++j )
59 : : {
60 [ # # ]: 0 : rval = iface.get_connectivity( conn[j], conn2, len2 );
61 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
62 [ # # ][ # # ]: 0 : rval = iface.get_coords( conn2, len2, &( coordverts[0][0] ) );
63 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
64 [ # # ]: 0 : for( int k = 0; k < len2; ++k )
65 : : {
66 [ # # ][ # # ]: 0 : update_min( coordverts[k].array() );
67 [ # # ][ # # ]: 0 : update_max( coordverts[k].array() );
68 : : }
69 : : }
70 : : }
71 : :
72 : : // sets
73 [ + - ]: 44 : BoundBox box;
74 [ # # ][ + - ]: 22 : for( i = set_begin; i != elems.end(); ++i )
[ + - ][ - + ]
75 : : {
76 [ # # ]: 0 : Range tmp_elems;
77 [ # # ][ # # ]: 0 : rval = iface.get_entities_by_handle( *i, tmp_elems );
78 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
79 [ # # ]: 0 : rval = box.update( iface, tmp_elems, spherical, radius );
80 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
81 : :
82 [ # # ][ # # ]: 0 : update( box );
83 : 0 : }
84 : :
85 : 44 : return MB_SUCCESS;
86 : : }
87 : 0 : void BoundBox::update_box_spherical_elem( const CartVect* verts, int len, double R )
88 : : {
89 : : // decide first the gnomonic plane, based on first coordinate
90 : 0 : int plane = -1;
91 [ # # ]: 0 : CartVect pos;
92 [ # # ]: 0 : IntxUtils::decide_gnomonic_plane( verts[0], plane );
93 : : double in_plane_positions[20]; // at most polygons with 10 edges; do we need to revise this?
94 [ # # ][ # # ]: 0 : for( int i = 0; i < len && i < 10; i++ )
95 [ # # ]: 0 : IntxUtils::gnomonic_projection( verts[i], R, plane, in_plane_positions[2 * i], in_plane_positions[2 * i + 1] );
96 : : // look for points on the intersection between edges in gnomonic plane and coordinate axes
97 : : // if yes, reverse projection of intx point, and update box
98 : 0 : double oriented_area2 = 0;
99 : : // if oriented area is != zero, it means the origin is inside the polygon, so we need to upgrade
100 : : // the box with the origin, so we do not miss anything (the top of the dome)
101 [ # # ]: 0 : for( int i = 0; i < len; i++ )
102 : : {
103 : 0 : int i1 = ( i + 1 ) % len; // next vertex in polygon
104 : : // check intx with x axis
105 : 0 : double ax = in_plane_positions[2 * i], ay = in_plane_positions[2 * i + 1];
106 : 0 : double bx = in_plane_positions[2 * i1], by = in_plane_positions[2 * i1 + 1];
107 [ # # ]: 0 : if( ay * by < 0 ) // it intersects with x axis
108 : : {
109 : 0 : double alfa = ay / ( ay - by );
110 : 0 : double xintx = ax + alfa * ( bx - ax );
111 [ # # ]: 0 : IntxUtils::reverse_gnomonic_projection( xintx, 0, R, plane, pos );
112 [ # # ][ # # ]: 0 : update_min( pos.array() );
113 [ # # ][ # # ]: 0 : update_max( pos.array() );
114 : : }
115 [ # # ]: 0 : if( ax * bx < 0 ) // it intersects with y axis
116 : : {
117 : 0 : double alfa = ax / ( ax - bx );
118 : 0 : double yintx = ay + alfa * ( by - ay );
119 [ # # ]: 0 : IntxUtils::reverse_gnomonic_projection( 0, yintx, R, plane, pos );
120 [ # # ][ # # ]: 0 : update_min( pos.array() );
121 [ # # ][ # # ]: 0 : update_max( pos.array() );
122 : : }
123 : 0 : oriented_area2 += ( ax * by - ay * bx );
124 : : }
125 [ # # ]: 0 : if( fabs( oriented_area2 ) > R * R * 1.e-6 ) // origin is in the interior, add the center
126 : : {
127 [ # # ]: 0 : IntxUtils::reverse_gnomonic_projection( 0, 0, R, plane, pos );
128 [ # # ][ # # ]: 0 : update_min( pos.array() );
129 [ # # ][ # # ]: 0 : update_max( pos.array() );
130 : : }
131 : 0 : }
132 [ + - ][ + - ]: 228 : } // namespace moab
|
