MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 <assert.h> 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 { MB_SET_ERR( MB_FAILURE, "Spectral order must be set or input before converting to spectral mesh" ); } 00077 00078 Range tmp_ents, ents; 00079 ErrorCode rval = mbImpl->get_entities_by_handle( spectral_set, tmp_ents ); 00080 if( MB_SUCCESS != rval || ents.empty() ) return rval; 00081 00082 // get the max-dimensional elements from it 00083 ents = tmp_ents.subset_by_dimension( 3 ); 00084 if( ents.empty() ) ents = tmp_ents.subset_by_dimension( 2 ); 00085 if( ents.empty() ) ents = tmp_ents.subset_by_dimension( 1 ); 00086 if( ents.empty() ) { MB_SET_ERR( MB_FAILURE, "Can't find any entities for conversion" ); } 00087 00088 // get a ptr to connectivity 00089 if( ents.psize() != 1 ) { MB_SET_ERR( MB_FAILURE, "Entities must be in one chunk for conversion" ); } 00090 EntityHandle* conn; 00091 int count, verts_per_e; 00092 rval = mbImpl->connect_iterate( ents.begin(), ents.end(), conn, verts_per_e, count ); 00093 if( MB_SUCCESS != rval || count != (int)ents.size() ) return rval; 00094 00095 Range tmp_range; 00096 return create_spectral_elems( conn, ents.size(), CN::Dimension( TYPE_FROM_HANDLE( *ents.begin() ) ), tmp_range ); 00097 } 00098 00099 template < class T > 00100 ErrorCode SpectralMeshTool::create_spectral_elems( const T* conn, int num_fine_elems, int dim, Range& output_range, 00101 int start_idx, Range* local_gids ) 00102 { 00103 assert( spectralOrder && num_fine_elems ); 00104 00105 // now create num_coarse_elems 00106 // compute the number of local elems, accounting for coarse or fine representation 00107 // spectral_unit is the # fine elems per coarse elem, or spectralOrder^2 00108 int spectral_unit = spectralOrder * spectralOrder; 00109 int num_coarse_elems = num_fine_elems / spectral_unit; 00110 00111 EntityHandle* new_conn; 00112 EntityHandle start_elem; 00113 ReadUtilIface* rmi; 00114 ErrorCode rval = mbImpl->query_interface( rmi ); 00115 if( MB_SUCCESS != rval ) return rval; 00116 00117 int verts_per_felem = spectralOrderp1 * spectralOrderp1, verts_per_celem = std::pow( (double)2.0, dim ); 00118 00119 rval = rmi->get_element_connect( num_coarse_elems, verts_per_celem, ( 2 == dim ? MBQUAD : MBHEX ), 0, start_elem, 00120 new_conn );MB_CHK_SET_ERR( rval, "Failed to create elems" ); 00121 00122 output_range.insert( start_elem, start_elem + num_coarse_elems - 1 ); 00123 00124 // read connectivity into that space 00125 00126 // permute_array takes a 4*order^2-long vector of integers, representing the connectivity of 00127 // order^2 elems (fine elems in a coarse elem), and picks out the ids of the vertices necessary 00128 // to get a lexicographically-ordered array of vertices in a spectral element of that order 00129 // assert(verts_per_felem == (sizeof(permute_array)/sizeof(unsigned int))); 00130 00131 // we're assuming here that elems was empty on input 00132 int count; 00133 EntityHandle* sv_ptr = NULL; 00134 rval = mbImpl->tag_iterate( spectral_vertices_tag( true ), output_range.begin(), output_range.end(), count, 00135 (void*&)sv_ptr );MB_CHK_SET_ERR( rval, "Failed to get SPECTRAL_VERTICES ptr" ); 00136 assert( count == num_coarse_elems ); 00137 int f = start_idx, fs = 0, fl = 0; 00138 for( int c = 0; c < num_coarse_elems; c++ ) 00139 { 00140 for( int i = 0; i < verts_per_celem; i++ ) 00141 new_conn[fl + i] = conn[f + lin_permute_array[i]]; 00142 fl += verts_per_celem; 00143 for( int i = 0; i < verts_per_felem; i++ ) 00144 sv_ptr[fs + i] = conn[f + permute_array[i]]; 00145 f += verts_per_celem * spectral_unit; 00146 fs += verts_per_felem; 00147 } 00148 if( local_gids ) std::copy( sv_ptr, sv_ptr + verts_per_felem * num_coarse_elems, range_inserter( *local_gids ) ); 00149 00150 return MB_SUCCESS; 00151 } 00152 00153 // force instantiation of a few specific types 00154 template ErrorCode SpectralMeshTool::create_spectral_elems< int >( const int* conn, int num_fine_elems, int dim, 00155 Range& output_range, int start_idx, 00156 Range* local_gids ); 00157 template ErrorCode SpectralMeshTool::create_spectral_elems< EntityHandle >( const EntityHandle* conn, 00158 int num_fine_elems, int dim, 00159 Range& output_range, int start_idx, 00160 Range* local_gids ); 00161 } // namespace moab