MOAB: Mesh Oriented datABase  (version 5.3.0)
ScdElementData.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines