Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ScdElementData.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 SCD_ELEMENT_DATA_HPP
17 #define SCD_ELEMENT_DATA_HPP
18 
19 //
20 // Class: ScdElementData
21 //
22 // Purpose: represent a rectangular element of mesh
23 //
24 // A ScdElement 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 "ScdVertexData.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 ScdElementData;
56 
57  VertexDataRef( const HomCoord& min, const HomCoord& max, const HomXform& tmp_xform, ScdVertexData* this_seq );
58 
59  bool contains( const HomCoord& coords ) const;
60  };
61 
62  private:
63  //! parameter min/max/stride for vertices, in homogeneous coords ijkh; elements
64  //! are one less than max
66 
67  //! difference between max and min params plus one (i.e. # VERTICES in
68  //! each parametric direction)
69  int dIJK[3];
70 
71  //! difference between max and min params (i.e. # ELEMENTS in
72  //! each parametric direction)
73  int dIJKm1[3];
74 
75  //! whether scd element rectangle is periodic in i and possibly j
76  int isPeriodic[2];
77 
78  //! bare constructor, so compiler doesn't create one for me
80 
81  //! list of bounding vertex blocks
82  std::vector< VertexDataRef > vertexSeqRefs;
83 
84  public:
85  //! constructor
87  const int imin,
88  const int jmin,
89  const int kmin,
90  const int imax,
91  const int jmax,
92  const int kmax,
93  int* is_periodic );
94 
95  virtual ~ScdElementData();
96 
97  //! get handle of vertex at homogeneous coords
98  inline EntityHandle get_vertex( const HomCoord& coords ) const;
99  inline EntityHandle get_vertex( int i, int j, int k ) const
100  {
101  return get_vertex( HomCoord( i, j, k ) );
102  }
103 
104  //! get handle of element at i, j, k
105  EntityHandle get_element( const int i, const int j, const int k ) const;
106 
107  //! get min params for this element
108  const HomCoord& min_params() const;
109 
110  //! get max params for this element
111  const HomCoord& max_params() const;
112 
113  //! get the number of vertices in each direction, inclusive
114  void param_extents( int& di, int& dj, int& dk ) const;
115 
116  //! given a handle, get the corresponding parameters
117  ErrorCode get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const;
118 
119  //! return whether rectangle is periodic in i
120  inline int is_periodic_i() const
121  {
122  return isPeriodic[0];
123  };
124 
125  //! return whether rectangle is periodic in j
126  inline int is_periodic_j() const
127  {
128  return isPeriodic[1];
129  };
130 
131  inline void is_periodic( int is_p[2] ) const
132  {
133  is_p[0] = isPeriodic[0];
134  is_p[1] = isPeriodic[1];
135  }
136 
137  //! convenience functions for parameter extents
138  int i_min() const
139  {
140  return ( boxParams[0].hom_coord() )[0];
141  }
142  int j_min() const
143  {
144  return ( boxParams[0].hom_coord() )[1];
145  }
146  int k_min() const
147  {
148  return ( boxParams[0].hom_coord() )[2];
149  }
150  int i_max() const
151  {
152  return ( boxParams[1].hom_coord() )[0];
153  }
154  int j_max() const
155  {
156  return ( boxParams[1].hom_coord() )[1];
157  }
158  int k_max() const
159  {
160  return ( boxParams[1].hom_coord() )[2];
161  }
162 
163  //! test the bounding vertex sequences and determine whether they fully
164  //! define the vertices covering this element block's parameter space
165  bool boundary_complete() const;
166 
167  //! test whether this sequence contains these parameters
168  inline bool contains( const HomCoord& coords ) const;
169 
170  //! test whether *vertex parameterization* in this sequence contains these parameters
171  inline bool contains_vertex( const HomCoord& coords ) const;
172 
173  //! get connectivity of an entity given entity's parameters
174  inline ErrorCode get_params_connectivity( const int i,
175  const int j,
176  const int k,
177  std::vector< EntityHandle >& connectivity ) const;
178 
179  //! add a vertex seq ref to this element sequence;
180  //! if bb_input is true, bounding box (in eseq-local coords) of vseq being added
181  //! is input in bb_min and bb_max (allows partial sharing of vseq rather than the whole
182  //! vseq); if it's false, the whole vseq is referenced and the eseq-local coordinates
183  //! is computed from the transformed bounding box of the vseq
185  const HomCoord& p1,
186  const HomCoord& q1,
187  const HomCoord& p2,
188  const HomCoord& q2,
189  const HomCoord& p3,
190  const HomCoord& q3,
191  bool bb_input = false,
192  const HomCoord& bb_min = HomCoord::unitv[0],
193  const HomCoord& bb_max = HomCoord::unitv[0] );
194 
196  EntityHandle end,
197  const int* sequence_data_sizes,
198  const int* tag_data_sizes ) const;
199 
201  int irange,
202  int jrange,
203  int krange,
204  int* is_periodic = NULL );
205 
206  unsigned long get_memory_use() const;
207 };
208 
209 inline EntityHandle ScdElementData::get_element( const int i, const int j, const int k ) const
210 {
211  return start_handle() + ( i - i_min() ) + ( j - j_min() ) * dIJKm1[0] + ( k - k_min() ) * dIJKm1[0] * dIJKm1[1];
212 }
213 
215 {
216  return boxParams[0];
217 }
218 
220 {
221  return boxParams[1];
222 }
223 
224 //! get the number of vertices in each direction, inclusive
225 inline void ScdElementData::param_extents( int& di, int& dj, int& dk ) const
226 {
227  di = dIJK[0];
228  dj = dIJK[1];
229  dk = dIJK[2];
230 }
231 
232 inline ErrorCode ScdElementData::get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const
233 {
234  if( TYPE_FROM_HANDLE( ehandle ) != TYPE_FROM_HANDLE( start_handle() ) ) return MB_FAILURE;
235 
236  int hdiff = ehandle - start_handle();
237 
238  // use double ?: test below because on some platforms, both sides of the : are
239  // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero
240  k = ( dIJKm1[1] > 0 ? hdiff / ( dIJKm1[1] > 0 ? dIJKm1[0] * dIJKm1[1] : 1 ) : 0 );
241  j = ( hdiff - ( k * dIJKm1[0] * dIJKm1[1] ) ) / dIJKm1[0];
242  i = hdiff % dIJKm1[0];
243 
244  k += boxParams[0].k();
245  j += boxParams[0].j();
246  i += boxParams[0].i();
247 
248  return ( ehandle >= start_handle() && ehandle < start_handle() + size() && i >= i_min() && i <= i_max() &&
249  j >= j_min() && j <= j_max() && k >= k_min() && k <= k_max() )
250  ? MB_SUCCESS
251  : MB_FAILURE;
252 }
253 
254 inline bool ScdElementData::contains( const HomCoord& temp ) const
255 {
256  // upper bound is < instead of <= because element params max is one less
257  // than vertex params max, except in case of 2d or 1d sequence
258  return ( ( dIJKm1[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJKm1[0] ) &&
259  ( ( !dIJKm1[1] && temp.j() == boxParams[1].j() ) ||
260  ( dIJKm1[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJKm1[1] ) ) &&
261  ( ( !dIJKm1[2] && temp.k() == boxParams[1].k() ) ||
262  ( dIJKm1[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJKm1[2] ) ) );
263 }
264 
265 inline bool ScdElementData::contains_vertex( const HomCoord& temp ) const
266 {
267  // upper bound is < instead of <= because element params max is one less
268  // than vertex params max, except in case of 2d or 1d sequence
269  return ( ( dIJK[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJK[0] ) &&
270  ( ( !dIJK[1] && temp.j() == boxParams[1].j() ) ||
271  ( dIJK[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJK[1] ) ) &&
272  ( ( !dIJK[2] && temp.k() == boxParams[1].k() ) ||
273  ( dIJK[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJK[2] ) ) );
274 }
275 
276 inline bool ScdElementData::VertexDataRef::contains( const HomCoord& coords ) const
277 {
278  return ( minmax[0] <= coords && minmax[1] >= coords );
279 }
280 
282  const HomCoord& this_max,
283  const HomXform& tmp_xform,
284  ScdVertexData* this_seq )
285  : xform( tmp_xform ), invXform( tmp_xform.inverse() ), srcSeq( this_seq )
286 {
287  minmax[0] = HomCoord( this_min );
288  minmax[1] = HomCoord( this_max );
289 }
290 
291 inline EntityHandle ScdElementData::get_vertex( const HomCoord& coords ) const
292 {
293  assert( boundary_complete() );
294  for( std::vector< VertexDataRef >::const_iterator it = vertexSeqRefs.begin(); it != vertexSeqRefs.end(); ++it )
295  {
296  if( ( *it ).minmax[0] <= coords && ( *it ).minmax[1] >= coords )
297  {
298  // first get the vertex block-local parameters
299  HomCoord local_coords = coords / ( *it ).xform;
300 
301  assert( ( *it ).srcSeq->contains( local_coords ) );
302 
303  // now get the vertex handle for those coords
304  return ( *it ).srcSeq->get_vertex( local_coords );
305  }
306  }
307 
308  // got here, it's an error
309  return 0;
310 }
311 
313  const HomCoord& p1,
314  const HomCoord& q1,
315  const HomCoord& p2,
316  const HomCoord& q2,
317  const HomCoord& p3,
318  const HomCoord& q3,
319  bool bb_input,
320  const HomCoord& bb_min,
321  const HomCoord& bb_max )
322 {
323  // compute the transform given the vseq-local parameters and the mapping to
324  // this element sequence's parameters passed in minmax
325  HomXform M;
326  M.three_pt_xform( p1, q1, p2, q2, p3, q3 );
327 
328  // min and max in element seq's parameter system may not be same as those in
329  // vseq's system, so need to take min/max
330 
331  HomCoord minmax[2];
332  if( bb_input )
333  {
334  minmax[0] = bb_min;
335  minmax[1] = bb_max;
336  }
337  else
338  {
339  minmax[0] = vseq->min_params() * M;
340  minmax[1] = vseq->max_params() * M;
341  }
342 
343  // check against other vseq's to make sure they don't overlap
344  for( std::vector< VertexDataRef >::const_iterator vsit = vertexSeqRefs.begin(); vsit != vertexSeqRefs.end();
345  ++vsit )
346  if( ( *vsit ).contains( minmax[0] ) || ( *vsit ).contains( minmax[1] ) ) return MB_FAILURE;
347 
348  HomCoord tmp_min( std::min( minmax[0].i(), minmax[1].i() ), std::min( minmax[0].j(), minmax[1].j() ),
349  std::min( minmax[0].k(), minmax[1].k() ) );
350  HomCoord tmp_max( std::max( minmax[0].i(), minmax[1].i() ), std::max( minmax[0].j(), minmax[1].j() ),
351  std::max( minmax[0].k(), minmax[1].k() ) );
352 
353  // set up a new vertex sequence reference
354  VertexDataRef tmp_seq_ref( tmp_min, tmp_max, M, vseq );
355 
356  // add to the list
357  vertexSeqRefs.push_back( tmp_seq_ref );
358 
359  return MB_SUCCESS;
360 }
361 
363  const int j,
364  const int k,
365  std::vector< EntityHandle >& connectivity ) const
366 {
367  if( contains( HomCoord( i, j, k ) ) == false ) return MB_FAILURE;
368 
369  int ip1 = ( isPeriodic[0] ? ( i + 1 ) % dIJKm1[0] : i + 1 ),
370  jp1 = ( isPeriodic[1] ? ( j + 1 ) % dIJKm1[1] : j + 1 );
371 
372  connectivity.push_back( get_vertex( i, j, k ) );
373  connectivity.push_back( get_vertex( ip1, j, k ) );
374  if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 2 ) return MB_SUCCESS;
375  connectivity.push_back( get_vertex( ip1, jp1, k ) );
376  connectivity.push_back( get_vertex( i, jp1, k ) );
377  if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 3 ) return MB_SUCCESS;
378  connectivity.push_back( get_vertex( i, j, k + 1 ) );
379  connectivity.push_back( get_vertex( ip1, j, k + 1 ) );
380  connectivity.push_back( get_vertex( ip1, jp1, k + 1 ) );
381  connectivity.push_back( get_vertex( i, jp1, k + 1 ) );
382  return MB_SUCCESS;
383 }
384 
385 } // namespace moab
386 
387 #endif