1 // --------------------------------------------------------------------- 2 // 3 // Copyright (C) 2005 - 2020 by the deal.II authors 4 // 5 // This file is part of the deal.II library. 6 // 7 // The deal.II library is free software; you can use it, redistribute 8 // it, and/or modify it under the terms of the GNU Lesser General 9 // Public License as published by the Free Software Foundation; either 10 // version 2.1 of the License, or (at your option) any later version. 11 // The full text of the license can be found in the file LICENSE.md at 12 // the top level directory of deal.II. 13 // 14 // --------------------------------------------------------------------- 15 16 #ifndef dealii_utilities_h 17 #define dealii_utilities_h 18 19 #include <deal.II/base/config.h> 20 21 #include <deal.II/base/exceptions.h> 22 23 #include <functional> 24 #include <string> 25 #include <tuple> 26 #include <type_traits> 27 #include <typeinfo> 28 #include <utility> 29 #include <vector> 30 31 #ifdef DEAL_II_WITH_TRILINOS 32 # include <Epetra_Comm.h> 33 # include <Epetra_Map.h> 34 # include <Teuchos_Comm.hpp> 35 # include <Teuchos_RCP.hpp> 36 # ifdef DEAL_II_WITH_MPI 37 # include <Epetra_MpiComm.h> 38 # else 39 # include <Epetra_SerialComm.h> 40 # endif 41 #endif 42 43 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS 44 45 #include <boost/archive/binary_iarchive.hpp> 46 #include <boost/archive/binary_oarchive.hpp> 47 #include <boost/core/demangle.hpp> 48 #include <boost/serialization/array.hpp> 49 #include <boost/serialization/complex.hpp> 50 #include <boost/serialization/vector.hpp> 51 52 #ifdef DEAL_II_WITH_ZLIB 53 # include <boost/iostreams/device/back_inserter.hpp> 54 # include <boost/iostreams/filter/gzip.hpp> 55 # include <boost/iostreams/filtering_stream.hpp> 56 # include <boost/iostreams/stream.hpp> 57 #endif 58 59 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS 60 61 DEAL_II_NAMESPACE_OPEN 62 63 // forward declare Point 64 #ifndef DOXYGEN 65 template <int dim, typename Number> 66 class Point; 67 #endif 68 69 /** 70 * A namespace for utility functions that are not particularly specific to 71 * finite element computing or numerical programs, but nevertheless are needed 72 * in various contexts when writing applications. 73 * 74 * @ingroup utilities 75 */ 76 namespace Utilities 77 { 78 /** 79 * Return a string of the form "deal.II version x.y.z" where "x.y.z" 80 * identifies the version of deal.II you are using. This information 81 * is also provided by the DEAL_II_PACKAGE_NAME and 82 * DEAL_II_PACKAGE_VERSION preprocessor variables. 83 */ 84 std::string 85 dealii_version_string(); 86 87 /** 88 * Assign to each point in @p points an index using the Hilbert space filling curve. 89 * To that end, a bounding box for @p points will be determined, based on which their 90 * integer coordinates are calculated. 91 * The linear index is given as a dim-collection of bits, from high to low. 92 * This is done in order to keep the maximum resolution in terms of bit depth 93 * along each axis. Note that this dim-integer index can still be easily used 94 * for sorting and ordering, for example using the lexicographic ordering of 95 * tuples of integers. 96 * 97 * The depth of the Hilbert curve (i.e. the number of bits per dimension) by 98 * default is equal to <code>64</code>. 99 * 100 * @note This function can also handle degenerate cases, e.g. when the bounding 101 * box has zero size in one of the dimensions. 102 */ 103 template <int dim, typename Number> 104 std::vector<std::array<std::uint64_t, dim>> 105 inverse_Hilbert_space_filling_curve( 106 const std::vector<Point<dim, Number>> &points, 107 const int bits_per_dim = 64); 108 109 /** 110 * Same as above, but for points in integer coordinates. 111 */ 112 template <int dim> 113 std::vector<std::array<std::uint64_t, dim>> 114 inverse_Hilbert_space_filling_curve( 115 const std::vector<std::array<std::uint64_t, dim>> &points, 116 const int bits_per_dim = 64); 117 118 /** 119 * Pack the least significant @p bits_per_dim bits from each element of @p index 120 * (starting from last) into a single unsigned integer. The last element 121 * of @p index will be used to set the first @p bits_per_dim bits in the 122 * resulting integer, the second to last element is used to set the next @p bits_per_dim bits, 123 * etc.. To fit all the data into the output, the following should hold 124 * <code>bits\_per\_dim * dim <= 64</code>. 125 * 126 * The function is useful in debugging and visualization of indices returned 127 * by inverse_Hilbert_space_filling_curve(). 128 * 129 * @note There is no need to use this function in order to compare indices 130 * returned by inverse_Hilbert_space_filling_curve(), as that can easily be 131 * done via <code>std::lexicographical_compare()</code>. 132 */ 133 template <int dim> 134 std::uint64_t 135 pack_integers(const std::array<std::uint64_t, dim> &index, 136 const int bits_per_dim); 137 138 /** 139 * If the library is configured with ZLIB, then this function compresses the 140 * input string and returns a non-zero terminated string containing the 141 * compressed input. 142 * 143 * If the library was not configured with ZLIB enabled, the returned string 144 * is identical to the input string. 145 * 146 * @param[in] input The string to compress 147 * 148 * @return A compressed version of the input string 149 */ 150 std::string 151 compress(const std::string &input); 152 153 /** 154 * If the library is configured with ZLIB, then this function assumes that the 155 * input string has been compressed using the compress() function, and returns 156 * the original decompressed string. 157 * 158 * If the library was not configured with ZLIB enabled, the returned string 159 * is identical to the input string. 160 * 161 * @param[in] compressed_input A compressed string, as returned by the 162 * function compress() 163 * 164 * @return The original uncompressed string. 165 */ 166 std::string 167 decompress(const std::string &compressed_input); 168 169 /** 170 * Encodes the binary input as a base64 string. 171 * 172 * Base64 is a group of binary-to-text encoding schemes that represent binary 173 * data in an ASCII string format by translating it into a radix-64 174 * representation. Base64 is designed to carry data stored in binary formats 175 * across channels that only reliably support text content. It is used also 176 * to store binary formats in a machine independent way. 177 * 178 * @param binary_input A vector of characters, representing your input as 179 * binary data. 180 * @return A string containing the binary input as a base64 string. 181 */ 182 std::string 183 encode_base64(const std::vector<unsigned char> &binary_input); 184 185 /** 186 * Decodes a base64 string into a binary output. 187 * 188 * This is the inverse of the encode_base64() function above. 189 * 190 * @param base64_input A string that contains the input in base64 format. 191 * @return A vector of characters that represents your input as binary data. 192 */ 193 std::vector<unsigned char> 194 decode_base64(const std::string &base64_input); 195 196 /** 197 * Convert a number @p value to a string, with as many digits as given to 198 * fill with leading zeros. 199 * 200 * If the second parameter is left at its default value, the number is not 201 * padded with leading zeros. The result is then the same as if the standard 202 * C++ `std::to_string` (or the older C function `itoa()`) had been called. 203 * 204 * This function takes an `unsigned int` as argument. As a consequence, 205 * if you call it with a `signed int` (which is of course the same 206 * type as `int`), the argument is implicitly converted to 207 * unsigned integers and negative numbers may not be printed as you had 208 * hoped. Similarly, if you call the function with a `long int`, the 209 * printed result might show the effects of an overflow upon conversion 210 * to `unsigned int`. 211 * 212 * @note The use of this function is discouraged and users should use 213 * <code>Utilities::to_string()</code> instead. In its current 214 * implementation the function simply calls <code>to_string@<unsigned 215 * int@>()</code>. 216 */ 217 std::string 218 int_to_string(const unsigned int value, 219 const unsigned int digits = numbers::invalid_unsigned_int); 220 221 /** 222 * Convert a number @p value to a string, with @p digits characters. The 223 * string is padded with leading zeros, after a possible minus sign. 224 * Therefore the total number of padding zeros is @p digits minus any signs, 225 * decimal points and digits of @p value. 226 * 227 * If the second parameter is left at its default value, the number is not 228 * padded with leading zeros. The result is then the same as if the C++ 229 * function `std::to_string()` had been called (for integral types), 230 * or if `boost::lexical_cast()` had been called (for all other types). 231 */ 232 template <typename number> 233 std::string 234 to_string(const number value, 235 const unsigned int digits = numbers::invalid_unsigned_int); 236 237 /** 238 * Determine how many digits are needed to represent numbers at most as 239 * large as the given number. 240 */ 241 unsigned int 242 needed_digits(const unsigned int max_number); 243 244 /** 245 * This function allows to cut off a floating point number @p number 246 * after @p n_digits of accuracy, i.e., after @p n_digits decimal places 247 * in scientific floating point notation. When interpreted as rounding 248 * operation, this function reduces the absolute value of a floating point 249 * number and always rounds towards zero, since decimal places are simply 250 * cut off. 251 */ 252 template <typename Number> 253 Number 254 truncate_to_n_digits(const Number number, const unsigned int n_digits); 255 256 /** 257 * Given a string, convert it to an integer. Throw an assertion if that is 258 * not possible. 259 */ 260 int 261 string_to_int(const std::string &s); 262 263 /** 264 * Return a string describing the dimensions of the object. Often, functions 265 * in the deal.II library as well as in user codes need to define a string 266 * containing the template dimensions of some objects defined using two 267 * template parameters: dim (the topological dimension of the object) and 268 * spacedim (the dimension of the embedding Euclidean space). Since in all 269 * deal.II classes, by default spacedim is equal to dimension, the above 270 * string is usually contracted to "<dim>", instead of "<dim,spacedim>". 271 * This function returns a string containing "dim" if dim is equal to 272 * spacedim, otherwise it returns "dim,spacedim". 273 */ 274 std::string 275 dim_string(const int dim, const int spacedim); 276 277 /** 278 * Given a list of strings, convert it to a list of integers. Throw an 279 * assertion if that is not possible. 280 */ 281 std::vector<int> 282 string_to_int(const std::vector<std::string> &s); 283 284 /** 285 * Given a string, convert it to an double. Throw an assertion if that is 286 * not possible. 287 */ 288 double 289 string_to_double(const std::string &s); 290 291 292 /** 293 * Given a list of strings, convert it to a list of doubles. Throw an 294 * assertion if that is not possible. 295 */ 296 std::vector<double> 297 string_to_double(const std::vector<std::string> &s); 298 299 300 /** 301 * Given a string that contains text separated by a @p delimiter, split it 302 * into its components; for each component, remove leading and trailing 303 * spaces. The default value of the delimiter is a comma, so that the 304 * function splits comma separated lists of strings. 305 * 306 * To make data input from tables simpler, if the input string ends in a 307 * delimiter (possibly followed by an arbitrary amount of whitespace), then 308 * this last delimiter is ignored. For example, 309 * @code 310 * Utilities::split_string_list("abc; def; ghi; ", ';'); 311 * @endcode 312 * yields the same 3-element list of output <code>{"abc","def","ghi"}</code> 313 * as you would get if the input had been 314 * @code 315 * Utilities::split_string_list("abc; def; ghi", ';'); 316 * @endcode 317 * or 318 * @code 319 * Utilities::split_string_list("abc; def; ghi;", ';'); 320 * @endcode 321 * As a consequence of this rule, a call like 322 * @code 323 * Utilities::split_string_list(" ; ", ';'); 324 * @endcode 325 * yields a one-element list. Because of the trimming of whitespace, the 326 * single element is the empty string. 327 * 328 * This function can digest the case that the delimiter is a space. In this 329 * case, it returns all words in the string. Combined with the rules above, 330 * this implies that 331 * @code 332 * Utilities::split_string_list("abc def ghi ", ' '); 333 * @endcode 334 * yields again the 3-element list of output 335 * <code>{"abc","def","ghi"}</code> from above despite the presence of space 336 * at the end of the string. Furthermore, 337 * @code 338 * Utilities::split_string_list(" ", ' '); 339 * @endcode 340 * yields an empty list regardless of the number of spaces in the string. 341 */ 342 std::vector<std::string> 343 split_string_list(const std::string &s, const std::string &delimiter = ","); 344 345 346 /** 347 * Specialization of split_string_list() for the case where the delimiter 348 * is a single char. 349 */ 350 std::vector<std::string> 351 split_string_list(const std::string &s, const char delimiter); 352 353 354 /** 355 * Take a text, usually a documentation or something, and try to break it 356 * into individual lines of text at most @p width characters wide, by 357 * breaking at positions marked by @p delimiter in the text. If this is not 358 * possible, return the shortest lines that are longer than @p width. The 359 * default value of the delimiter is a space character. If original_text 360 * contains newline characters (\n), the string is split at these locations, 361 * too. 362 */ 363 std::vector<std::string> 364 break_text_into_lines(const std::string &original_text, 365 const unsigned int width, 366 const char delimiter = ' '); 367 368 /** 369 * Return true if the given pattern string appears in the first position of 370 * the string. 371 */ 372 bool 373 match_at_string_start(const std::string &name, const std::string &pattern); 374 375 /** 376 * Read a (signed) integer starting at the position in @p name indicated by 377 * the second argument, and return this integer as a pair together with how 378 * many characters it takes up in the string. 379 * 380 * If no integer can be read at the indicated position, return 381 * (-1,numbers::invalid_unsigned_int) 382 */ 383 std::pair<int, unsigned int> 384 get_integer_at_position(const std::string &name, const unsigned int position); 385 386 /** 387 * Return a string with all occurrences of @p from in @p input replaced by 388 * @p to. 389 */ 390 std::string 391 replace_in_string(const std::string &input, 392 const std::string &from, 393 const std::string &to); 394 395 /** 396 * Return a string with all standard whitespace characters (including 397 * '<tt>\\t</tt>', '<tt>\\n</tt>', and '<tt>\\r</tt>') at the beginning and 398 * end of @p input removed. 399 */ 400 std::string 401 trim(const std::string &input); 402 403 /** 404 * Generate a random number from a normalized Gaussian probability 405 * distribution centered around @p a and with standard deviation @p sigma. 406 * The returned number will be different every time the function is called. 407 * 408 * This function is reentrant, i.e., it can safely be called from multiple 409 * threads at the same time. In addition, each thread will get the same 410 * sequence of numbers every time. On the other hand, if you run 411 * Threads::Task objects via the Threading Building Blocks, then tasks will 412 * be assigned to mostly random threads, and may get a different sequence of 413 * random numbers in different runs of the program, since a previous task 414 * may already have consumed the first few random numbers generated for the 415 * thread you're on. If this is a problem, you need to create your own 416 * random number generator objects every time you want to start from a 417 * defined point. 418 * 419 * @note Like the system function rand(), this function produces the same 420 * sequence of random numbers every time a program is started. This is an 421 * important property for debugging codes, but it makes it impossible to 422 * really verify statistical properties of a code. For `rand()`, you can call 423 * `srand()` to "seed" the random number generator to get different sequences 424 * of random numbers every time a program is called. However, this function 425 * does not allow seeding the random number generator. If you need this, as 426 * above, use one of the C++ or BOOST facilities. 427 */ 428 double 429 generate_normal_random_number(const double a, const double sigma); 430 431 /** 432 * Return a string description of the type of the variable @p t. 433 * 434 * In general, C++ uses mangled names to identify types. This function 435 * uses boost::core::demangle to return a human readable string describing 436 * the type of the variable passed as argument. 437 */ 438 template <class T> 439 std::string 440 type_to_string(const T &t); 441 442 /** 443 * Calculate a fixed power, provided as a template argument, of a number. 444 * 445 * This function provides an efficient way to calculate things like 446 * <code>t^N</code> where <code>N</code> is a known number at compile time. 447 * 448 * Use this function as in <code>fixed_power@<dim@> (n)</code>. 449 */ 450 template <int N, typename T> 451 T 452 fixed_power(const T t); 453 454 /** 455 * A replacement for <code>std::pow</code> that allows compile-time 456 * calculations for constant expression arguments. The @p base must 457 * be an integer type and the exponent @p iexp must not be negative. 458 */ 459 template <typename T> 460 constexpr T pow(const T base,const int iexp)461 pow(const T base, const int iexp) 462 { 463 #if defined(DEBUG) && !defined(DEAL_II_CXX14_CONSTEXPR_BUG) 464 // Up to __builtin_expect this is the same code as in the 'Assert' macro. 465 // The call to __builtin_expect turns out to be problematic. 466 if (!(iexp >= 0)) 467 ::dealii::deal_II_exceptions::internals::issue_error_noreturn( 468 ::dealii::deal_II_exceptions::internals::abort_or_throw_on_exception, 469 __FILE__, 470 __LINE__, 471 __PRETTY_FUNCTION__, 472 "iexp>=0", 473 "ExcMessage(\"The exponent must not be negative!\")", 474 ExcMessage("The exponent must not be negative!")); 475 #endif 476 // The "exponentiation by squaring" algorithm used below has to be 477 // compressed to one statement due to C++11's restrictions on constexpr 478 // functions. A more descriptive version would be: 479 // 480 // <code> 481 // if (iexp <= 0) 482 // return 1; 483 // 484 // // avoid overflow of one additional recursion with pow(base * base, 0) 485 // if (iexp == 1) 486 // return base; 487 // 488 // // if the current exponent is not divisible by two, 489 // // we need to account for that. 490 // const unsigned int prefactor = (iexp % 2 == 1) ? base : 1; 491 // 492 // // a^b = (a*a)^(b/2) for b even 493 // // a^b = a*(a*a)^((b-1)/2 for b odd 494 // return prefactor * dealii::Utilities::pow(base*base, iexp/2); 495 // </code> 496 497 static_assert(std::is_integral<T>::value, "Only integral types supported"); 498 499 return iexp <= 0 ? 500 1 : 501 (iexp == 1 ? base : 502 (((iexp % 2 == 1) ? base : 1) * 503 dealii::Utilities::pow(base * base, iexp / 2))); 504 } 505 506 /** 507 * Optimized replacement for <tt>std::lower_bound</tt> for searching within 508 * the range of column indices. Slashes execution time by approximately one 509 * half for the present application, partly because the binary search is 510 * replaced by a linear search for small loop lengths. 511 * 512 * Another reason for this function is rather obscure: when using the GCC 513 * libstdc++ function std::lower_bound, complexity is O(log(N)) as required. 514 * However, when using the debug version of the GCC libstdc++ as we do when 515 * running the testsuite, then std::lower_bound tests whether the sequence 516 * is in fact partitioned with respect to the pivot 'value' (i.e. in essence 517 * that the sequence is sorted as required for binary search to work). 518 * However, verifying this means that the complexity of std::lower_bound 519 * jumps to O(N); we call this function O(N) times below, making the overall 520 * complexity O(N**2). The consequence is that a few tests with big meshes 521 * completely run off the wall time limit for tests and fail with the 522 * libstdc++ debug mode 523 * 524 * This function simply makes the assumption that the sequence is sorted, 525 * and we simply don't do the additional check. 526 */ 527 template <typename Iterator, typename T> 528 Iterator 529 lower_bound(Iterator first, Iterator last, const T &val); 530 531 532 /** 533 * The same function as above, but taking an argument that is used to 534 * compare individual elements of the sequence of objects pointed to by the 535 * iterators. 536 */ 537 template <typename Iterator, typename T, typename Comp> 538 Iterator 539 lower_bound(Iterator first, Iterator last, const T &val, const Comp comp); 540 541 /** 542 * Given a permutation vector (i.e. a vector $p_0\ldots p_{N-1}$ where each 543 * $p_i\in [0,N)$ and $p_i\neq p_j$ for $i\neq j$), produce the reverse 544 * permutation $q_i=N-1-p_i$. 545 */ 546 template <typename Integer> 547 std::vector<Integer> 548 reverse_permutation(const std::vector<Integer> &permutation); 549 550 /** 551 * Given a permutation vector (i.e. a vector $p_0\ldots p_{N-1}$ where each 552 * $p_i\in [0,N)$ and $p_i\neq p_j$ for $i\neq j$), produce the inverse 553 * permutation $q_0\ldots q_{N-1}$ so that $q_{p_i}=p_{q_i}=i$. 554 */ 555 template <typename Integer> 556 std::vector<Integer> 557 invert_permutation(const std::vector<Integer> &permutation); 558 559 /** 560 * Given an arbitrary object of type T, use boost::serialization utilities 561 * to pack the object into a vector of characters and append it to the 562 * given buffer. The number of elements that have been added to the buffer 563 * will be returned. The object can be unpacked using the Utilities::unpack 564 * function below. 565 * 566 * If the library has been compiled with ZLIB enabled, then the output buffer 567 * can be compressed. This can be triggered with the parameter 568 * @p allow_compression, and is only of effect if ZLIB is enabled. 569 * 570 * If many consecutive calls with the same buffer are considered, it is 571 * recommended for reasons of performance to ensure that its capacity is 572 * sufficient. 573 */ 574 template <typename T> 575 size_t 576 pack(const T & object, 577 std::vector<char> &dest_buffer, 578 const bool allow_compression = true); 579 580 /** 581 * Creates and returns a buffer solely for the given object, using the 582 * above mentioned pack function. 583 * 584 * If the library has been compiled with ZLIB enabled, then the output buffer 585 * can be compressed. This can be triggered with the parameter 586 * @p allow_compression, and is only of effect if ZLIB is enabled. 587 */ 588 template <typename T> 589 std::vector<char> 590 pack(const T &object, const bool allow_compression = true); 591 592 /** 593 * Given a vector of characters, obtained through a call to the function 594 * Utilities::pack, restore its content in an object of type T. 595 * 596 * This function uses boost::serialization utilities to unpack the object 597 * from a vector of characters, and it is the inverse of the function 598 * Utilities::pack(). 599 * 600 * The @p allow_compression parameter denotes if the buffer to 601 * read from could have been previously compressed with ZLIB, and 602 * is only of effect if ZLIB is enabled. 603 * 604 * @note Since no arguments to this function depend on the template type 605 * @p T, you must manually specify the template argument when calling 606 * this function. 607 * 608 * @note If you want to pack() or unpack() arrays of objects, then the 609 * following works: 610 * @code 611 * double array[3] = {1,2,3}; 612 * std::vector<char> buffer = Utilities::pack(array); 613 * @endcode 614 * However, the converse does not: 615 * @code 616 * array = Utilities::unpack<double[3]>(buffer); 617 * @endcode 618 * This is because C++ does not allow functions to return arrays. 619 * Consequently, there is a separate unpack() function for arrays, see 620 * below. 621 */ 622 template <typename T> 623 T 624 unpack(const std::vector<char> &buffer, const bool allow_compression = true); 625 626 /** 627 * Same unpack function as above, but takes constant iterators on 628 * (a fraction of) a given packed buffer of type std::vector<char> instead. 629 * 630 * The @p allow_compression parameter denotes if the buffer to 631 * read from could have been previously compressed with ZLIB, and 632 * is only of effect if ZLIB is enabled. 633 */ 634 template <typename T> 635 T 636 unpack(const std::vector<char>::const_iterator &cbegin, 637 const std::vector<char>::const_iterator &cend, 638 const bool allow_compression = true); 639 640 /** 641 * Given a vector of characters, obtained through a call to the function 642 * Utilities::pack, restore its content in an array of type T. 643 * 644 * This function uses boost::serialization utilities to unpack the object 645 * from a vector of characters, and it is the inverse of the function 646 * Utilities::pack(). 647 * 648 * The @p allow_compression parameter denotes if the buffer to 649 * read from could have been previously compressed with ZLIB, and 650 * is only of effect if ZLIB is enabled. 651 * 652 * @note This function exists due to a quirk of C++. Specifically, 653 * if you want to pack() or unpack() arrays of objects, then the 654 * following works: 655 * @code 656 * double array[3] = {1,2,3}; 657 * std::vector<char> buffer = Utilities::pack(array); 658 * @endcode 659 * However, the converse does not: 660 * @code 661 * array = Utilities::unpack<double[3]>(buffer); 662 * @endcode 663 * This is because C++ does not allow functions to return arrays. 664 * The current function therefore allows to write 665 * @code 666 * Utilities::unpack(buffer, array); 667 * @endcode 668 * Note that unlike the other unpack() function, it is not necessary 669 * to explicitly specify the template arguments since they can be 670 * deduced from the second argument. 671 */ 672 template <typename T, int N> 673 void 674 unpack(const std::vector<char> &buffer, 675 T (&unpacked_object)[N], 676 const bool allow_compression = true); 677 678 /** 679 * Same unpack function as above, but takes constant iterators on 680 * (a fraction of) a given packed buffer of type std::vector<char> instead. 681 * 682 * The @p allow_compression parameter denotes if the buffer to 683 * read from could have been previously compressed with ZLIB, and 684 * is only of effect if ZLIB is enabled. 685 */ 686 template <typename T, int N> 687 void 688 unpack(const std::vector<char>::const_iterator &cbegin, 689 const std::vector<char>::const_iterator &cend, 690 T (&unpacked_object)[N], 691 const bool allow_compression = true); 692 693 /** 694 * Convert an object of type `std::unique_ptr<From>` to an object of 695 * type `std::unique_ptr<To>`, where it is assumed that we can cast 696 * the pointer to `From` to a pointer to `To` using a `dynamic_cast` 697 * -- in other words, we assume that `From` and `To` are connected 698 * through a class hierarchy, and that the object pointed to is in 699 * fact of a type that contains both a `From` and a `To`. An example 700 * is if either `To` is derived from `From` or the other way around. 701 * 702 * The function throws an exception of type `std::bad_cast` if the 703 * `dynamic_cast` does not succeed. This is the same exception you 704 * would get if a regular `dynamic_cast` between object types (but not 705 * pointer types) does not succeed. 706 * 707 * An example of how this function works is as follows: 708 * @code 709 * // A base class. Assume that it has virtual 710 * // functions so that dynamic_cast can work. 711 * class B 712 * { 713 * ... 714 * }; 715 * 716 * // A derived class 717 * class D : public B 718 * { 719 * ... 720 * }; 721 * 722 * // A factory function 723 * std::unique_ptr<B> create_object (...) 724 * { 725 * ... 726 * } 727 * 728 * void foo (...) 729 * { 730 * std::unique_ptr<B> b = create_object (...); 731 * 732 * // Assume that we know for some reason that the object above must 733 * // have created a D object but returned it as a std::unique_ptr<B>. 734 * // In order to access the D functionality, we need to cast the 735 * // pointer. Use the equivalent to dynamic_cast: 736 * std::unique_ptr<D> d = dynamic_unique_cast<D>(std::move(b)); 737 * 738 * // If the object really was a D, then 'd' now points to it. Note 739 * // also that in accordance with the semantics of std::unique_ptr, 740 * // it was necessary to std::move the 'b' object, and indeed 'b' 741 * // now no longer points to anything -- ownership has been 742 * // transferred to 'd'! 743 * @endcode 744 * 745 * @note This function does not try to convert the `Deleter` objects stored 746 * by `std::unique_ptr` objects. The function therefore only works if the 747 * deleter objects are at their defaults, i.e., if they are of type 748 * `std::default_delete<To>` and `std::default_delete<From>`. 749 */ 750 template <typename To, typename From> 751 std::unique_ptr<To> dynamic_unique_cast(std::unique_ptr<From> && p)752 dynamic_unique_cast(std::unique_ptr<From> &&p) 753 { 754 // Let's see if we can cast from 'From' to 'To'. If so, do the cast, 755 // and then release the pointer from the old 756 // owner 757 if (To *cast = dynamic_cast<To *>(p.get())) 758 { 759 std::unique_ptr<To> result(cast); 760 p.release(); 761 return result; 762 } 763 else 764 throw std::bad_cast(); 765 } 766 767 /** 768 * A namespace for utility functions that probe system properties. 769 * 770 * @ingroup utilities 771 */ 772 namespace System 773 { 774 /** 775 * Return the CPU load as returned by "uptime". Note that the 776 * interpretation of this number depends on the actual number of 777 * processors in the machine. This is presently only implemented on Linux, 778 * using the /proc/loadavg pseudo-file, on other systems we simply return 779 * zero. 780 */ 781 double 782 get_cpu_load(); 783 784 /** 785 * Return the instruction set extension for vectorization as described by 786 * DEAL_II_VECTORIZATION_WIDTH_IN_BITS in vectorization.h as a string. The 787 * list of possible return values is: 788 * 789 * <table> 790 * <tr> 791 * <td><tt>VECTORIZATION_LEVEL</tt></td> 792 * <td>Return Value</td> 793 * <td>Width in bits</td> 794 * </tr> 795 * <tr> 796 * <td>0</td> 797 * <td>disabled</td> 798 * <td>64</td> 799 * </tr> 800 * <tr> 801 * <td>1</td> 802 * <td>SSE2/AltiVec</td> 803 * <td>128</td> 804 * </tr> 805 * <tr> 806 * <td>2</td> 807 * <td>AVX</td> 808 * <td>256</td> 809 * </tr> 810 * <tr> 811 * <td>3</td> 812 * <td>AVX512</td> 813 * <td>512</td> 814 * </tr> 815 * </table> 816 */ 817 const std::string 818 get_current_vectorization_level(); 819 820 /** 821 * Structure that hold information about memory usage in kB. Used by 822 * get_memory_stats(). See man 5 proc entry /status for details. 823 */ 824 struct MemoryStats 825 { 826 /** 827 * Peak virtual memory size in kB. 828 */ 829 unsigned long int VmPeak; 830 831 /** 832 * Current virtual memory size in kB. 833 */ 834 unsigned long int VmSize; 835 836 /** 837 * Peak resident memory size in kB. Also known as "high water mark" (HWM). 838 */ 839 unsigned long int VmHWM; 840 841 /** 842 * Current resident memory size in kB. Also known as "resident set size" 843 * (RSS). 844 */ 845 unsigned long int VmRSS; 846 }; 847 848 849 /** 850 * Fill the @p stats structure with information about the memory 851 * consumption of this process. This is only implemented on Linux. 852 */ 853 void 854 get_memory_stats(MemoryStats &stats); 855 856 857 /** 858 * Return the name of the host this process runs on. 859 */ 860 std::string 861 get_hostname(); 862 863 864 /** 865 * Return the present time as HH:MM:SS. 866 */ 867 std::string 868 get_time(); 869 870 /** 871 * Return the present date as YYYY/MM/DD. MM and DD may be either one or 872 * two digits. 873 */ 874 std::string 875 get_date(); 876 877 /** 878 * Call the system function posix_memalign, or a replacement function if 879 * not available, to allocate memory with a certain minimal alignment. The 880 * first argument will then return a pointer to this memory block that can 881 * be released later on through a standard <code>free</code> call. 882 * 883 * @param memptr The address of a pointer variable that will after this 884 * call point to the allocated memory. 885 * @param alignment The minimal alignment of the memory block, in bytes. 886 * @param size The size of the memory block to be allocated, in bytes. 887 * 888 * @note This function checks internally for error codes, rather than 889 * leaving this task to the calling site. 890 */ 891 void 892 posix_memalign(void **memptr, std::size_t alignment, std::size_t size); 893 } // namespace System 894 895 896 #ifdef DEAL_II_WITH_TRILINOS 897 /** 898 * This namespace provides some of the basic structures used in the 899 * initialization of the Trilinos objects (e.g., matrices, vectors, and 900 * preconditioners). 901 */ 902 namespace Trilinos 903 { 904 /** 905 * Return a Trilinos Epetra_Comm object needed for creation of 906 * Epetra_Maps. 907 * 908 * If deal.II has been configured to use a compiler that does not support 909 * MPI then the resulting communicator will be a serial one. Otherwise, 910 * the communicator will correspond to MPI_COMM_WORLD, i.e. a communicator 911 * that encompasses all processes within this MPI universe. 912 */ 913 const Epetra_Comm & 914 comm_world(); 915 916 /** 917 * Return a Trilinos Epetra_Comm object needed for creation of 918 * Epetra_Maps. 919 * 920 * If deal.II has been configured to use a compiler that does not support 921 * MPI then the resulting communicator will be a serial one. Otherwise, 922 * the communicator will correspond to MPI_COMM_SELF, i.e. a communicator 923 * that comprises only this one processor. 924 */ 925 const Epetra_Comm & 926 comm_self(); 927 928 /** 929 * Return a Teuchos::Comm object needed for creation of Tpetra::Maps. 930 * 931 * If deal.II has been configured to use a compiler that does not support 932 * MPI then the resulting communicator will be a serial one. Otherwise, 933 * the communicator will correspond to MPI_COMM_SELF, i.e. a communicator 934 * that comprises only this one processor. 935 */ 936 const Teuchos::RCP<const Teuchos::Comm<int>> & 937 tpetra_comm_self(); 938 939 /** 940 * Given a communicator, duplicate it. If the given communicator is 941 * serial, that means to just return a copy of itself. On the other hand, 942 * if it is %parallel, we duplicate the underlying MPI_Comm object: we 943 * create a separate MPI communicator that contains the same processors 944 * and in the same order but has a separate identifier distinct from the 945 * given communicator. The function returns a pointer to a new object of a 946 * class derived from Epetra_Comm. The caller of this function needs to 947 * assume ownership of this function. The returned object should be 948 * destroyed using the destroy_communicator() function. 949 * 950 * This facility is used to separate streams of communication. For 951 * example, a program could simply use MPI_Comm_World for everything. But 952 * it is easy to come up with scenarios where sometimes not all processors 953 * participate in a communication that is intended to be global -- for 954 * example if we assemble a matrix on a coarse mesh with fewer cells than 955 * there are processors, some processors may not sync their matrices with 956 * the rest because they haven't written into it because they own no 957 * cells. That's clearly a bug. However, if these processors just continue 958 * their work, and the next %parallel operation happens to be a sync on a 959 * different matrix, then the sync could succeed -- by accident, since 960 * different processors are talking about different matrices. 961 * 962 * This kind of situation can be avoided if we use different communicators 963 * for different matrices which reduces the likelihood that communications 964 * meant to be separate aren't recognized as such just because they happen 965 * on the same communicator. In addition, it is conceivable that some MPI 966 * operations can be parallelized using multiple threads because their 967 * communicators identifies the communication in question, not their 968 * relative timing as is the case in a sequential program that just uses a 969 * single communicator. 970 */ 971 Epetra_Comm * 972 duplicate_communicator(const Epetra_Comm &communicator); 973 974 /** 975 * Given an Epetra communicator that was created by the 976 * duplicate_communicator() function, destroy the underlying MPI 977 * communicator object and reset the Epetra_Comm object to a the result of 978 * comm_self(). 979 * 980 * It is necessary to call this function at the time when the result of 981 * duplicate_communicator() is no longer needed. The reason is that in 982 * that function, we first create a new MPI_Comm object and then create an 983 * Epetra_Comm around it. While we can take care of destroying the latter, 984 * it doesn't destroy the communicator since it can only assume that it 985 * may also be still used by other objects in the program. Consequently, 986 * we have to take care of destroying it ourselves, explicitly. 987 * 988 * This function does exactly that. Because this has to happen while the 989 * Epetra_Comm object is still around, it first resets the latter and then 990 * destroys the communicator object. 991 * 992 * @note If you call this function on an Epetra_Comm object that is not 993 * created by duplicate_communicator(), you are likely doing something 994 * quite wrong. Don't do this. 995 */ 996 void 997 destroy_communicator(Epetra_Comm &communicator); 998 999 /** 1000 * Return the number of MPI processes there exist in the given 1001 * @ref GlossMPICommunicator "communicator" 1002 * object. If this is a sequential job (i.e., the program 1003 * is not using MPI at all, or is using MPI but has been started with 1004 * only one MPI process), then the communicator necessarily involves 1005 * only one process and the function returns 1. 1006 */ 1007 unsigned int 1008 get_n_mpi_processes(const Epetra_Comm &mpi_communicator); 1009 1010 /** 1011 * Return the number of the present MPI process in the space of processes 1012 * described by the given communicator. This will be a unique value for 1013 * each process between zero and (less than) the number of all processes 1014 * (given by get_n_mpi_processes()). 1015 */ 1016 unsigned int 1017 get_this_mpi_process(const Epetra_Comm &mpi_communicator); 1018 1019 /** 1020 * Given a Trilinos Epetra map, create a new map that has the same 1021 * subdivision of elements to processors but uses the given communicator 1022 * object instead of the one stored in the first argument. In essence, 1023 * this means that we create a map that communicates among the same 1024 * processors in the same way, but using a separate channel. 1025 * 1026 * This function is typically used with a communicator that has been 1027 * obtained by the duplicate_communicator() function. 1028 */ 1029 Epetra_Map 1030 duplicate_map(const Epetra_BlockMap &map, const Epetra_Comm &comm); 1031 } // namespace Trilinos 1032 1033 #endif 1034 1035 1036 } // namespace Utilities 1037 1038 1039 // --------------------- inline functions 1040 1041 namespace Utilities 1042 { 1043 template <int N, typename T> 1044 inline T fixed_power(const T x)1045 fixed_power(const T x) 1046 { 1047 Assert( 1048 !std::is_integral<T>::value || (N >= 0), 1049 ExcMessage( 1050 "The non-type template parameter N must be a non-negative integer for integral type T")); 1051 1052 if (N == 0) 1053 return T(1.); 1054 else if (N < 0) 1055 return T(1.) / fixed_power<-N>(x); 1056 else 1057 // Use exponentiation by squaring: 1058 return ((N % 2 == 1) ? x * fixed_power<N / 2>(x * x) : 1059 fixed_power<N / 2>(x * x)); 1060 } 1061 1062 1063 1064 template <class T> 1065 inline std::string type_to_string(const T & t)1066 type_to_string(const T &t) 1067 { 1068 return boost::core::demangle(typeid(t).name()); 1069 } 1070 1071 1072 1073 template <typename Iterator, typename T> 1074 inline Iterator lower_bound(Iterator first,Iterator last,const T & val)1075 lower_bound(Iterator first, Iterator last, const T &val) 1076 { 1077 return Utilities::lower_bound(first, last, val, std::less<T>()); 1078 } 1079 1080 1081 1082 template <typename Iterator, typename T, typename Comp> 1083 inline Iterator lower_bound(Iterator first,Iterator last,const T & val,const Comp comp)1084 lower_bound(Iterator first, Iterator last, const T &val, const Comp comp) 1085 { 1086 // verify that the two iterators are properly ordered. since 1087 // we need operator- for the iterator type anyway, do the 1088 // test as follows, rather than via 'last >= first' 1089 Assert(last - first >= 0, 1090 ExcMessage( 1091 "The given iterators do not satisfy the proper ordering.")); 1092 1093 unsigned int len = static_cast<unsigned int>(last - first); 1094 1095 if (len == 0) 1096 return first; 1097 1098 while (true) 1099 { 1100 // if length equals 8 or less, 1101 // then do a rolled out 1102 // search. use a switch without 1103 // breaks for that and roll-out 1104 // the loop somehow 1105 if (len < 8) 1106 { 1107 switch (len) 1108 { 1109 case 7: 1110 if (!comp(*first, val)) 1111 return first; 1112 ++first; 1113 DEAL_II_FALLTHROUGH; 1114 case 6: 1115 if (!comp(*first, val)) 1116 return first; 1117 ++first; 1118 DEAL_II_FALLTHROUGH; 1119 case 5: 1120 if (!comp(*first, val)) 1121 return first; 1122 ++first; 1123 DEAL_II_FALLTHROUGH; 1124 case 4: 1125 if (!comp(*first, val)) 1126 return first; 1127 ++first; 1128 DEAL_II_FALLTHROUGH; 1129 case 3: 1130 if (!comp(*first, val)) 1131 return first; 1132 ++first; 1133 DEAL_II_FALLTHROUGH; 1134 case 2: 1135 if (!comp(*first, val)) 1136 return first; 1137 ++first; 1138 DEAL_II_FALLTHROUGH; 1139 case 1: 1140 if (!comp(*first, val)) 1141 return first; 1142 return first + 1; 1143 default: 1144 // indices seem 1145 // to not be 1146 // sorted 1147 // correctly!? or 1148 // did len 1149 // become==0 1150 // somehow? that 1151 // shouldn't have 1152 // happened 1153 Assert(false, ExcInternalError()); 1154 } 1155 } 1156 1157 1158 1159 const unsigned int half = len >> 1; 1160 const Iterator middle = first + half; 1161 1162 // if the value is larger than 1163 // that pointed to by the 1164 // middle pointer, then the 1165 // insertion point must be 1166 // right of it 1167 if (comp(*middle, val)) 1168 { 1169 first = middle + 1; 1170 len -= half + 1; 1171 } 1172 else 1173 len = half; 1174 } 1175 } 1176 1177 1178 // --------------------- non-inline functions 1179 1180 template <typename T> 1181 size_t pack(const T & object,std::vector<char> & dest_buffer,const bool allow_compression)1182 pack(const T & object, 1183 std::vector<char> &dest_buffer, 1184 const bool allow_compression) 1185 { 1186 // the data is never compressed when we can't use zlib. 1187 (void)allow_compression; 1188 1189 std::size_t size = 0; 1190 1191 // see if the object is small and copyable via memcpy. if so, use 1192 // this fast path. otherwise, we have to go through the BOOST 1193 // serialization machinery 1194 // 1195 // we have to work around the fact that GCC 4.8.x claims to be C++ 1196 // conforming, but is not actually as it does not implement 1197 // std::is_trivially_copyable. 1198 #if __GNUG__ && __GNUC__ < 5 1199 if (__has_trivial_copy(T) && sizeof(T) < 256) 1200 #else 1201 # ifdef DEAL_II_HAVE_CXX17 1202 if constexpr (std::is_trivially_copyable<T>() && sizeof(T) < 256) 1203 # else 1204 if (std::is_trivially_copyable<T>() && sizeof(T) < 256) 1205 # endif 1206 #endif 1207 { 1208 const std::size_t previous_size = dest_buffer.size(); 1209 dest_buffer.resize(previous_size + sizeof(T)); 1210 1211 std::memcpy(dest_buffer.data() + previous_size, &object, sizeof(T)); 1212 1213 size = sizeof(T); 1214 } 1215 else 1216 { 1217 // use buffer as the target of a compressing 1218 // stream into which we serialize the current object 1219 const std::size_t previous_size = dest_buffer.size(); 1220 #ifdef DEAL_II_WITH_ZLIB 1221 if (allow_compression) 1222 { 1223 boost::iostreams::filtering_ostream out; 1224 out.push( 1225 boost::iostreams::gzip_compressor(boost::iostreams::gzip_params( 1226 boost::iostreams::gzip::default_compression))); 1227 out.push(boost::iostreams::back_inserter(dest_buffer)); 1228 1229 boost::archive::binary_oarchive archive(out); 1230 archive << object; 1231 out.flush(); 1232 } 1233 else 1234 #endif 1235 { 1236 std::ostringstream out; 1237 boost::archive::binary_oarchive archive(out); 1238 archive << object; 1239 1240 const std::string s = out.str(); 1241 dest_buffer.reserve(dest_buffer.size() + s.size()); 1242 std::move(s.begin(), s.end(), std::back_inserter(dest_buffer)); 1243 } 1244 1245 size = dest_buffer.size() - previous_size; 1246 } 1247 1248 return size; 1249 } 1250 1251 1252 template <typename T> 1253 std::vector<char> pack(const T & object,const bool allow_compression)1254 pack(const T &object, const bool allow_compression) 1255 { 1256 std::vector<char> buffer; 1257 pack<T>(object, buffer, allow_compression); 1258 return buffer; 1259 } 1260 1261 1262 template <typename T> 1263 T unpack(const std::vector<char>::const_iterator & cbegin,const std::vector<char>::const_iterator & cend,const bool allow_compression)1264 unpack(const std::vector<char>::const_iterator &cbegin, 1265 const std::vector<char>::const_iterator &cend, 1266 const bool allow_compression) 1267 { 1268 T object; 1269 1270 // the data is never compressed when we can't use zlib. 1271 (void)allow_compression; 1272 1273 // see if the object is small and copyable via memcpy. if so, use 1274 // this fast path. otherwise, we have to go through the BOOST 1275 // serialization machinery 1276 // 1277 // we have to work around the fact that GCC 4.8.x claims to be C++ 1278 // conforming, but is not actually as it does not implement 1279 // std::is_trivially_copyable. 1280 #if __GNUG__ && __GNUC__ < 5 1281 if (__has_trivial_copy(T) && sizeof(T) < 256) 1282 #else 1283 # ifdef DEAL_II_HAVE_CXX17 1284 if constexpr (std::is_trivially_copyable<T>() && sizeof(T) < 256) 1285 # else 1286 if (std::is_trivially_copyable<T>() && sizeof(T) < 256) 1287 # endif 1288 #endif 1289 { 1290 Assert(std::distance(cbegin, cend) == sizeof(T), ExcInternalError()); 1291 std::memcpy(&object, &*cbegin, sizeof(T)); 1292 } 1293 else 1294 { 1295 std::string decompressed_buffer; 1296 1297 // first decompress the buffer 1298 #ifdef DEAL_II_WITH_ZLIB 1299 if (allow_compression) 1300 { 1301 boost::iostreams::filtering_ostream decompressing_stream; 1302 decompressing_stream.push(boost::iostreams::gzip_decompressor()); 1303 decompressing_stream.push( 1304 boost::iostreams::back_inserter(decompressed_buffer)); 1305 decompressing_stream.write(&*cbegin, std::distance(cbegin, cend)); 1306 } 1307 else 1308 #endif 1309 { 1310 decompressed_buffer.assign(cbegin, cend); 1311 } 1312 1313 // then restore the object from the buffer 1314 std::istringstream in(decompressed_buffer); 1315 boost::archive::binary_iarchive archive(in); 1316 1317 archive >> object; 1318 } 1319 1320 return object; 1321 } 1322 1323 1324 template <typename T> 1325 T unpack(const std::vector<char> & buffer,const bool allow_compression)1326 unpack(const std::vector<char> &buffer, const bool allow_compression) 1327 { 1328 return unpack<T>(buffer.cbegin(), buffer.cend(), allow_compression); 1329 } 1330 1331 1332 template <typename T, int N> 1333 void unpack(const std::vector<char>::const_iterator & cbegin,const std::vector<char>::const_iterator & cend,T (& unpacked_object)[N],const bool allow_compression)1334 unpack(const std::vector<char>::const_iterator &cbegin, 1335 const std::vector<char>::const_iterator &cend, 1336 T (&unpacked_object)[N], 1337 const bool allow_compression) 1338 { 1339 // see if the object is small and copyable via memcpy. if so, use 1340 // this fast path. otherwise, we have to go through the BOOST 1341 // serialization machinery 1342 // 1343 // we have to work around the fact that GCC 4.8.x claims to be C++ 1344 // conforming, but is not actually as it does not implement 1345 // std::is_trivially_copyable. 1346 if ( 1347 #if __GNUG__ && __GNUC__ < 5 1348 __has_trivial_copy(T) 1349 #else 1350 std::is_trivially_copyable<T>() 1351 #endif 1352 && sizeof(T) * N < 256) 1353 { 1354 Assert(std::distance(cbegin, cend) == sizeof(T) * N, 1355 ExcInternalError()); 1356 std::memcpy(unpacked_object, &*cbegin, sizeof(T) * N); 1357 } 1358 else 1359 { 1360 std::string decompressed_buffer; 1361 1362 // first decompress the buffer 1363 (void)allow_compression; 1364 #ifdef DEAL_II_WITH_ZLIB 1365 if (allow_compression) 1366 { 1367 boost::iostreams::filtering_ostream decompressing_stream; 1368 decompressing_stream.push(boost::iostreams::gzip_decompressor()); 1369 decompressing_stream.push( 1370 boost::iostreams::back_inserter(decompressed_buffer)); 1371 decompressing_stream.write(&*cbegin, std::distance(cbegin, cend)); 1372 } 1373 else 1374 #endif 1375 { 1376 decompressed_buffer.assign(cbegin, cend); 1377 } 1378 1379 // then restore the object from the buffer 1380 std::istringstream in(decompressed_buffer); 1381 boost::archive::binary_iarchive archive(in); 1382 1383 archive >> unpacked_object; 1384 } 1385 } 1386 1387 1388 template <typename T, int N> 1389 void unpack(const std::vector<char> & buffer,T (& unpacked_object)[N],const bool allow_compression)1390 unpack(const std::vector<char> &buffer, 1391 T (&unpacked_object)[N], 1392 const bool allow_compression) 1393 { 1394 unpack<T, N>(buffer.cbegin(), 1395 buffer.cend(), 1396 unpacked_object, 1397 allow_compression); 1398 } 1399 1400 1401 1402 template <typename Integer> 1403 std::vector<Integer> reverse_permutation(const std::vector<Integer> & permutation)1404 reverse_permutation(const std::vector<Integer> &permutation) 1405 { 1406 const std::size_t n = permutation.size(); 1407 1408 std::vector<Integer> out(n); 1409 for (std::size_t i = 0; i < n; ++i) 1410 out[i] = n - 1 - permutation[i]; 1411 1412 return out; 1413 } 1414 1415 1416 1417 template <typename Integer> 1418 std::vector<Integer> invert_permutation(const std::vector<Integer> & permutation)1419 invert_permutation(const std::vector<Integer> &permutation) 1420 { 1421 const std::size_t n = permutation.size(); 1422 1423 std::vector<Integer> out(n, numbers::invalid_unsigned_int); 1424 1425 for (std::size_t i = 0; i < n; ++i) 1426 { 1427 AssertIndexRange(permutation[i], n); 1428 out[permutation[i]] = i; 1429 } 1430 1431 // check that we have actually reached 1432 // all indices 1433 for (std::size_t i = 0; i < n; ++i) 1434 Assert(out[i] != numbers::invalid_unsigned_int, 1435 ExcMessage("The given input permutation had duplicate entries!")); 1436 1437 return out; 1438 } 1439 } // namespace Utilities 1440 1441 1442 DEAL_II_NAMESPACE_CLOSE 1443 1444 #ifndef DOXYGEN 1445 namespace boost 1446 { 1447 namespace serialization 1448 { 1449 // Provides boost and c++11 with a way to serialize tuples and pairs 1450 // automatically 1451 template <int N> 1452 struct Serialize 1453 { 1454 template <class Archive, typename... Args> 1455 static void serializeSerialize1456 serialize(Archive &ar, std::tuple<Args...> &t, const unsigned int version) 1457 { 1458 ar &std::get<N - 1>(t); 1459 Serialize<N - 1>::serialize(ar, t, version); 1460 } 1461 }; 1462 1463 template <> 1464 struct Serialize<0> 1465 { 1466 template <class Archive, typename... Args> 1467 static void 1468 serialize(Archive &ar, std::tuple<Args...> &t, const unsigned int version) 1469 { 1470 (void)ar; 1471 (void)t; 1472 (void)version; 1473 } 1474 }; 1475 1476 template <class Archive, typename... Args> 1477 void 1478 serialize(Archive &ar, std::tuple<Args...> &t, const unsigned int version) 1479 { 1480 Serialize<sizeof...(Args)>::serialize(ar, t, version); 1481 } 1482 } // namespace serialization 1483 } // namespace boost 1484 #endif 1485 1486 #endif 1487