MOAB: Mesh Oriented datABase  (version 5.4.1)
SpectralMeshTool.cpp
Go to the documentation of this file.
00001 #include "moab/SpectralMeshTool.hpp"
00002 #include "moab/Interface.hpp"
00003 #include "Internals.hpp"
00004 #include "moab/ReadUtilIface.hpp"
00005 #include "moab/CN.hpp"
00006 
00007 #include <cmath>
00008 #include <cassert>
00009 
00010 namespace moab
00011 {
00012 const short int SpectralMeshTool::permute_array[] = { 0, 1, 13, 25, 3, 2, 14, 26, 7, 6, 18, 30, 11, 10, 22, 34 };
00013 
00014 // lin_permute_array does the same to get the linear vertices of the coarse quad
00015 const short int SpectralMeshTool::lin_permute_array[] = { 0, 25, 34, 11 };
00016 
00017 Tag SpectralMeshTool::spectral_vertices_tag( const bool create_if_missing )
00018 {
00019     ErrorCode rval = MB_SUCCESS;
00020     if( !svTag && create_if_missing )
00021     {
00022         if( !spectralOrder )
00023         {
00024             // should already be a spectral order tag...
00025             MB_SET_ERR_RET_VAL( "Spectral order must be set before creating spectral vertices tag", 0 );
00026         }
00027 
00028         // create it
00029         std::vector< EntityHandle > dum_val( spectralOrderp1 * spectralOrderp1, 0 );
00030         rval = mbImpl->tag_get_handle( "SPECTRAL_VERTICES", spectralOrderp1 * spectralOrderp1, MB_TYPE_HANDLE, svTag,
00031                                        MB_TAG_DENSE | MB_TAG_CREAT, &dum_val[0] );
00032     }
00033 
00034     return ( rval == MB_SUCCESS ? svTag : 0 );
00035 }
00036 
00037 /** \brief Return tag used to store spectral order
00038  * \param create_if_missing If true, will create this tag if it doesn't exist already
00039  */
00040 Tag SpectralMeshTool::spectral_order_tag( const bool create_if_missing )
00041 {
00042     ErrorCode rval = MB_SUCCESS;
00043     if( !soTag && create_if_missing )
00044     {
00045 
00046         // create it
00047         int dum = 0;
00048         rval = mbImpl->tag_get_handle( "SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, soTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
00049     }
00050 
00051     return ( rval == MB_SUCCESS ? soTag : 0 );
00052 }
00053 
00054 /** \brief Convert representation from coarse to fine
00055  * Each element in set, or in interface if set is not input, is converted to fine elements, using
00056  * vertices in SPECTRAL_VERTICES tagged array
00057  * \param spectral_set Set containing spectral elements
00058  */
00059 ErrorCode SpectralMeshTool::convert_to_fine( EntityHandle /* spectral_set */ )
00060 {
00061     return MB_NOT_IMPLEMENTED;
00062 }
00063 
00064 /** \brief Convert representation from fine to coarse
00065  * Each element in set, or in interface if set is not input, is converted to coarse elements, with
00066  * fine vertices put into SPECTRAL_VERTICES tagged array.  NOTE: This function assumes that each
00067  * order^d (fine) elements comprise each coarse element, and are in order of fine elements in each
00068  * coarse element.  If order is input as 0, looks for a SPECTRAL_ORDER tag on the mesh.
00069  * \param order Order of the spectral mesh
00070  * \param spectral_set Set containing spectral elements
00071  */
00072 ErrorCode SpectralMeshTool::convert_to_coarse( int order, EntityHandle spectral_set )
00073 {
00074     if( order ) spectralOrder = order;
00075     if( !spectralOrder )
00076     {
00077         MB_SET_ERR( MB_FAILURE, "Spectral order must be set or input before converting to spectral mesh" );
00078     }
00079 
00080     Range tmp_ents, ents;
00081     ErrorCode rval = mbImpl->get_entities_by_handle( spectral_set, tmp_ents );
00082     if( MB_SUCCESS != rval || ents.empty() ) return rval;
00083 
00084     // get the max-dimensional elements from it
00085     ents = tmp_ents.subset_by_dimension( 3 );
00086     if( ents.empty() ) ents = tmp_ents.subset_by_dimension( 2 );
00087     if( ents.empty() ) ents = tmp_ents.subset_by_dimension( 1 );
00088     if( ents.empty() )
00089     {
00090         MB_SET_ERR( MB_FAILURE, "Can't find any entities for conversion" );
00091     }
00092 
00093     // get a ptr to connectivity
00094     if( ents.psize() != 1 )
00095     {
00096         MB_SET_ERR( MB_FAILURE, "Entities must be in one chunk for conversion" );
00097     }
00098     EntityHandle* conn;
00099     int count, verts_per_e;
00100     rval = mbImpl->connect_iterate( ents.begin(), ents.end(), conn, verts_per_e, count );
00101     if( MB_SUCCESS != rval || count != (int)ents.size() ) return rval;
00102 
00103     Range tmp_range;
00104     return create_spectral_elems( conn, ents.size(), CN::Dimension( TYPE_FROM_HANDLE( *ents.begin() ) ), tmp_range );
00105 }
00106 
00107 template < class T >
00108 ErrorCode SpectralMeshTool::create_spectral_elems( const T* conn,
00109                                                    int num_fine_elems,
00110                                                    int dim,
00111                                                    Range& output_range,
00112                                                    int start_idx,
00113                                                    Range* local_gids )
00114 {
00115     assert( spectralOrder && num_fine_elems );
00116 
00117     // now create num_coarse_elems
00118     // compute the number of local elems, accounting for coarse or fine representation
00119     // spectral_unit is the # fine elems per coarse elem, or spectralOrder^2
00120     int spectral_unit    = spectralOrder * spectralOrder;
00121     int num_coarse_elems = num_fine_elems / spectral_unit;
00122 
00123     EntityHandle* new_conn;
00124     EntityHandle start_elem;
00125     ReadUtilIface* rmi;
00126     ErrorCode rval = mbImpl->query_interface( rmi );
00127     if( MB_SUCCESS != rval ) return rval;
00128 
00129     int verts_per_felem = spectralOrderp1 * spectralOrderp1, verts_per_celem = std::pow( (double)2.0, dim );
00130 
00131     rval = rmi->get_element_connect( num_coarse_elems, verts_per_celem, ( 2 == dim ? MBQUAD : MBHEX ), 0, start_elem,
00132                                      new_conn );MB_CHK_SET_ERR( rval, "Failed to create elems" );
00133 
00134     output_range.insert( start_elem, start_elem + num_coarse_elems - 1 );
00135 
00136     // read connectivity into that space
00137 
00138     // permute_array takes a 4*order^2-long vector of integers, representing the connectivity of
00139     // order^2 elems (fine elems in a coarse elem), and picks out the ids of the vertices necessary
00140     // to get a lexicographically-ordered array of vertices in a spectral element of that order
00141     // assert(verts_per_felem == (sizeof(permute_array)/sizeof(unsigned int)));
00142 
00143     // we're assuming here that elems was empty on input
00144     int count;
00145     EntityHandle* sv_ptr = NULL;
00146     rval = mbImpl->tag_iterate( spectral_vertices_tag( true ), output_range.begin(), output_range.end(), count,
00147                                 (void*&)sv_ptr );MB_CHK_SET_ERR( rval, "Failed to get SPECTRAL_VERTICES ptr" );
00148     assert( count == num_coarse_elems );
00149     int f = start_idx, fs = 0, fl = 0;
00150     for( int c = 0; c < num_coarse_elems; c++ )
00151     {
00152         for( int i = 0; i < verts_per_celem; i++ )
00153             new_conn[fl + i] = conn[f + lin_permute_array[i]];
00154         fl += verts_per_celem;
00155         for( int i = 0; i < verts_per_felem; i++ )
00156             sv_ptr[fs + i] = conn[f + permute_array[i]];
00157         f += verts_per_celem * spectral_unit;
00158         fs += verts_per_felem;
00159     }
00160     if( local_gids ) std::copy( sv_ptr, sv_ptr + verts_per_felem * num_coarse_elems, range_inserter( *local_gids ) );
00161 
00162     return MB_SUCCESS;
00163 }
00164 
00165 // force instantiation of a few specific types
00166 template ErrorCode SpectralMeshTool::create_spectral_elems< int >( const int* conn,
00167                                                                    int num_fine_elems,
00168                                                                    int dim,
00169                                                                    Range& output_range,
00170                                                                    int start_idx,
00171                                                                    Range* local_gids );
00172 template ErrorCode SpectralMeshTool::create_spectral_elems< EntityHandle >( const EntityHandle* conn,
00173                                                                             int num_fine_elems,
00174                                                                             int dim,
00175                                                                             Range& output_range,
00176                                                                             int start_idx,
00177                                                                             Range* local_gids );
00178 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines