1 #include "moab/SpectralMeshTool.hpp"
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 <assert.h>
9
10 namespace moab
11 {
12 const short int SpectralMeshTool::permute_array[] =
13 {0, 1, 13, 25, 3, 2, 14, 26, 7, 6, 18, 30, 11, 10, 22, 34};
14
15 // lin_permute_array does the same to get the linear vertices of the coarse quad
16 const short int SpectralMeshTool::lin_permute_array[] =
17 {0, 25, 34, 11};
18
spectral_vertices_tag(const bool create_if_missing)19 Tag SpectralMeshTool::spectral_vertices_tag(const bool create_if_missing)
20 {
21 ErrorCode rval = MB_SUCCESS;
22 if (!svTag && create_if_missing) {
23 if (!spectralOrder) {
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,
31 svTag, 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 */
spectral_order_tag(const bool create_if_missing)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 // create it
46 int dum = 0;
47 rval = mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER,
48 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 */
convert_to_fine(EntityHandle)59 ErrorCode SpectralMeshTool::convert_to_fine(EntityHandle /* spectral_set */)
60 {
61 return MB_NOT_IMPLEMENTED;
62 }
63
64
65 /** \brief Convert representation from fine to coarse
66 * Each element in set, or in interface if set is not input, is converted to coarse elements, with
67 * fine vertices put into SPECTRAL_VERTICES tagged array. NOTE: This function assumes that each
68 * order^d (fine) elements comprise each coarse element, and are in order of fine elements in each
69 * coarse element. If order is input as 0, looks for a SPECTRAL_ORDER tag on the mesh.
70 * \param order Order of the spectral mesh
71 * \param spectral_set Set containing spectral elements
72 */
convert_to_coarse(int order,EntityHandle spectral_set)73 ErrorCode SpectralMeshTool::convert_to_coarse(int order, EntityHandle spectral_set)
74 {
75 if (order) spectralOrder = order;
76 if (!spectralOrder) {
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 MB_SET_ERR(MB_FAILURE, "Can't find any entities for conversion");
90 }
91
92 // get a ptr to connectivity
93 if (ents.psize() != 1) {
94 MB_SET_ERR(MB_FAILURE, "Entities must be in one chunk for conversion");
95 }
96 EntityHandle *conn;
97 int count, verts_per_e;
98 rval = mbImpl->connect_iterate(ents.begin(), ents.end(), conn, verts_per_e, count);
99 if (MB_SUCCESS != rval || count != (int)ents.size()) return rval;
100
101 Range tmp_range;
102 return create_spectral_elems(conn, ents.size(), CN::Dimension(TYPE_FROM_HANDLE(*ents.begin())), tmp_range);
103 }
104
105 template <class T>
create_spectral_elems(const T * conn,int num_fine_elems,int dim,Range & output_range,int start_idx,Range * local_gids)106 ErrorCode SpectralMeshTool::create_spectral_elems(const T *conn, int num_fine_elems, int dim,
107 Range &output_range, int start_idx, Range *local_gids)
108 {
109 assert(spectralOrder && num_fine_elems);
110
111 // now create num_coarse_elems
112 // compute the number of local elems, accounting for coarse or fine representation
113 // spectral_unit is the # fine elems per coarse elem, or spectralOrder^2
114 int spectral_unit = spectralOrder*spectralOrder;
115 int num_coarse_elems = num_fine_elems / spectral_unit;
116
117 EntityHandle *new_conn;
118 EntityHandle start_elem;
119 ReadUtilIface *rmi;
120 ErrorCode rval = mbImpl->query_interface(rmi);
121 if (MB_SUCCESS != rval) return rval;
122
123 int verts_per_felem = spectralOrderp1*spectralOrderp1,
124 verts_per_celem = std::pow((double)2.0, dim);
125
126 rval = rmi->get_element_connect(num_coarse_elems, verts_per_celem,
127 (2 == dim ? MBQUAD : MBHEX), 0,
128 start_elem, new_conn);MB_CHK_SET_ERR(rval, "Failed to create elems");
129
130 output_range.insert(start_elem, start_elem + num_coarse_elems - 1);
131
132 // read connectivity into that space
133
134 // permute_array takes a 4*order^2-long vector of integers, representing the connectivity of order^2
135 // elems (fine elems in a coarse elem), and picks out the ids of the vertices necessary
136 // to get a lexicographically-ordered array of vertices in a spectral element of that order
137 //assert(verts_per_felem == (sizeof(permute_array)/sizeof(unsigned int)));
138
139 // we're assuming here that elems was empty on input
140 int count;
141 EntityHandle *sv_ptr = NULL;
142 rval = mbImpl->tag_iterate(spectral_vertices_tag(true), output_range.begin(), output_range.end(), count,
143 (void*&)sv_ptr);MB_CHK_SET_ERR(rval, "Failed to get SPECTRAL_VERTICES ptr");
144 assert(count == num_coarse_elems);
145 int f = start_idx, fs = 0, fl = 0;
146 for (int c = 0; c < num_coarse_elems; c++) {
147 for (int i = 0; i < verts_per_celem; i++)
148 new_conn[fl+i] = conn[f+lin_permute_array[i]];
149 fl += verts_per_celem;
150 for (int i = 0; i < verts_per_felem; i++)
151 sv_ptr[fs+i] = conn[f+permute_array[i]];
152 f += verts_per_celem*spectral_unit;
153 fs += verts_per_felem;
154 }
155 if (local_gids)
156 std::copy(sv_ptr, sv_ptr+verts_per_felem*num_coarse_elems, range_inserter(*local_gids));
157
158 return MB_SUCCESS;
159 }
160
161 // force instantiation of a few specific types
162 template ErrorCode SpectralMeshTool::create_spectral_elems<int>(const int *conn, int num_fine_elems, int dim,
163 Range &output_range, int start_idx, Range *local_gids);
164 template ErrorCode SpectralMeshTool::create_spectral_elems<EntityHandle>(const EntityHandle *conn, int num_fine_elems, int dim,
165 Range &output_range, int start_idx, Range *local_gids);
166 } // namespace moab
167
168