Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
BoundBox.cpp
Go to the documentation of this file.
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
6 
7 namespace moab
8 {
9 ErrorCode BoundBox::update( Interface& iface, const Range& elems, bool spherical, double radius )
10 {
11  ErrorCode rval;
12  bMin = CartVect( HUGE_VAL );
13  bMax = CartVect( -HUGE_VAL );
14 
15  CartVect coords;
16  EntityHandle const *conn = NULL, *conn2 = NULL;
17  int len = 0, len2 = 0;
19  CartVect coordverts[27]; // maximum nodes per element supported by MOAB
20 
21  // vertices
22  const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
23  for( i = elems.begin(); i != elem_begin; ++i )
24  {
25  rval = iface.get_coords( &*i, 1, coords.array() );
26  if( MB_SUCCESS != rval ) return rval;
27  update_min( coords.array() );
28  update_max( coords.array() );
29  }
30 
31  // elements with vertex-handle connectivity list
32  const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
33  std::vector< EntityHandle > dum_vector;
34  for( i = elem_begin; i != poly_begin; ++i )
35  {
36  rval = iface.get_connectivity( *i, conn, len, true, &dum_vector );
37  if( MB_SUCCESS != rval ) return rval;
38 
39  rval = iface.get_coords( conn, len, &( coordverts[0][0] ) );
40  if( MB_SUCCESS != rval ) return rval;
41  for( int j = 0; j < len; ++j )
42  {
43  update_min( coordverts[j].array() );
44  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  if( spherical ) update_box_spherical_elem( coordverts, len, radius );
49  }
50 
51  // polyhedra
52  const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
53  for( i = poly_begin; i != set_begin; ++i )
54  {
55  rval = iface.get_connectivity( *i, conn, len, true );
56  if( MB_SUCCESS != rval ) return rval;
57 
58  for( int j = 0; j < len; ++j )
59  {
60  rval = iface.get_connectivity( conn[j], conn2, len2 );
61  if( MB_SUCCESS != rval ) return rval;
62  rval = iface.get_coords( conn2, len2, &( coordverts[0][0] ) );
63  if( MB_SUCCESS != rval ) return rval;
64  for( int k = 0; k < len2; ++k )
65  {
66  update_min( coordverts[k].array() );
67  update_max( coordverts[k].array() );
68  }
69  }
70  }
71 
72  // sets
73  BoundBox box;
74  for( i = set_begin; i != elems.end(); ++i )
75  {
76  Range tmp_elems;
77  rval = iface.get_entities_by_handle( *i, tmp_elems );
78  if( MB_SUCCESS != rval ) return rval;
79  rval = box.update( iface, tmp_elems, spherical, radius );
80  if( MB_SUCCESS != rval ) return rval;
81 
82  update( box );
83  }
84 
85  return MB_SUCCESS;
86 }
87 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  int plane = -1;
91  CartVect pos;
92  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  for( int i = 0; i < len && i < 10; i++ )
95  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  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  for( int i = 0; i < len; i++ )
102  {
103  int i1 = ( i + 1 ) % len; // next vertex in polygon
104  // check intx with x axis
105  double ax = in_plane_positions[2 * i], ay = in_plane_positions[2 * i + 1];
106  double bx = in_plane_positions[2 * i1], by = in_plane_positions[2 * i1 + 1];
107  if( ay * by < 0 ) // it intersects with x axis
108  {
109  double alfa = ay / ( ay - by );
110  double xintx = ax + alfa * ( bx - ax );
111  IntxUtils::reverse_gnomonic_projection( xintx, 0, R, plane, pos );
112  update_min( pos.array() );
113  update_max( pos.array() );
114  }
115  if( ax * bx < 0 ) // it intersects with y axis
116  {
117  double alfa = ax / ( ax - bx );
118  double yintx = ay + alfa * ( by - ay );
119  IntxUtils::reverse_gnomonic_projection( 0, yintx, R, plane, pos );
120  update_min( pos.array() );
121  update_max( pos.array() );
122  }
123  oriented_area2 += ( ax * by - ay * bx );
124  }
125  if( fabs( oriented_area2 ) > R * R * 1.e-6 ) // origin is in the interior, add the center
126  {
127  IntxUtils::reverse_gnomonic_projection( 0, 0, R, plane, pos );
128  update_min( pos.array() );
129  update_max( pos.array() );
130  }
131 }
132 } // namespace moab