Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00087                     const int imin,
00088                     const int jmin,
00089                     const int kmin,
00090                     const int imax,
00091                     const int jmax,
00092                     const int kmax,
00093                     int* is_periodic );
00094 
00095     virtual ~ScdElementData();
00096 
00097     //! get handle of vertex at homogeneous coords
00098     inline EntityHandle get_vertex( const HomCoord& coords ) const;
00099     inline EntityHandle get_vertex( int i, int j, int k ) const
00100     {
00101         return get_vertex( HomCoord( i, j, k ) );
00102     }
00103 
00104     //! get handle of element at i, j, k
00105     EntityHandle get_element( const int i, const int j, const int k ) const;
00106 
00107     //! get min params for this element
00108     const HomCoord& min_params() const;
00109 
00110     //! get max params for this element
00111     const HomCoord& max_params() const;
00112 
00113     //! get the number of vertices in each direction, inclusive
00114     void param_extents( int& di, int& dj, int& dk ) const;
00115 
00116     //! given a handle, get the corresponding parameters
00117     ErrorCode get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const;
00118 
00119     //! return whether rectangle is periodic in i
00120     inline int is_periodic_i() const
00121     {
00122         return isPeriodic[0];
00123     };
00124 
00125     //! return whether rectangle is periodic in j
00126     inline int is_periodic_j() const
00127     {
00128         return isPeriodic[1];
00129     };
00130 
00131     inline void is_periodic( int is_p[2] ) const
00132     {
00133         is_p[0] = isPeriodic[0];
00134         is_p[1] = isPeriodic[1];
00135     }
00136 
00137     //! convenience functions for parameter extents
00138     int i_min() const
00139     {
00140         return ( boxParams[0].hom_coord() )[0];
00141     }
00142     int j_min() const
00143     {
00144         return ( boxParams[0].hom_coord() )[1];
00145     }
00146     int k_min() const
00147     {
00148         return ( boxParams[0].hom_coord() )[2];
00149     }
00150     int i_max() const
00151     {
00152         return ( boxParams[1].hom_coord() )[0];
00153     }
00154     int j_max() const
00155     {
00156         return ( boxParams[1].hom_coord() )[1];
00157     }
00158     int k_max() const
00159     {
00160         return ( boxParams[1].hom_coord() )[2];
00161     }
00162 
00163     //! test the bounding vertex sequences and determine whether they fully
00164     //! define the vertices covering this element block's parameter space
00165     bool boundary_complete() const;
00166 
00167     //! test whether this sequence contains these parameters
00168     inline bool contains( const HomCoord& coords ) const;
00169 
00170     //! test whether *vertex parameterization* in this sequence contains these parameters
00171     inline bool contains_vertex( const HomCoord& coords ) const;
00172 
00173     //! get connectivity of an entity given entity's parameters
00174     inline ErrorCode get_params_connectivity( const int i,
00175                                               const int j,
00176                                               const int k,
00177                                               std::vector< EntityHandle >& connectivity ) const;
00178 
00179     //! add a vertex seq ref to this element sequence;
00180     //! if bb_input is true, bounding box (in eseq-local coords) of vseq being added
00181     //! is input in bb_min and bb_max (allows partial sharing of vseq rather than the whole
00182     //! vseq); if it's false, the whole vseq is referenced and the eseq-local coordinates
00183     //! is computed from the transformed bounding box of the vseq
00184     ErrorCode add_vsequence( ScdVertexData* vseq,
00185                              const HomCoord& p1,
00186                              const HomCoord& q1,
00187                              const HomCoord& p2,
00188                              const HomCoord& q2,
00189                              const HomCoord& p3,
00190                              const HomCoord& q3,
00191                              bool bb_input          = false,
00192                              const HomCoord& bb_min = HomCoord::unitv[0],
00193                              const HomCoord& bb_max = HomCoord::unitv[0] );
00194 
00195     SequenceData* subset( EntityHandle start,
00196                           EntityHandle end,
00197                           const int* sequence_data_sizes,
00198                           const int* tag_data_sizes ) const;
00199 
00200     static EntityID calc_num_entities( EntityHandle start_handle,
00201                                        int irange,
00202                                        int jrange,
00203                                        int krange,
00204                                        int* is_periodic = NULL );
00205 
00206     unsigned long get_memory_use() const;
00207 };
00208 
00209 inline EntityHandle ScdElementData::get_element( const int i, const int j, const int k ) const
00210 {
00211     return start_handle() + ( i - i_min() ) + ( j - j_min() ) * dIJKm1[0] + ( k - k_min() ) * dIJKm1[0] * dIJKm1[1];
00212 }
00213 
00214 inline const HomCoord& ScdElementData::min_params() const
00215 {
00216     return boxParams[0];
00217 }
00218 
00219 inline const HomCoord& ScdElementData::max_params() const
00220 {
00221     return boxParams[1];
00222 }
00223 
00224 //! get the number of vertices in each direction, inclusive
00225 inline void ScdElementData::param_extents( int& di, int& dj, int& dk ) const
00226 {
00227     di = dIJK[0];
00228     dj = dIJK[1];
00229     dk = dIJK[2];
00230 }
00231 
00232 inline ErrorCode ScdElementData::get_params( const EntityHandle ehandle, int& i, int& j, int& k ) const
00233 {
00234     if( TYPE_FROM_HANDLE( ehandle ) != TYPE_FROM_HANDLE( start_handle() ) ) return MB_FAILURE;
00235 
00236     int hdiff = ehandle - start_handle();
00237 
00238     // use double ?: test below because on some platforms, both sides of the : are
00239     // evaluated, and if dIJKm1[1] is zero, that'll generate a divide-by-zero
00240     k = ( dIJKm1[1] > 0 ? hdiff / ( dIJKm1[1] > 0 ? dIJKm1[0] * dIJKm1[1] : 1 ) : 0 );
00241     j = ( hdiff - ( k * dIJKm1[0] * dIJKm1[1] ) ) / dIJKm1[0];
00242     i = hdiff % dIJKm1[0];
00243 
00244     k += boxParams[0].k();
00245     j += boxParams[0].j();
00246     i += boxParams[0].i();
00247 
00248     return ( ehandle >= start_handle() && ehandle < start_handle() + size() && i >= i_min() && i <= i_max() &&
00249              j >= j_min() && j <= j_max() && k >= k_min() && k <= k_max() )
00250                ? MB_SUCCESS
00251                : MB_FAILURE;
00252 }
00253 
00254 inline bool ScdElementData::contains( const HomCoord& temp ) const
00255 {
00256     // upper bound is < instead of <= because element params max is one less
00257     // than vertex params max, except in case of 2d or 1d sequence
00258     return ( ( dIJKm1[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJKm1[0] ) &&
00259              ( ( !dIJKm1[1] && temp.j() == boxParams[1].j() ) ||
00260                ( dIJKm1[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJKm1[1] ) ) &&
00261              ( ( !dIJKm1[2] && temp.k() == boxParams[1].k() ) ||
00262                ( dIJKm1[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJKm1[2] ) ) );
00263 }
00264 
00265 inline bool ScdElementData::contains_vertex( const HomCoord& temp ) const
00266 {
00267     // upper bound is < instead of <= because element params max is one less
00268     // than vertex params max, except in case of 2d or 1d sequence
00269     return ( ( dIJK[0] && temp.i() >= boxParams[0].i() && temp.i() < boxParams[0].i() + dIJK[0] ) &&
00270              ( ( !dIJK[1] && temp.j() == boxParams[1].j() ) ||
00271                ( dIJK[1] && temp.j() >= boxParams[0].j() && temp.j() < boxParams[0].j() + dIJK[1] ) ) &&
00272              ( ( !dIJK[2] && temp.k() == boxParams[1].k() ) ||
00273                ( dIJK[2] && temp.k() >= boxParams[0].k() && temp.k() < boxParams[0].k() + dIJK[2] ) ) );
00274 }
00275 
00276 inline bool ScdElementData::VertexDataRef::contains( const HomCoord& coords ) const
00277 {
00278     return ( minmax[0] <= coords && minmax[1] >= coords );
00279 }
00280 
00281 inline ScdElementData::VertexDataRef::VertexDataRef( const HomCoord& this_min,
00282                                                      const HomCoord& this_max,
00283                                                      const HomXform& tmp_xform,
00284                                                      ScdVertexData* this_seq )
00285     : xform( tmp_xform ), invXform( tmp_xform.inverse() ), srcSeq( this_seq )
00286 {
00287     minmax[0] = HomCoord( this_min );
00288     minmax[1] = HomCoord( this_max );
00289 }
00290 
00291 inline EntityHandle ScdElementData::get_vertex( const HomCoord& coords ) const
00292 {
00293     assert( boundary_complete() );
00294     for( std::vector< VertexDataRef >::const_iterator it = vertexSeqRefs.begin(); it != vertexSeqRefs.end(); ++it )
00295     {
00296         if( ( *it ).minmax[0] <= coords && ( *it ).minmax[1] >= coords )
00297         {
00298             // first get the vertex block-local parameters
00299             HomCoord local_coords = coords / ( *it ).xform;
00300 
00301             assert( ( *it ).srcSeq->contains( local_coords ) );
00302 
00303             // now get the vertex handle for those coords
00304             return ( *it ).srcSeq->get_vertex( local_coords );
00305         }
00306     }
00307 
00308     // got here, it's an error
00309     return 0;
00310 }
00311 
00312 inline ErrorCode ScdElementData::add_vsequence( ScdVertexData* vseq,
00313                                                 const HomCoord& p1,
00314                                                 const HomCoord& q1,
00315                                                 const HomCoord& p2,
00316                                                 const HomCoord& q2,
00317                                                 const HomCoord& p3,
00318                                                 const HomCoord& q3,
00319                                                 bool bb_input,
00320                                                 const HomCoord& bb_min,
00321                                                 const HomCoord& bb_max )
00322 {
00323     // compute the transform given the vseq-local parameters and the mapping to
00324     // this element sequence's parameters passed in minmax
00325     HomXform M;
00326     M.three_pt_xform( p1, q1, p2, q2, p3, q3 );
00327 
00328     // min and max in element seq's parameter system may not be same as those in
00329     // vseq's system, so need to take min/max
00330 
00331     HomCoord minmax[2];
00332     if( bb_input )
00333     {
00334         minmax[0] = bb_min;
00335         minmax[1] = bb_max;
00336     }
00337     else
00338     {
00339         minmax[0] = vseq->min_params() * M;
00340         minmax[1] = vseq->max_params() * M;
00341     }
00342 
00343     // check against other vseq's to make sure they don't overlap
00344     for( std::vector< VertexDataRef >::const_iterator vsit = vertexSeqRefs.begin(); vsit != vertexSeqRefs.end();
00345          ++vsit )
00346         if( ( *vsit ).contains( minmax[0] ) || ( *vsit ).contains( minmax[1] ) ) return MB_FAILURE;
00347 
00348     HomCoord tmp_min( std::min( minmax[0].i(), minmax[1].i() ), std::min( minmax[0].j(), minmax[1].j() ),
00349                       std::min( minmax[0].k(), minmax[1].k() ) );
00350     HomCoord tmp_max( std::max( minmax[0].i(), minmax[1].i() ), std::max( minmax[0].j(), minmax[1].j() ),
00351                       std::max( minmax[0].k(), minmax[1].k() ) );
00352 
00353     // set up a new vertex sequence reference
00354     VertexDataRef tmp_seq_ref( tmp_min, tmp_max, M, vseq );
00355 
00356     // add to the list
00357     vertexSeqRefs.push_back( tmp_seq_ref );
00358 
00359     return MB_SUCCESS;
00360 }
00361 
00362 inline ErrorCode ScdElementData::get_params_connectivity( const int i,
00363                                                           const int j,
00364                                                           const int k,
00365                                                           std::vector< EntityHandle >& connectivity ) const
00366 {
00367     if( contains( HomCoord( i, j, k ) ) == false ) return MB_FAILURE;
00368 
00369     int ip1 = ( isPeriodic[0] ? ( i + 1 ) % dIJKm1[0] : i + 1 ),
00370         jp1 = ( isPeriodic[1] ? ( j + 1 ) % dIJKm1[1] : j + 1 );
00371 
00372     connectivity.push_back( get_vertex( i, j, k ) );
00373     connectivity.push_back( get_vertex( ip1, j, k ) );
00374     if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 2 ) return MB_SUCCESS;
00375     connectivity.push_back( get_vertex( ip1, jp1, k ) );
00376     connectivity.push_back( get_vertex( i, jp1, k ) );
00377     if( CN::Dimension( TYPE_FROM_HANDLE( start_handle() ) ) < 3 ) return MB_SUCCESS;
00378     connectivity.push_back( get_vertex( i, j, k + 1 ) );
00379     connectivity.push_back( get_vertex( ip1, j, k + 1 ) );
00380     connectivity.push_back( get_vertex( ip1, jp1, k + 1 ) );
00381     connectivity.push_back( get_vertex( i, jp1, k + 1 ) );
00382     return MB_SUCCESS;
00383 }
00384 
00385 }  // namespace moab
00386 
00387 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines