Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Mesh/Reorder.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Mesh reordering routines for unstructured meshes
9 : : \details Mesh reordering routines for unstructured meshes.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <algorithm>
14 : : #include <iterator>
15 : : #include <unordered_map>
16 : : #include <map>
17 : : #include <tuple>
18 : : #include <cstddef>
19 : :
20 : : #include "Reorder.hpp"
21 : : #include "Exception.hpp"
22 : : #include "ContainerUtil.hpp"
23 : : #include "Vector.hpp"
24 : :
25 : : namespace tk {
26 : :
27 : : std::size_t
28 : 91 : shiftToZero( std::vector< std::size_t >& inpoel )
29 : : // *****************************************************************************
30 : : // Shift node IDs to start with zero in element connectivity
31 : : //! \param[inout] inpoel Inteconnectivity of points and elements
32 : : //! \return Amount shifted
33 : : //! \details This function implements a simple reordering of the node ids of the
34 : : //! element connectivity in inpoel by shifting the node ids so that the
35 : : //! smallest is zero.
36 : : //! \note It is okay to call this function with an empty container; it will
37 : : //! simply return without throwing an exception.
38 : : // *****************************************************************************
39 : : {
40 [ + + ]: 91 : if (inpoel.empty()) return 0;
41 : :
42 : : // find smallest node id
43 : 55 : auto minId = *std::min_element( begin(inpoel), end(inpoel) );
44 : :
45 : : // shift node ids to start from zero
46 : : // cppcheck-suppress useStlAlgorithm
47 [ + + ]: 2107374 : for (auto& n : inpoel) n -= minId;
48 : :
49 : 55 : return minId;
50 : : }
51 : :
52 : : void
53 : 3624 : remap( std::vector< std::size_t >& ids, const std::vector< std::size_t >& map )
54 : : // *****************************************************************************
55 : : // Apply new maping to vector of indices
56 : : //! \param[inout] ids Vector of integer IDs to remap
57 : : //! \param[in] map Array of indices creating a new order
58 : : //! \details This function applies a mapping (reordering) to the integer IDs
59 : : //! passed in using the map passed in. The mapping is expressed between the
60 : : //! array index and its value. The function overwrites every value, i, of
61 : : //! vector ids with map[i].
62 : : //! \note The sizes of ids and map need not equal. Only the maximum index in ids
63 : : //! must be lower than the size of map.
64 : : //! \note It is okay to call this function with either of the containers empty;
65 : : //! it will simply return without throwing an exception.
66 : : // *****************************************************************************
67 : : {
68 [ + - ][ - + ]: 3624 : if (ids.empty() || map.empty()) return;
[ - + ]
69 : :
70 [ - + ][ - - ]: 3624 : Assert( *max_element( begin(ids), end(ids) ) < map.size(),
[ - - ][ - - ]
71 : : "Indexing out of bounds" );
72 : :
73 : : // remap integer IDs in vector ids
74 : : // cppcheck-suppress useStlAlgorithm
75 [ + + ]: 7033088 : for (auto& i : ids) i = map[i];
76 : : }
77 : :
78 : : void
79 : 3 : remap( std::vector< tk::real >& r, const std::vector< std::size_t >& map )
80 : : // *****************************************************************************
81 : : // Apply new maping to vector of real numbers
82 : : //! \param[inout] r Vector of real numbers to remap
83 : : //! \param[in] map Array of indices creating a new order
84 : : //! \details This function applies a mapping (reordering) to the real values
85 : : //! passed in using the map passed in. The mapping is expressed between the
86 : : //! array index and its value. The function moves every value r[i] to
87 : : //! r[ map[i] ].
88 : : //! \note The sizes of r and map must be equal and the maximum index in map must
89 : : //! be lower than the size of map.
90 : : //! \note It is okay to call this function with either of the containers empty;
91 : : //! it will simply return without throwing an exception.
92 : : // *****************************************************************************
93 : : {
94 [ + - ][ - + ]: 3 : if (r.empty() || map.empty()) return;
[ - + ]
95 : :
96 [ - + ][ - - ]: 3 : Assert( r.size() == map.size(), "Size mismatch" );
[ - - ][ - - ]
97 [ + - ][ - + ]: 3 : Assert( *max_element( begin(map), end(map) ) < map.size(),
[ - - ][ - - ]
[ - - ]
98 : : "Indexing out of bounds" );
99 : :
100 : : // remap real numbers in vector
101 [ + - ]: 6 : auto m = r;
102 [ + + ]: 52731 : for (std::size_t i=0; i<map.size(); ++i) r[ map[i] ] = m[ i ];
103 : : }
104 : :
105 : : std::vector< std::size_t >
106 : 2872 : remap( const std::vector< std::size_t >& ids,
107 : : const std::vector< std::size_t >& map )
108 : : // *****************************************************************************
109 : : // Create remapped vector of indices using a vector
110 : : //! \param[in] ids Vector of integer IDs to remap
111 : : //! \param[in] map Array of indices creating a new order
112 : : //! \return Remapped vector of ids
113 : : //! \details This function applies a mapping (reordering) to the integer IDs
114 : : //! passed in using the map passed in. The mapping is expressed between the
115 : : //! array index and its value. The function creates and returns a new container
116 : : //! with remapped ids of identical size of the origin ids container.
117 : : //! \note The sizes of ids and map must be equal and the maximum index in map
118 : : //! must be lower than the size of map.
119 : : //! \note It is okay to call this function with either of the containers empty;
120 : : //! if ids is empty, it returns an empty container; if map is empty, it will
121 : : //! return the original container.
122 : : // *****************************************************************************
123 : : {
124 [ + + ]: 2872 : if (ids.empty()) return {};
125 [ - + ][ - - ]: 2856 : if (map.empty()) return ids;
126 : :
127 [ + - ][ - + ]: 2856 : Assert( *max_element( begin(ids), end(ids) ) < map.size(),
[ - - ][ - - ]
[ - - ]
128 : : "Indexing out of bounds" );
129 : :
130 : : // in terms of the in-place remap of a vector usinga vector
131 [ + - ]: 5712 : auto newids = ids;
132 [ + - ]: 2856 : remap( newids, map );
133 : :
134 : 2856 : return newids;
135 : : }
136 : :
137 : : void
138 : 13304 : remap( std::vector< std::size_t >& ids,
139 : : const std::unordered_map< std::size_t, std::size_t >& map )
140 : : // *****************************************************************************
141 : : // In-place remap vector of indices using a map
142 : : //! \param[in] ids Vector of integer IDs to remap
143 : : //! \param[in] map Hash-map of key->value creating a new order
144 : : //! \details This function applies a mapping (reordering) to the integer IDs
145 : : //! passed in using the map passed in. The mapping is expressed as a hash-map
146 : : //! of key->value pairs, where the key is the original and the value is the
147 : : //! new ids of the mapping. The function overwrites the ids container with the
148 : : //! remapped ids of identical size.
149 : : //! \note All ids in the input ids container must have a key in the map.
150 : : //! Otherwise an exception is thrown.
151 : : //! \note It is okay to call this function with the ids container empty but not
152 : : //! okay to pass an empty map.
153 : : // *****************************************************************************
154 : : {
155 [ - + ][ - - ]: 13304 : Assert( !map.empty(), "Map must not be empty" );
[ - - ][ - - ]
156 : :
157 [ + + ][ + - ]: 6899927 : for (auto& i : ids) i = tk::cref_find( map, i );
158 : 13304 : }
159 : :
160 : : std::vector< std::size_t >
161 : 5165 : remap( const std::vector< std::size_t >& ids,
162 : : const std::unordered_map< std::size_t, std::size_t >& map )
163 : : // *****************************************************************************
164 : : // Create remapped vector of indices using a map
165 : : //! \param[in] ids Vector of integer IDs to create new container of ids from
166 : : //! \param[in] map Hash-map of key->value creating a new order
167 : : //! \return Remapped vector of ids
168 : : //! \details This function applies a mapping (reordering) to the integer IDs
169 : : //! passed in using the map passed in. The mapping is expressed as a hash-map
170 : : //! of key->value pairs, where the key is the original and the value is the
171 : : //! new ids of the mapping. The function creates and returns a new container
172 : : //! with the remapped ids of identical size of the original ids container.
173 : : //! \note All ids in the input ids container must have a key in the map.
174 : : //! Otherwise an exception is thrown.
175 : : //! \note It is okay to call this function with the ids container empty but not
176 : : //! okay to pass an empty map.
177 : : // *****************************************************************************
178 : : {
179 [ - + ][ - - ]: 5165 : Assert( !map.empty(), "Map must not be empty" );
[ - - ][ - - ]
180 : :
181 : : // in terms of the in-place remap of a vector using a map
182 : 5165 : auto newids = ids;
183 [ + - ]: 5165 : remap( newids, map );
184 : :
185 : 5165 : return newids;
186 : : }
187 : :
188 : : std::map< int, std::vector< std::size_t > >
189 : 3806 : remap( const std::map< int, std::vector< std::size_t > >& ids,
190 : : const std::unordered_map< std::size_t, std::size_t >& map )
191 : : // *****************************************************************************
192 : : // Create remapped map of vector of indices using a map
193 : : //! \param[in] ids Map of vector of integer IDs to create new container of ids
194 : : //! from
195 : : //! \param[in] map Hash-map of key->value creating a new order
196 : : //! \return Remapped vector of ids
197 : : //! \details This function applies a mapping (reordering) to the map of integer
198 : : //! IDs passed in using the map passed in by applying remap(vector,map) on
199 : : //! each vector of ids. The keys in the returned map will be the same as in
200 : : //! ids.
201 : : // *****************************************************************************
202 : : {
203 [ - + ][ - - ]: 3806 : Assert( !map.empty(), "Map must not be empty" );
[ - - ][ - - ]
204 : :
205 : : // in terms of the in-place remap of a vector using a map
206 : 3806 : auto newids = ids;
207 [ + + ][ + - ]: 11429 : for (auto& m : newids) remap( m.second, map );
208 : :
209 : 3806 : return newids;
210 : : }
211 : :
212 : : std::vector< std::size_t >
213 : 1 : renumber( const std::pair< std::vector< std::size_t >,
214 : : std::vector< std::size_t > >& psup )
215 : : // *****************************************************************************
216 : : // Reorder mesh points with the advancing front technique
217 : : //! \param[in] psup Points surrounding points
218 : : //! \return Mapping created by renumbering (reordering)
219 : : // *****************************************************************************
220 : : {
221 : : // Find out number of nodes in graph
222 : 1 : auto npoin = psup.second.size()-1;
223 : :
224 : : // Construct mapping using advancing front
225 [ + - ][ + - ]: 2 : std::vector< int > hpoin( npoin, -1 ), lpoin( npoin, 0 );
226 [ + - ]: 1 : std::vector< std::size_t > map( npoin, 0 );
227 : 1 : hpoin[0] = 0;
228 : 1 : lpoin[0] = 1;
229 : 1 : std::size_t num = 1;
230 [ + + ]: 26 : while (num < npoin) {
231 : 25 : std::size_t cnt = 0;
232 : 25 : std::size_t i = 0;
233 [ + - ]: 50 : std::vector< int > kpoin( npoin, -1 );
234 : : int p;
235 [ + + ]: 15604 : while ((p = hpoin[i]) != -1) {
236 : 15579 : ++i;
237 : 15579 : auto P = static_cast< std::size_t >( p );
238 [ + + ]: 226609 : for (auto j=psup.second[P]+1; j<=psup.second[P+1]; ++j) {
239 : 211030 : auto q = psup.first[j];
240 [ + + ]: 211030 : if (lpoin[q] != 1) { // consider points not yet counted
241 : 17575 : map[q] = num++;
242 : 17575 : kpoin[cnt] = static_cast< int >( q ); // register point as counted
243 : 17575 : lpoin[q] = 1; // register the point as counted
244 : 17575 : ++cnt;
245 : : }
246 : : }
247 : : }
248 [ + - ]: 25 : hpoin = kpoin;
249 : : }
250 : :
251 : : // // Construct new->old id map
252 : : // std::size_t i = 0;
253 : : // std::vector< std::size_t > oldmap( npoin );
254 : : // for (auto n : map) oldmap[n] = i++;
255 : :
256 : : // Return old->new and new->old maps
257 : 2 : return map;
258 : : }
259 : :
260 : : std::unordered_map< std::size_t, std::size_t >
261 : 8095 : assignLid( const std::vector< std::size_t >& gid )
262 : : // *****************************************************************************
263 : : // Assign local ids to global ids
264 : : //! \param[in] gid Global ids
265 : : //! \return Map associating global ids to local ids
266 : : // *****************************************************************************
267 : : {
268 : 8095 : std::unordered_map< std::size_t, std::size_t > lid;
269 : 8095 : std::size_t l = 0;
270 [ + + ][ + - ]: 3498612 : for (auto p : gid) lid[p] = l++;
271 : 8095 : return lid;
272 : : }
273 : :
274 : : std::tuple< std::vector< std::size_t >,
275 : : std::vector< std::size_t >,
276 : : std::unordered_map< std::size_t, std::size_t > >
277 : 4963 : global2local( const std::vector< std::size_t >& ginpoel )
278 : : // *****************************************************************************
279 : : // Generate element connectivity of local node IDs from connectivity of global
280 : : // node IDs also returning the mapping between local to global IDs
281 : : //! \param[in] ginpoel Element connectivity with global node IDs
282 : : //! \return Tuple of (1) element connectivity with local node IDs, (2) the
283 : : //! vector of unique global node IDs (i.e., the mapping between local to
284 : : //! global node IDs), and (3) mapping between global to local node IDs.
285 : : // *****************************************************************************
286 : : {
287 : : // Make a copy of the element connectivity with global node ids
288 [ + - ]: 9926 : auto gid = ginpoel;
289 : :
290 : : // Generate a vector that holds only the unique global mesh node ids
291 [ + - ]: 4963 : tk::unique( gid );
292 : :
293 : : // Assign local node ids to global node ids
294 [ + - ]: 9926 : const auto lid = tk::assignLid( gid );
295 : :
296 [ - + ][ - - ]: 4963 : Assert( gid.size() == lid.size(), "Size mismatch" );
[ - - ][ - - ]
297 : :
298 : : // Generate element connectivity using local node ids
299 [ + - ]: 9926 : std::vector< std::size_t > inpoel( ginpoel.size() );
300 : 4963 : std::size_t j = 0;
301 [ + + ][ + - ]: 28529507 : for (auto p : ginpoel) inpoel[ j++ ] = tk::cref_find( lid, p );
302 : :
303 : : // Return element connectivty with local node IDs
304 [ + - ]: 9926 : return std::make_tuple( inpoel, gid, lid );
305 : : }
306 : :
307 : : bool
308 : 4595 : positiveJacobians( const std::vector< std::size_t >& inpoel,
309 : : const std::array< std::vector< real >, 3 >& coord )
310 : : // *****************************************************************************
311 : : // Test for positivity of the Jacobian for all cells in mesh
312 : : //! \param[in] inpoel Element connectivity (zero-based, i.e., local if parallel)
313 : : //! \param[in] coord Node coordinates
314 : : //! \return True if Jacobians of all mesh cells are positive
315 : : // *****************************************************************************
316 : : {
317 [ - + ][ - - ]: 4595 : Assert( !inpoel.empty(), "Mesh connectivity empty" );
[ - - ][ - - ]
318 [ - + ][ - - ]: 4595 : Assert( inpoel.size() % 4 == 0,
[ - - ][ - - ]
319 : : "Mesh connectivity size must be divisible by 4 " );
320 [ - + ][ - - ]: 4595 : Assert( tk::uniquecopy(inpoel).size() == coord[0].size(), "Number of unique "
[ - - ][ - - ]
321 : : "nodes in mesh connectivity must equal the number of nodes to which "
322 : : "coordinates have been supplied" );
323 [ - + ][ - - ]: 4595 : Assert( tk::uniquecopy(inpoel).size() == coord[1].size(), "Number of unique "
[ - - ][ - - ]
324 : : "nodes in mesh connectivity must equal the number of nodes to which "
325 : : "coordinates have been supplied" );
326 [ - + ][ - - ]: 4595 : Assert( tk::uniquecopy(inpoel).size() == coord[2].size(), "Number of unique "
[ - - ][ - - ]
327 : : "nodes in mesh connectivity must equal the number of nodes to which "
328 : : "coordinates have been supplied" );
329 [ - + ][ - - ]: 4595 : Assert( *std::minmax_element( begin(inpoel), end(inpoel) ).first == 0,
[ - - ][ - - ]
330 : : "node ids should start from zero" );
331 : :
332 : 4595 : const auto& x = coord[0];
333 : 4595 : const auto& y = coord[1];
334 : 4595 : const auto& z = coord[2];
335 : :
336 [ + + ]: 5604566 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
337 : 5599971 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
338 : 5599971 : inpoel[e*4+2], inpoel[e*4+3] }};
339 : : // compute element Jacobi determinant / (5/120) = element volume * 4
340 : : const std::array< tk::real, 3 >
341 : 5599971 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
342 : 5599971 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
343 : 5599971 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
344 [ - + ]: 5599971 : if (tk::triple( ba, ca, da ) < 0) return false;
345 : : }
346 : :
347 : 4595 : return true;
348 : : }
349 : :
350 : : std::map< int, std::vector< std::size_t > >
351 : 3619 : bfacenodes( const std::map< int, std::vector< std::size_t > >& bface,
352 : : const std::vector< std::size_t >& triinpoel )
353 : : // *****************************************************************************
354 : : // Generate nodes of side set faces
355 : : //! \param[in] bface Boundary-faces mapped to side set ids
356 : : //! \param[in] triinpoel Boundary-face connectivity
357 : : //! \return Nodes of side set faces for each side set passed in
358 : : // *****************************************************************************
359 : : {
360 : 3619 : auto bfn = bface;
361 : :
362 [ + + ]: 11446 : for (auto& [s,b] : bfn) {
363 : 15654 : std::vector< std::size_t > nodes;
364 [ + + ]: 385237 : for (auto f : b) {
365 [ + - ]: 377410 : nodes.push_back( triinpoel[f*3+0] );
366 [ + - ]: 377410 : nodes.push_back( triinpoel[f*3+1] );
367 [ + - ]: 377410 : nodes.push_back( triinpoel[f*3+2] );
368 : : }
369 [ + - ]: 7827 : tk::unique( nodes );
370 : 7827 : b = std::move( nodes );
371 : : }
372 : :
373 : 3619 : return bfn;
374 : : }
375 : :
376 : : } // tk::
|
