Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
SpectralMeshTool.cpp
Go to the documentation of this file.
2 #include "moab/Interface.hpp"
3 #include "Internals.hpp"
4 #include "moab/ReadUtilIface.hpp"
5 #include "moab/CN.hpp"
6 
7 #include <cmath>
8 #include <cassert>
9 
10 namespace moab
11 {
12 const short int SpectralMeshTool::permute_array[] = { 0, 1, 13, 25, 3, 2, 14, 26, 7, 6, 18, 30, 11, 10, 22, 34 };
13 
14 // lin_permute_array does the same to get the linear vertices of the coarse quad
15 const short int SpectralMeshTool::lin_permute_array[] = { 0, 25, 34, 11 };
16 
17 Tag SpectralMeshTool::spectral_vertices_tag( const bool create_if_missing )
18 {
19  ErrorCode rval = MB_SUCCESS;
20  if( !svTag && create_if_missing )
21  {
22  if( !spectralOrder )
23  {
24  // should already be a spectral order tag...
25  MB_SET_ERR_RET_VAL( "Spectral order must be set before creating spectral vertices tag", 0 );
26  }
27 
28  // create it
29  std::vector< EntityHandle > dum_val( spectralOrderp1 * spectralOrderp1, 0 );
30  rval = mbImpl->tag_get_handle( "SPECTRAL_VERTICES", spectralOrderp1 * spectralOrderp1, MB_TYPE_HANDLE, svTag,
31  MB_TAG_DENSE | MB_TAG_CREAT, &dum_val[0] );
32  }
33 
34  return ( rval == MB_SUCCESS ? svTag : 0 );
35 }
36 
37 /** \brief Return tag used to store spectral order
38  * \param create_if_missing If true, will create this tag if it doesn't exist already
39  */
40 Tag SpectralMeshTool::spectral_order_tag( const bool create_if_missing )
41 {
42  ErrorCode rval = MB_SUCCESS;
43  if( !soTag && create_if_missing )
44  {
45 
46  // create it
47  int dum = 0;
48  rval = mbImpl->tag_get_handle( "SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, soTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
49  }
50 
51  return ( rval == MB_SUCCESS ? soTag : 0 );
52 }
53 
54 /** \brief Convert representation from coarse to fine
55  * Each element in set, or in interface if set is not input, is converted to fine elements, using
56  * vertices in SPECTRAL_VERTICES tagged array
57  * \param spectral_set Set containing spectral elements
58  */
60 {
61  return MB_NOT_IMPLEMENTED;
62 }
63 
64 /** \brief Convert representation from fine to coarse
65  * Each element in set, or in interface if set is not input, is converted to coarse elements, with
66  * fine vertices put into SPECTRAL_VERTICES tagged array. NOTE: This function assumes that each
67  * order^d (fine) elements comprise each coarse element, and are in order of fine elements in each
68  * coarse element. If order is input as 0, looks for a SPECTRAL_ORDER tag on the mesh.
69  * \param order Order of the spectral mesh
70  * \param spectral_set Set containing spectral elements
71  */
73 {
74  if( order ) spectralOrder = order;
75  if( !spectralOrder )
76  {
77  MB_SET_ERR( MB_FAILURE, "Spectral order must be set or input before converting to spectral mesh" );
78  }
79 
80  Range tmp_ents, ents;
81  ErrorCode rval = mbImpl->get_entities_by_handle( spectral_set, tmp_ents );
82  if( MB_SUCCESS != rval || ents.empty() ) return rval;
83 
84  // get the max-dimensional elements from it
85  ents = tmp_ents.subset_by_dimension( 3 );
86  if( ents.empty() ) ents = tmp_ents.subset_by_dimension( 2 );
87  if( ents.empty() ) ents = tmp_ents.subset_by_dimension( 1 );
88  if( ents.empty() )
89  {
90  MB_SET_ERR( MB_FAILURE, "Can't find any entities for conversion" );
91  }
92 
93  // get a ptr to connectivity
94  if( ents.psize() != 1 )
95  {
96  MB_SET_ERR( MB_FAILURE, "Entities must be in one chunk for conversion" );
97  }
98  EntityHandle* conn;
99  int count, verts_per_e;
100  rval = mbImpl->connect_iterate( ents.begin(), ents.end(), conn, verts_per_e, count );
101  if( MB_SUCCESS != rval || count != (int)ents.size() ) return rval;
102 
103  Range tmp_range;
104  return create_spectral_elems( conn, ents.size(), CN::Dimension( TYPE_FROM_HANDLE( *ents.begin() ) ), tmp_range );
105 }
106 
107 template < class T >
109  int num_fine_elems,
110  int dim,
111  Range& output_range,
112  int start_idx,
113  Range* local_gids )
114 {
115  assert( spectralOrder && num_fine_elems );
116 
117  // now create num_coarse_elems
118  // compute the number of local elems, accounting for coarse or fine representation
119  // spectral_unit is the # fine elems per coarse elem, or spectralOrder^2
120  int spectral_unit = spectralOrder * spectralOrder;
121  int num_coarse_elems = num_fine_elems / spectral_unit;
122 
123  EntityHandle* new_conn;
124  EntityHandle start_elem;
125  ReadUtilIface* rmi;
126  ErrorCode rval = mbImpl->query_interface( rmi );
127  if( MB_SUCCESS != rval ) return rval;
128 
129  int verts_per_felem = spectralOrderp1 * spectralOrderp1, verts_per_celem = std::pow( (double)2.0, dim );
130 
131  rval = rmi->get_element_connect( num_coarse_elems, verts_per_celem, ( 2 == dim ? MBQUAD : MBHEX ), 0, start_elem,
132  new_conn );MB_CHK_SET_ERR( rval, "Failed to create elems" );
133 
134  output_range.insert( start_elem, start_elem + num_coarse_elems - 1 );
135 
136  // read connectivity into that space
137 
138  // permute_array takes a 4*order^2-long vector of integers, representing the connectivity of
139  // order^2 elems (fine elems in a coarse elem), and picks out the ids of the vertices necessary
140  // to get a lexicographically-ordered array of vertices in a spectral element of that order
141  // assert(verts_per_felem == (sizeof(permute_array)/sizeof(unsigned int)));
142 
143  // we're assuming here that elems was empty on input
144  int count;
145  EntityHandle* sv_ptr = NULL;
146  rval = mbImpl->tag_iterate( spectral_vertices_tag( true ), output_range.begin(), output_range.end(), count,
147  (void*&)sv_ptr );MB_CHK_SET_ERR( rval, "Failed to get SPECTRAL_VERTICES ptr" );
148  assert( count == num_coarse_elems );
149  int f = start_idx, fs = 0, fl = 0;
150  for( int c = 0; c < num_coarse_elems; c++ )
151  {
152  for( int i = 0; i < verts_per_celem; i++ )
153  new_conn[fl + i] = conn[f + lin_permute_array[i]];
154  fl += verts_per_celem;
155  for( int i = 0; i < verts_per_felem; i++ )
156  sv_ptr[fs + i] = conn[f + permute_array[i]];
157  f += verts_per_celem * spectral_unit;
158  fs += verts_per_felem;
159  }
160  if( local_gids ) std::copy( sv_ptr, sv_ptr + verts_per_felem * num_coarse_elems, range_inserter( *local_gids ) );
161 
162  return MB_SUCCESS;
163 }
164 
165 // force instantiation of a few specific types
166 template ErrorCode SpectralMeshTool::create_spectral_elems< int >( const int* conn,
167  int num_fine_elems,
168  int dim,
169  Range& output_range,
170  int start_idx,
171  Range* local_gids );
172 template ErrorCode SpectralMeshTool::create_spectral_elems< EntityHandle >( const EntityHandle* conn,
173  int num_fine_elems,
174  int dim,
175  Range& output_range,
176  int start_idx,
177  Range* local_gids );
178 } // namespace moab