MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /** 00002 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 00003 * storing and accessing finite element mesh data. 00004 * 00005 * Copyright 2004 Sandia Corporation. Under the terms of Contract 00006 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00007 * retains certain rights in this software. 00008 * 00009 * This library is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU Lesser General Public 00011 * License as published by the Free Software Foundation; either 00012 * version 2.1 of the License, or (at your option) any later version. 00013 * 00014 */ 00015 00016 #ifndef SCD_ELEMENT_DATA_HPP 00017 #define SCD_ELEMENT_DATA_HPP 00018 00019 // 00020 // Class: ScdElementData 00021 // 00022 // Purpose: represent a rectangular element of mesh 00023 // 00024 // A ScdElement represents a rectangular element of mesh, including both vertices and 00025 // elements, and the parametric space used to address that element. Vertex data, 00026 // i.e. coordinates, may not be stored directly in the element, but the element returns 00027 // information about the vertex handles of vertices in the element. Vertex and element 00028 // handles associated with the element are each contiguous. 00029 00030 #include "SequenceData.hpp" 00031 #include "moab/HomXform.hpp" 00032 #include "moab/CN.hpp" 00033 #include "ScdVertexData.hpp" 00034 #include "Internals.hpp" 00035 #include "moab/Range.hpp" 00036 00037 #include <vector> 00038 #include <algorithm> 00039 00040 namespace moab 00041 { 00042 00043 class ScdElementData : public SequenceData 00044 { 00045 00046 //! structure to hold references to bounding vertex blocks 00047 class VertexDataRef 00048 { 00049 private: 00050 HomCoord minmax[2]; 00051 HomXform xform, invXform; 00052 ScdVertexData* srcSeq; 00053 00054 public: 00055 friend class ScdElementData; 00056 00057 VertexDataRef( const HomCoord& min, const HomCoord& max, const HomXform& tmp_xform, ScdVertexData* this_seq ); 00058 00059 bool contains( const HomCoord& coords ) const; 00060 }; 00061 00062 private: 00063 //! parameter min/max/stride for vertices, in homogeneous coords ijkh; elements 00064 //! are one less than max 00065 HomCoord boxParams[3]; 00066 00067 //! difference between max and min params plus one (i.e. # VERTICES in 00068 //! each parametric direction) 00069 int dIJK[3]; 00070 00071 //! difference between max and min params (i.e. # ELEMENTS in 00072 //! each parametric direction) 00073 int dIJKm1[3]; 00074 00075 //! whether scd element rectangle is periodic in i and possibly j 00076 int isPeriodic[2]; 00077 00078 //! bare constructor, so compiler doesn't create one for me 00079 ScdElementData(); 00080 00081 //! list of bounding vertex blocks 00082 std::vector< VertexDataRef > vertexSeqRefs; 00083 00084 public: 00085 //! constructor 00086 ScdElementData( EntityHandle start_handle, const int imin, const int jmin, const int kmin, const int imax, 00087 const int jmax, const int kmax, int* is_periodic ); 00088 00089 virtual ~ScdElementData(); 00090 00091 //! get handle of vertex at homogeneous coords 00092 inline EntityHandle get_vertex( const HomCoord& coords ) const; 00093 inline EntityHandle get_vertex( int i, int j, int k ) const 00094 { 00095 return get_vertex( HomCoord( i, j, k ) ); 00096 } 00097 00098 //! get handle of element at i, j, k 00099 EntityHandle get_element( const int i, const int j, const int k ) const; 00100 00101 //! get min params for this element 00102 const HomCoord& min_params() const; 00103 00104 //! get max params for this element 00105 const HomCoord& max_params() const; 00106 00107 //! get the number of vertices in each direction, inclusive 00108 void param_extents( int& di, int& dj, int& dk ) const; 00109 00110 //! given a handle, get the corresponding parameters 00111 ErrorCode get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const; 00112 00113 //! return whether rectangle is periodic in i 00114 inline int is_periodic_i() const 00115 { 00116 return isPeriodic[0]; 00117 }; 00118 00119 //! return whether rectangle is periodic in j 00120 inline int is_periodic_j() const 00121 { 00122 return isPeriodic[1]; 00123 }; 00124 00125 inline void is_periodic( int is_p[2] ) const 00126 { 00127 is_p[0] = isPeriodic[0]; 00128 is_p[1] = isPeriodic[1]; 00129 } 00130 00131 //! convenience functions for parameter extents 00132 int i_min() const 00133 { 00134 return ( boxParams[0].hom_coord() )[0]; 00135 } 00136 int j_min() const 00137 { 00138 return ( boxParams[0].hom_coord() )[1]; 00139 } 00140 int k_min() const 00141 { 00142 return ( boxParams[0].hom_coord() )[2]; 00143 } 00144 int i_max() const 00145 { 00146 return ( boxParams[1].hom_coord() )[0]; 00147 } 00148 int j_max() const 00149 { 00150 return ( boxParams[1].hom_coord() )[1]; 00151 } 00152 int k_max() const 00153 { 00154 return ( boxParams[1].hom_coord() )[2]; 00155 } 00156 00157 //! test the bounding vertex sequences and determine whether they fully 00158 //! define the vertices covering this element block's parameter space 00159 bool boundary_complete() const; 00160 00161 //! test whether this sequence contains these parameters 00162 inline bool contains( const HomCoord& coords ) const; 00163 00164 //! test whether *vertex parameterization* in this sequence contains these parameters 00165 inline bool contains_vertex( const HomCoord& coords ) const; 00166 00167 //! get connectivity of an entity given entity's parameters 00168 inline ErrorCode get_params_connectivity( const int i, const int j, const int k, 00169 std::vector< EntityHandle >& connectivity ) const; 00170 00171 //! add a vertex seq ref to this element sequence; 00172 //! if bb_input is true, bounding box (in eseq-local coords) of vseq being added 00173 //! is input in bb_min and bb_max (allows partial sharing of vseq rather than the whole 00174 //! vseq); if it's false, the whole vseq is referenced and the eseq-local coordinates 00175 //! is computed from the transformed bounding box of the vseq 00176 ErrorCode add_vsequence( ScdVertexData* vseq, const HomCoord& p1, const HomCoord& q1, const HomCoord& p2, 00177 const HomCoord& q2, const HomCoord& p3, const HomCoord& q3, bool bb_input = false, 00178 const HomCoord& bb_min = HomCoord::unitv[0], const HomCoord& bb_max = HomCoord::unitv[0] ); 00179 00180 SequenceData* subset( EntityHandle start, EntityHandle end, const int* sequence_data_sizes, 00181 const int* tag_data_sizes ) const; 00182 00183 static EntityID calc_num_entities( EntityHandle start_handle, int irange, int jrange, int krange, 00184 int* is_periodic = NULL ); 00185 00186 unsigned long get_memory_use() const; 00187 }; 00188 00189 inline EntityHandle ScdElementData::get_element( const int i, const int j, const int k ) const 00190 { 00191 return start_handle() + ( i - i_min() ) + ( j - j_min() ) * dIJKm1[0] + ( k - k_min() ) * dIJKm1[0] * dIJKm1[1]; 00192 } 00193 00194 inline const HomCoord& ScdElementData::min_params() const 00195 { 00196 return boxParams[0]; 00197 } 00198 00199 inline const HomCoord& ScdElementData::max_params() const 00200 { 00201 return boxParams[1]; 00202 } 00203 00204 //! get the number of vertices in each direction, inclusive 00205 inline void ScdElementData::param_extents( int& di, int& dj, int& dk ) const 00206 { 00207 di = dIJK[0]; 00208 dj = dIJK[1]; 00209 dk = dIJK[2]; 00210 } 00211 00212 inline ErrorCode ScdElementData::get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const 00213 { 00214 if( TYPE_FROM_HANDLE( ehandle ) != TYPE_FROM_HANDLE( start_handle() ) ) return MB_FAILURE; 00215 00216 int hdiff = ehandle - start_handle(); 00217 00218 // use double ?: test below because on some platforms, both sides of the : are 00219 // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero 00220 k = ( dIJKm1[1] > 0 ? hdiff / ( dIJKm1[1] > 0 ? dIJKm1[0] * dIJKm1[1] : 1 ) : 0 ); 00221 j = ( hdiff - ( k * dIJKm1[0] * dIJKm1[1] ) ) / dIJKm1[0]; 00222 i = hdiff % dIJKm1[0]; 00223 00224 k += boxParams[0].k(); 00225 j += boxParams[0].j(); 00226 i += boxParams[0].i(); 00227 00228 return ( ehandle >= start_handle() && ehandle < start_handle() + size() && i >= i_min() && i <= i_max() && 00229 j >= j_min() && j <= j_max() && k >= k_min() && k <= k_max() ) 00230 ? MB_SUCCESS 00231 : MB_FAILURE; 00232 } 00233 00234 inline bool ScdElementData::contains( const HomCoord& temp ) const 00235 { 00236 // upper bound is < instead of <= because element params max is one less 00237 // than vertex params max, except in case of 2d or 1d sequence 00238 return ( ( dIJKm1[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJKm1[0] ) && 00239 ( ( !dIJKm1[1] && temp.j() == boxParams[1].j() ) || 00240 ( dIJKm1[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJKm1[1] ) ) && 00241 ( ( !dIJKm1[2] && temp.k() == boxParams[1].k() ) || 00242 ( dIJKm1[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJKm1[2] ) ) ); 00243 } 00244 00245 inline bool ScdElementData::contains_vertex( const HomCoord& temp ) const 00246 { 00247 // upper bound is < instead of <= because element params max is one less 00248 // than vertex params max, except in case of 2d or 1d sequence 00249 return ( ( dIJK[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJK[0] ) && 00250 ( ( !dIJK[1] && temp.j() == boxParams[1].j() ) || 00251 ( dIJK[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJK[1] ) ) && 00252 ( ( !dIJK[2] && temp.k() == boxParams[1].k() ) || 00253 ( dIJK[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJK[2] ) ) ); 00254 } 00255 00256 inline bool ScdElementData::VertexDataRef::contains( const HomCoord& coords ) const 00257 { 00258 return ( minmax[0] <= coords && minmax[1] >= coords ); 00259 } 00260 00261 inline ScdElementData::VertexDataRef::VertexDataRef( const HomCoord& this_min, const HomCoord& this_max, 00262 const HomXform& tmp_xform, ScdVertexData* this_seq ) 00263 : xform( tmp_xform ), invXform( tmp_xform.inverse() ), srcSeq( this_seq ) 00264 { 00265 minmax[0] = HomCoord( this_min ); 00266 minmax[1] = HomCoord( this_max ); 00267 } 00268 00269 inline EntityHandle ScdElementData::get_vertex( const HomCoord& coords ) const 00270 { 00271 assert( boundary_complete() ); 00272 for( std::vector< VertexDataRef >::const_iterator it = vertexSeqRefs.begin(); it != vertexSeqRefs.end(); ++it ) 00273 { 00274 if( ( *it ).minmax[0] <= coords && ( *it ).minmax[1] >= coords ) 00275 { 00276 // first get the vertex block-local parameters 00277 HomCoord local_coords = coords / ( *it ).xform; 00278 00279 assert( ( *it ).srcSeq->contains( local_coords ) ); 00280 00281 // now get the vertex handle for those coords 00282 return ( *it ).srcSeq->get_vertex( local_coords ); 00283 } 00284 } 00285 00286 // got here, it's an error 00287 return 0; 00288 } 00289 00290 inline ErrorCode ScdElementData::add_vsequence( ScdVertexData* vseq, const HomCoord& p1, const HomCoord& q1, 00291 const HomCoord& p2, const HomCoord& q2, const HomCoord& p3, 00292 const HomCoord& q3, bool bb_input, const HomCoord& bb_min, 00293 const HomCoord& bb_max ) 00294 { 00295 // compute the transform given the vseq-local parameters and the mapping to 00296 // this element sequence's parameters passed in minmax 00297 HomXform M; 00298 M.three_pt_xform( p1, q1, p2, q2, p3, q3 ); 00299 00300 // min and max in element seq's parameter system may not be same as those in 00301 // vseq's system, so need to take min/max 00302 00303 HomCoord minmax[2]; 00304 if( bb_input ) 00305 { 00306 minmax[0] = bb_min; 00307 minmax[1] = bb_max; 00308 } 00309 else 00310 { 00311 minmax[0] = vseq->min_params() * M; 00312 minmax[1] = vseq->max_params() * M; 00313 } 00314 00315 // check against other vseq's to make sure they don't overlap 00316 for( std::vector< VertexDataRef >::const_iterator vsit = vertexSeqRefs.begin(); vsit != vertexSeqRefs.end(); 00317 ++vsit ) 00318 if( ( *vsit ).contains( minmax[0] ) || ( *vsit ).contains( minmax[1] ) ) return MB_FAILURE; 00319 00320 HomCoord tmp_min( std::min( minmax[0].i(), minmax[1].i() ), std::min( minmax[0].j(), minmax[1].j() ), 00321 std::min( minmax[0].k(), minmax[1].k() ) ); 00322 HomCoord tmp_max( std::max( minmax[0].i(), minmax[1].i() ), std::max( minmax[0].j(), minmax[1].j() ), 00323 std::max( minmax[0].k(), minmax[1].k() ) ); 00324 00325 // set up a new vertex sequence reference 00326 VertexDataRef tmp_seq_ref( tmp_min, tmp_max, M, vseq ); 00327 00328 // add to the list 00329 vertexSeqRefs.push_back( tmp_seq_ref ); 00330 00331 return MB_SUCCESS; 00332 } 00333 00334 inline ErrorCode ScdElementData::get_params_connectivity( const int i, const int j, const int k, 00335 std::vector< EntityHandle >& connectivity ) const 00336 { 00337 if( contains( HomCoord( i, j, k ) ) == false ) return MB_FAILURE; 00338 00339 int ip1 = ( isPeriodic[0] ? ( i + 1 ) % dIJKm1[0] : i + 1 ), 00340 jp1 = ( isPeriodic[1] ? ( j + 1 ) % dIJKm1[1] : j + 1 ); 00341 00342 connectivity.push_back( get_vertex( i, j, k ) ); 00343 connectivity.push_back( get_vertex( ip1, j, k ) ); 00344 if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 2 ) return MB_SUCCESS; 00345 connectivity.push_back( get_vertex( ip1, jp1, k ) ); 00346 connectivity.push_back( get_vertex( i, jp1, k ) ); 00347 if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 3 ) return MB_SUCCESS; 00348 connectivity.push_back( get_vertex( i, j, k + 1 ) ); 00349 connectivity.push_back( get_vertex( ip1, j, k + 1 ) ); 00350 connectivity.push_back( get_vertex( ip1, jp1, k + 1 ) ); 00351 connectivity.push_back( get_vertex( i, jp1, k + 1 ) ); 00352 return MB_SUCCESS; 00353 } 00354 00355 } // namespace moab 00356 00357 #endif