Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
SweptElementData.hpp
Go to the documentation of this file.
1 /**
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 #ifndef SWEPT_ELEMENT_DATA_HPP
17 #define SWEPT_ELEMENT_DATA_HPP
18 
19 //
20 // Class: SweptElementData
21 //
22 // Purpose: represent a rectangular element of mesh
23 //
24 // A SweptElement represents a rectangular element of mesh, including both vertices and
25 // elements, and the parametric space used to address that element. Vertex data,
26 // i.e. coordinates, may not be stored directly in the element, but the element returns
27 // information about the vertex handles of vertices in the element. Vertex and element
28 // handles associated with the element are each contiguous.
29 
30 #include "SequenceData.hpp"
31 #include "moab/HomXform.hpp"
32 #include "moab/CN.hpp"
33 #include "SweptVertexData.hpp"
34 #include "Internals.hpp"
35 #include "moab/Range.hpp"
36 
37 #include <vector>
38 #include <algorithm>
39 
40 namespace moab
41 {
42 
44 {
45 
46  //! structure to hold references to bounding vertex blocks
48  {
49  private:
53 
54  public:
55  friend class SweptElementData;
56 
57  VertexDataRef( const HomCoord& min, const HomCoord& max, const HomXform& tmp_xform, SweptVertexData* this_seq );
58 
59  bool contains( const HomCoord& coords ) const;
60  };
61 
62  private:
63  //! parameter min/max/stride, in homogeneous coords ijkh
65 
66  //! difference between max and min params plus one (i.e. # VERTICES in
67  //! each parametric direction)
68  int dIJK[3];
69 
70  //! difference between max and min params (i.e. # ELEMENTS in
71  //! each parametric direction)
72  int dIJKm1[3];
73 
74  //! bare constructor, so compiler doesn't create one for me
76 
77  //! list of bounding vertex blocks
78  std::vector< VertexDataRef > vertexSeqRefs;
79 
80  public:
81  //! constructor
83  const int imin,
84  const int jmin,
85  const int kmin,
86  const int imax,
87  const int jmax,
88  const int kmax,
89  const int* Cq );
90 
91  virtual ~SweptElementData();
92 
93  //! get handle of vertex at homogeneous coords
94  inline EntityHandle get_vertex( const HomCoord& coords ) const;
95  inline EntityHandle get_vertex( int i, int j, int k ) const
96  {
97  return get_vertex( HomCoord( i, j, k ) );
98  }
99 
100  //! get handle of element at i, j, k
101  EntityHandle get_element( const int i, const int j, const int k ) const;
102 
103  //! get min params for this element
104  const HomCoord& min_params() const;
105 
106  //! get max params for this element
107  const HomCoord& max_params() const;
108 
109  //! get the number of vertices in each direction, inclusive
110  void param_extents( int& di, int& dj, int& dk ) const;
111 
112  //! given a handle, get the corresponding parameters
113  ErrorCode get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const;
114 
115  //! convenience functions for parameter extents
116  int i_min() const
117  {
118  return ( elementParams[0].hom_coord() )[0];
119  }
120  int j_min() const
121  {
122  return ( elementParams[0].hom_coord() )[1];
123  }
124  int k_min() const
125  {
126  return ( elementParams[0].hom_coord() )[2];
127  }
128  int i_max() const
129  {
130  return ( elementParams[1].hom_coord() )[0];
131  }
132  int j_max() const
133  {
134  return ( elementParams[1].hom_coord() )[1];
135  }
136  int k_max() const
137  {
138  return ( elementParams[1].hom_coord() )[2];
139  }
140 
141  //! test the bounding vertex sequences and determine whether they fully
142  //! define the vertices covering this element block's parameter space
143  bool boundary_complete() const;
144 
145  //! test whether this sequence contains these parameters
146  inline bool contains( const HomCoord& coords ) const;
147 
148  //! get connectivity of an entity given entity's parameters
149  inline ErrorCode get_params_connectivity( const int i,
150  const int j,
151  const int k,
152  std::vector< EntityHandle >& connectivity ) const;
153 
154  //! add a vertex seq ref to this element sequence;
155  //! if bb_input is true, bounding box (in eseq-local coords) of vseq being added
156  //! is input in bb_min and bb_max (allows partial sharing of vseq rather than the whole
157  //! vseq); if it's false, the whole vseq is referenced and the eseq-local coordinates
158  //! is computed from the transformed bounding box of the vseq
160  const HomCoord& p1,
161  const HomCoord& q1,
162  const HomCoord& p2,
163  const HomCoord& q2,
164  const HomCoord& p3,
165  const HomCoord& q3,
166  bool bb_input = false,
167  const HomCoord& bb_min = HomCoord::unitv[0],
168  const HomCoord& bb_max = HomCoord::unitv[0] );
169 
171  EntityHandle end,
172  const int* sequence_data_sizes,
173  const int* tag_data_sizes ) const;
174 
175  static EntityID calc_num_entities( EntityHandle start_handle, int irange, int jrange, int krange );
176 
177  unsigned long get_memory_use() const;
178 };
179 
180 inline EntityHandle SweptElementData::get_element( const int i, const int j, const int k ) const
181 {
182  return start_handle() + ( i - i_min() ) + ( j - j_min() ) * dIJKm1[0] + ( k - k_min() ) * dIJKm1[0] * dIJKm1[1];
183 }
184 
186 {
187  return elementParams[0];
188 }
189 
191 {
192  return elementParams[1];
193 }
194 
195 //! get the number of vertices in each direction, inclusive
196 inline void SweptElementData::param_extents( int& di, int& dj, int& dk ) const
197 {
198  di = dIJK[0];
199  dj = dIJK[1];
200  dk = dIJK[2];
201 }
202 
203 inline ErrorCode SweptElementData::get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const
204 {
205  if( TYPE_FROM_HANDLE( ehandle ) != TYPE_FROM_HANDLE( start_handle() ) ) return MB_FAILURE;
206 
207  int hdiff = ehandle - start_handle();
208 
209  // use double ?: test below because on some platforms, both sides of the : are
210  // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero
211  k = ( dIJKm1[1] > 0 ? hdiff / ( dIJKm1[1] > 0 ? dIJKm1[0] * dIJKm1[1] : 1 ) : 0 );
212  j = ( hdiff - ( k * dIJKm1[0] * dIJKm1[1] ) ) / dIJKm1[0];
213  i = hdiff % dIJKm1[0];
214 
215  k += elementParams[0].k();
216  j += elementParams[0].j();
217  i += elementParams[0].i();
218 
219  return ( ehandle >= start_handle() && ehandle < start_handle() + size() && i >= i_min() && i <= i_max() &&
220  j >= j_min() && j <= j_max() && k >= k_min() && k <= k_max() )
221  ? MB_SUCCESS
222  : MB_FAILURE;
223 }
224 
225 inline bool SweptElementData::contains( const HomCoord& temp ) const
226 {
227  // upper bound is < instead of <= because element params max is one less
228  // than vertex params max
229  return ( temp >= elementParams[0] && temp < elementParams[1] );
230 }
231 
232 inline bool SweptElementData::VertexDataRef::contains( const HomCoord& coords ) const
233 {
234  return ( minmax[0] <= coords && minmax[1] >= coords );
235 }
236 
238  const HomCoord& this_max,
239  const HomXform& tmp_xform,
240  SweptVertexData* this_seq )
241  : xform( tmp_xform ), invXform( tmp_xform.inverse() ), srcSeq( this_seq )
242 {
243  minmax[0] = HomCoord( this_min );
244  minmax[1] = HomCoord( this_max );
245 }
246 
248 {
249  assert( boundary_complete() );
250  for( std::vector< VertexDataRef >::const_iterator it = vertexSeqRefs.begin(); it != vertexSeqRefs.end(); ++it )
251  {
252  if( ( *it ).minmax[0] <= coords && ( *it ).minmax[1] >= coords )
253  {
254  // first get the vertex block-local parameters
255  HomCoord local_coords = coords / ( *it ).xform;
256 
257  // now get the vertex handle for those coords
258  return ( *it ).srcSeq->get_vertex( local_coords );
259  }
260  }
261 
262  // got here, it's an error
263  return 0;
264 }
265 
267  const HomCoord& p1,
268  const HomCoord& q1,
269  const HomCoord& p2,
270  const HomCoord& q2,
271  const HomCoord& p3,
272  const HomCoord& q3,
273  bool bb_input,
274  const HomCoord& bb_min,
275  const HomCoord& bb_max )
276 {
277  // compute the transform given the vseq-local parameters and the mapping to
278  // this element sequence's parameters passed in minmax
279  HomXform M;
280  M.three_pt_xform( p1, q1, p2, q2, p3, q3 );
281 
282  // min and max in element seq's parameter system may not be same as those in
283  // vseq's system, so need to take min/max
284 
285  HomCoord minmax[2];
286  if( bb_input )
287  {
288  minmax[0] = bb_min;
289  minmax[1] = bb_max;
290  }
291  else
292  {
293  minmax[0] = vseq->min_params() * M;
294  minmax[1] = vseq->max_params() * M;
295  }
296 
297  // check against other vseq's to make sure they don't overlap
298  for( std::vector< VertexDataRef >::const_iterator vsit = vertexSeqRefs.begin(); vsit != vertexSeqRefs.end();
299  ++vsit )
300  if( ( *vsit ).contains( minmax[0] ) || ( *vsit ).contains( minmax[1] ) ) return MB_FAILURE;
301 
302  HomCoord tmp_min( std::min( minmax[0].i(), minmax[1].i() ), std::min( minmax[0].j(), minmax[1].j() ),
303  std::min( minmax[0].k(), minmax[1].k() ) );
304  HomCoord tmp_max( std::max( minmax[0].i(), minmax[1].i() ), std::max( minmax[0].j(), minmax[1].j() ),
305  std::max( minmax[0].k(), minmax[1].k() ) );
306 
307  // set up a new vertex sequence reference
308  VertexDataRef tmp_seq_ref( tmp_min, tmp_max, M, vseq );
309 
310  // add to the list
311  vertexSeqRefs.push_back( tmp_seq_ref );
312 
313  return MB_SUCCESS;
314 }
315 
317  const int j,
318  const int k,
319  std::vector< EntityHandle >& connectivity ) const
320 {
321  if( contains( HomCoord( i, j, k ) ) == false ) return MB_FAILURE;
322 
323  connectivity.push_back( get_vertex( i, j, k ) );
324  connectivity.push_back( get_vertex( i + 1, j, k ) );
325  if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 2 ) return MB_SUCCESS;
326  connectivity.push_back( get_vertex( i + 1, j + 1, k ) );
327  connectivity.push_back( get_vertex( i, j + 1, k ) );
328  if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 3 ) return MB_SUCCESS;
329  connectivity.push_back( get_vertex( i, j, k + 1 ) );
330  connectivity.push_back( get_vertex( i + 1, j, k + 1 ) );
331  connectivity.push_back( get_vertex( i + 1, j + 1, k + 1 ) );
332  connectivity.push_back( get_vertex( i, j + 1, k + 1 ) );
333  return MB_SUCCESS;
334 }
335 
336 } // namespace moab
337 
338 #endif