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