Reference documentation for deal.II version 9.4.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
utilities.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2005 - 2022 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#include <deal.II/base/config.h>
17
18// It's necessary to include winsock2.h before thread_local_storage.h,
19// because Intel implementation of TBB includes winsock.h,
20// and we'll get a conflict between winsock.h and winsock2.h otherwise.
21#ifdef DEAL_II_MSVC
22# include <winsock2.h>
23#endif
24
26#include <deal.II/base/mpi.h>
27#include <deal.II/base/point.h>
30
32#define BOOST_BIND_GLOBAL_PLACEHOLDERS
33#include <boost/archive/iterators/base64_from_binary.hpp>
34#include <boost/archive/iterators/binary_from_base64.hpp>
35#include <boost/archive/iterators/transform_width.hpp>
36#include <boost/iostreams/copy.hpp>
37#include <boost/lexical_cast.hpp>
38#include <boost/random.hpp>
39#undef BOOST_BIND_GLOBAL_PLACEHOLDERS
41
42#include <algorithm>
43#include <bitset>
44#include <cctype>
45#include <cerrno>
46#include <cmath>
47#include <cstddef>
48#include <cstdio>
49#include <ctime>
50#include <fstream>
51#include <iomanip>
52#include <iostream>
53#include <limits>
54#include <sstream>
55#include <string>
56
57#if defined(DEAL_II_HAVE_UNISTD_H) && defined(DEAL_II_HAVE_GETHOSTNAME)
58# include <unistd.h>
59#endif
60
61#ifndef DEAL_II_MSVC
62# include <cstdlib>
63#endif
64
65
66#ifdef DEAL_II_WITH_TRILINOS
67# ifdef DEAL_II_WITH_MPI
71
72# include <Epetra_MpiComm.h>
73# include <Teuchos_DefaultComm.hpp>
74# endif
75# include <Epetra_SerialComm.h>
76# include <Teuchos_RCP.hpp>
77#endif
78
80
81
82namespace Utilities
83{
85 unsigned int,
86 unsigned int,
87 << "When trying to convert " << arg1 << " to a string with "
88 << arg2 << " digits");
89 DeclException1(ExcInvalidNumber, unsigned int, << "Invalid number " << arg1);
91 std::string,
92 << "Can't convert the string " << arg1
93 << " to the desired type");
94
95
96 std::string
98 {
100 }
101
102
103 namespace
104 {
105 template <int dim,
106 typename Number,
107 int effective_dim,
108 typename LongDouble,
109 typename Integer>
110 std::vector<std::array<std::uint64_t, effective_dim>>
111 inverse_Hilbert_space_filling_curve_effective(
112 const std::vector<Point<dim, Number>> &points,
113 const Point<dim, Number> & bl,
114 const std::array<LongDouble, dim> & extents,
115 const std::bitset<dim> & valid_extents,
116 const int min_bits,
117 const Integer max_int)
118 {
119 std::vector<std::array<Integer, effective_dim>> int_points(points.size());
120
121 for (unsigned int i = 0; i < points.size(); ++i)
122 {
123 // convert into integers:
124 unsigned int eff_d = 0;
125 for (unsigned int d = 0; d < dim; ++d)
126 if (valid_extents[d])
127 {
128 Assert(extents[d] > 0, ExcInternalError());
129 const LongDouble v = (static_cast<LongDouble>(points[i][d]) -
130 static_cast<LongDouble>(bl[d])) /
131 extents[d];
132 Assert(v >= 0. && v <= 1., ExcInternalError());
133 AssertIndexRange(eff_d, effective_dim);
134 int_points[i][eff_d] =
135 static_cast<Integer>(v * static_cast<LongDouble>(max_int));
136 ++eff_d;
137 }
138 }
139
140 // note that we call this with "min_bits"
141 return inverse_Hilbert_space_filling_curve<effective_dim>(int_points,
142 min_bits);
143 }
144 } // namespace
145
146 template <int dim, typename Number>
147 std::vector<std::array<std::uint64_t, dim>>
149 const std::vector<Point<dim, Number>> &points,
150 const int bits_per_dim)
151 {
152 using Integer = std::uint64_t;
153 // take floating point number hopefully with mantissa >= 64bit
154 using LongDouble = long double;
155
156 // return if there is nothing to do
157 if (points.size() == 0)
158 return std::vector<std::array<std::uint64_t, dim>>();
159
160 // get bounding box:
161 Point<dim, Number> bl = points[0], tr = points[0];
162 for (const auto &p : points)
163 for (unsigned int d = 0; d < dim; ++d)
164 {
165 const double cid = p[d];
166 bl[d] = std::min(cid, bl[d]);
167 tr[d] = std::max(cid, tr[d]);
168 }
169
170 std::array<LongDouble, dim> extents;
171 std::bitset<dim> valid_extents;
172 for (unsigned int i = 0; i < dim; ++i)
173 {
174 extents[i] =
175 static_cast<LongDouble>(tr[i]) - static_cast<LongDouble>(bl[i]);
176 valid_extents[i] = (extents[i] > 0.);
177 }
178
179 // make sure our conversion from fractional coordinates to
180 // Integers work as expected, namely our cast (LongDouble)max_int
181 const int min_bits =
182 std::min(bits_per_dim,
183 std::min(std::numeric_limits<Integer>::digits,
184 std::numeric_limits<LongDouble>::digits));
185
186 // based on that get the maximum integer:
187 const Integer max_int = (min_bits == std::numeric_limits<Integer>::digits ?
189 (Integer(1) << min_bits) - 1);
190
191 const unsigned int effective_dim = valid_extents.count();
192 if (effective_dim == dim)
193 {
194 return inverse_Hilbert_space_filling_curve_effective<dim,
195 Number,
196 dim,
197 LongDouble,
198 Integer>(
199 points, bl, extents, valid_extents, min_bits, max_int);
200 }
201
202 // various degenerate cases
203 std::array<std::uint64_t, dim> zero_ind;
204 for (unsigned int d = 0; d < dim; ++d)
205 zero_ind[d] = 0;
206
207 std::vector<std::array<std::uint64_t, dim>> ind(points.size(), zero_ind);
208 // manually check effective_dim == 1 and effective_dim == 2
209 if (dim == 3 && effective_dim == 2)
210 {
211 const auto ind2 =
212 inverse_Hilbert_space_filling_curve_effective<dim,
213 Number,
214 2,
215 LongDouble,
216 Integer>(
217 points, bl, extents, valid_extents, min_bits, max_int);
218
219 for (unsigned int i = 0; i < ind.size(); ++i)
220 for (unsigned int d = 0; d < 2; ++d)
221 ind[i][d + 1] = ind2[i][d];
222
223 return ind;
224 }
225 else if (effective_dim == 1)
226 {
227 const auto ind1 =
228 inverse_Hilbert_space_filling_curve_effective<dim,
229 Number,
230 1,
231 LongDouble,
232 Integer>(
233 points, bl, extents, valid_extents, min_bits, max_int);
234
235 for (unsigned int i = 0; i < ind.size(); ++i)
236 ind[i][dim - 1] = ind1[i][0];
237
238 return ind;
239 }
240
241 // we should get here only if effective_dim == 0
242 Assert(effective_dim == 0, ExcInternalError());
243
244 // if the bounding box is degenerate in all dimensions,
245 // can't do much but exit gracefully by setting index according
246 // to the index of each point so that there is no re-ordering
247 for (unsigned int i = 0; i < points.size(); ++i)
248 ind[i][dim - 1] = i;
249
250 return ind;
251 }
252
253
254
255 template <int dim>
256 std::vector<std::array<std::uint64_t, dim>>
258 const std::vector<std::array<std::uint64_t, dim>> &points,
259 const int bits_per_dim)
260 {
261 using Integer = std::uint64_t;
262
263 std::vector<std::array<Integer, dim>> int_points(points);
264
265 std::vector<std::array<Integer, dim>> res(int_points.size());
266
267 // follow
268 // J. Skilling, Programming the Hilbert curve, AIP Conf. Proc. 707, 381
269 // (2004); http://dx.doi.org/10.1063/1.1751381 also see
270 // https://stackoverflow.com/questions/499166/mapping-n-dimensional-value-to-a-point-on-hilbert-curve
271 // https://gitlab.com/octopus-code/octopus/blob/develop/src/grid/hilbert.c
272 // https://github.com/trilinos/Trilinos/blob/master/packages/zoltan/src/hsfc/hsfc_hilbert.c
273 // (Zoltan_HSFC_InvHilbertXd)
274 // https://github.com/aditi137/Hilbert/blob/master/Hilbert/hilbert.cpp
275
276 // now we can map to 1D coordinate stored in Transpose format
277 // adopt AxestoTranspose function from the paper, that
278 // transforms in-place between geometrical axes and Hilbert transpose.
279 // Example: b=5 bits for each of n=3 coordinates.
280 // 15-bit Hilbert integer = A B C D E F G H I J K L M N O is
281 // stored as its Transpose
282 // X[0] = A D G J M X[2]|
283 // X[1] = B E H K N <-------> | /X[1]
284 // X[2] = C F I L O axes |/
285 // high low 0------ X[0]
286
287 // Depth of the Hilbert curve
288 Assert(bits_per_dim <= std::numeric_limits<Integer>::digits,
289 ExcMessage("This integer type can not hold " +
290 std::to_string(bits_per_dim) + " bits."));
291
292 const Integer M = Integer(1) << (bits_per_dim - 1); // largest bit
293
294 for (unsigned int index = 0; index < int_points.size(); ++index)
295 {
296 auto &X = int_points[index];
297 auto &L = res[index];
298
299 // Inverse undo
300 for (Integer q = M; q > 1; q >>= 1)
301 {
302 const Integer p = q - 1;
303 for (unsigned int i = 0; i < dim; ++i)
304 {
305 // invert
306 if (X[i] & q)
307 {
308 X[0] ^= p;
309 }
310 // exchange
311 else
312 {
313 const Integer t = (X[0] ^ X[i]) & p;
314 X[0] ^= t;
315 X[i] ^= t;
316 }
317 }
318 }
319
320 // Gray encode (inverse of decode)
321 for (unsigned int i = 1; i < dim; ++i)
322 X[i] ^= X[i - 1];
323
324 Integer t = 0;
325 for (Integer q = M; q > 1; q >>= 1)
326 if (X[dim - 1] & q)
327 t ^= q - 1;
328 for (unsigned int i = 0; i < dim; ++i)
329 X[i] ^= t;
330
331 // now we need to go from index stored in transpose format to
332 // consecutive format, which is better suited for comparators.
333 // we could interleave into some big unsigned int...
334 // https://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/
335 // https://stackoverflow.com/questions/4431522/given-2-16-bit-ints-can-i-interleave-those-bits-to-form-a-single-32-bit-int
336 // ...but we would loose spatial resolution!
337
338 // interleave using brute force, follow TransposetoLine from
339 // https://github.com/aditi137/Hilbert/blob/master/Hilbert/hilbert.cpp
340 {
341 Integer p = M;
342 unsigned int j = 0;
343 for (unsigned int i = 0; i < dim; ++i)
344 {
345 L[i] = 0;
346 // go through bits using a mask q
347 for (Integer q = M; q > 0; q >>= 1)
348 {
349 if (X[j] & p)
350 L[i] |= q;
351 if (++j == dim)
352 {
353 j = 0;
354 p >>= 1;
355 }
356 }
357 }
358 }
359
360 } // end of the loop over points
361
362 return res;
363 }
364
365
366
367 template <int dim>
368 std::uint64_t
369 pack_integers(const std::array<std::uint64_t, dim> &index,
370 const int bits_per_dim)
371 {
372 using Integer = std::uint64_t;
373
374 AssertIndexRange(bits_per_dim * dim, 65);
375 Assert(bits_per_dim > 0, ExcMessage("bits_per_dim should be positive"));
376
377 const Integer mask = (Integer(1) << bits_per_dim) - 1;
378
379 Integer res = 0;
380 for (unsigned int i = 0; i < dim; ++i)
381 {
382 // take bits_per_dim from each integer and shift them
383 const Integer v = (mask & index[dim - 1 - i]) << (bits_per_dim * i);
384 res |= v;
385 }
386 return res;
387 }
388
389
390
391 std::string
392 compress(const std::string &input)
393 {
394#ifdef DEAL_II_WITH_ZLIB
395 namespace bio = boost::iostreams;
396
397 std::stringstream compressed;
398 std::stringstream origin(input);
399
400 bio::filtering_streambuf<bio::input> out;
401 out.push(bio::gzip_compressor());
402 out.push(origin);
403 bio::copy(out, compressed);
404
405 return compressed.str();
406#else
407 return input;
408#endif
409 }
410
411
412
413 std::string
414 decompress(const std::string &compressed_input)
415 {
416#ifdef DEAL_II_WITH_ZLIB
417 namespace bio = boost::iostreams;
418
419 std::stringstream compressed(compressed_input);
420 std::stringstream decompressed;
421
422 bio::filtering_streambuf<bio::input> out;
423 out.push(bio::gzip_decompressor());
424 out.push(compressed);
425 bio::copy(out, decompressed);
426
427 return decompressed.str();
428#else
429 return compressed_input;
430#endif
431 }
432
433
434
435 std::string
436 encode_base64(const std::vector<unsigned char> &binary_input)
437 {
438 using namespace boost::archive::iterators;
439 using It = base64_from_binary<
440 transform_width<std::vector<unsigned char>::const_iterator, 6, 8>>;
441 auto base64 = std::string(It(binary_input.begin()), It(binary_input.end()));
442 // Add padding.
443 return base64.append((3 - binary_input.size() % 3) % 3, '=');
444 }
445
446
447
448 std::vector<unsigned char>
449 decode_base64(const std::string &base64_input)
450 {
451 using namespace boost::archive::iterators;
452 using It =
453 transform_width<binary_from_base64<std::string::const_iterator>, 8, 6>;
454 auto binary = std::vector<unsigned char>(It(base64_input.begin()),
455 It(base64_input.end()));
456 // Remove padding.
457 auto length = base64_input.size();
458 if (binary.size() > 2 && base64_input[length - 1] == '=' &&
459 base64_input[length - 2] == '=')
460 {
461 binary.erase(binary.end() - 2, binary.end());
462 }
463 else if (binary.size() > 1 && base64_input[length - 1] == '=')
464 {
465 binary.erase(binary.end() - 1, binary.end());
466 }
467 return binary;
468 }
469
470
471
472 std::string
473 int_to_string(const unsigned int value, const unsigned int digits)
474 {
475 return to_string(value, digits);
476 }
477
478
479
480 template <typename number>
481 std::string
482 to_string(const number value, const unsigned int digits)
483 {
484 // For integer data types, use the standard std::to_string()
485 // function. On the other hand, that function is defined in terms
486 // of std::sprintf, which does not use the usual std::iostream
487 // interface and tries to render floating point numbers in awkward
488 // ways (see
489 // https://en.cppreference.com/w/cpp/string/basic_string/to_string). So
490 // resort to boost::lexical_cast for all other types (in
491 // particular for floating point types.
492 std::string lc_string = (std::is_integral<number>::value ?
493 std::to_string(value) :
494 boost::lexical_cast<std::string>(value));
495
496 if ((digits != numbers::invalid_unsigned_int) &&
497 (lc_string.size() < digits))
498 {
499 // We have to add the padding zeroes in front of the number
500 const unsigned int padding_position = (lc_string[0] == '-') ? 1 : 0;
501
502 const std::string padding(digits - lc_string.size(), '0');
503 lc_string.insert(padding_position, padding);
504 }
505
506 return lc_string;
507 }
508
509
510
511 std::string
512 replace_in_string(const std::string &input,
513 const std::string &from,
514 const std::string &to)
515 {
516 if (from.empty())
517 return input;
518
519 std::string out = input;
520 std::string::size_type pos = out.find(from);
521
522 while (pos != std::string::npos)
523 {
524 out.replace(pos, from.size(), to);
525 pos = out.find(from, pos + to.size());
526 }
527 return out;
528 }
529
530 std::string
531 trim(const std::string &input)
532 {
533 std::string::size_type left = 0;
534 std::string::size_type right = input.size() > 0 ? input.size() - 1 : 0;
535
536 for (; left < input.size(); ++left)
537 {
538 if (std::isspace(input[left]) == 0)
539 {
540 break;
541 }
542 }
543
544 for (; right >= left; --right)
545 {
546 if (std::isspace(input[right]) == 0)
547 {
548 break;
549 }
550 }
551
552 return std::string(input, left, right - left + 1);
553 }
554
555
556
557 std::string
558 dim_string(const int dim, const int spacedim)
559 {
560 if (dim == spacedim)
561 return int_to_string(dim);
562 else
563 return int_to_string(dim) + "," + int_to_string(spacedim);
564 }
565
566
567 unsigned int
568 needed_digits(const unsigned int max_number)
569 {
570 if (max_number > 0)
571 return static_cast<int>(
572 std::ceil(std::log10(std::fabs(max_number + 0.1))));
573
574 return 1;
575 }
576
577
578
579 template <typename Number>
580 Number
581 truncate_to_n_digits(const Number number, const unsigned int n_digits)
582 {
583 AssertThrow(n_digits >= 1, ExcMessage("invalid parameter."));
584
586 return number;
587
588 const int order =
589 static_cast<int>(std::floor(std::log10(std::fabs(number))));
590
591 const int shift = -order + static_cast<int>(n_digits) - 1;
592
593 Assert(shift <= static_cast<int>(std::floor(
596 "Overflow. Use a smaller value for n_digits and/or make sure "
597 "that the absolute value of 'number' does not become too small."));
598
599 const Number factor = std::pow(10.0, static_cast<Number>(shift));
600
601 const Number number_cutoff = std::trunc(number * factor) / factor;
602
603 return number_cutoff;
604 }
605
606
607 int
608 string_to_int(const std::string &s_)
609 {
610 // trim whitespace on either side of the text if necessary
611 std::string s = s_;
612 while ((s.size() > 0) && (s[0] == ' '))
613 s.erase(s.begin());
614 while ((s.size() > 0) && (s.back() == ' '))
615 s.erase(s.end() - 1);
616
617 // Now convert and see whether we succeed. Note that strtol only
618 // touches errno if an error occurred, so if we want to check
619 // whether an error happened, we need to make sure that errno==0
620 // before calling strtol since otherwise it may be that the
621 // conversion succeeds and that errno remains at the value it
622 // was before, whatever that was.
623 char *p;
624 errno = 0;
625 const int i = std::strtol(s.c_str(), &p, 10);
626
627 // We have an error if one of the following conditions is true:
628 // - strtol sets errno != 0
629 // - The original string was empty (we could have checked that
630 // earlier already)
631 // - The string has non-zero length and strtol converted the
632 // first part to something useful, but stopped converting short
633 // of the terminating '\0' character. This happens, for example,
634 // if the given string is "1234 abc".
635 AssertThrow(!((errno != 0) || (s.size() == 0) ||
636 ((s.size() > 0) && (*p != '\0'))),
637 ExcMessage("Can't convert <" + s + "> to an integer."));
638
639 return i;
640 }
641
642
643
644 std::vector<int>
645 string_to_int(const std::vector<std::string> &s)
646 {
647 std::vector<int> tmp(s.size());
648 for (unsigned int i = 0; i < s.size(); ++i)
649 tmp[i] = string_to_int(s[i]);
650 return tmp;
651 }
652
653
654
655 double
656 string_to_double(const std::string &s_)
657 {
658 // trim whitespace on either side of the text if necessary
659 std::string s = s_;
660 while ((s.size() > 0) && (s[0] == ' '))
661 s.erase(s.begin());
662 while ((s.size() > 0) && (s.back() == ' '))
663 s.erase(s.end() - 1);
664
665 // Now convert and see whether we succeed. Note that strtol only
666 // touches errno if an error occurred, so if we want to check
667 // whether an error happened, we need to make sure that errno==0
668 // before calling strtol since otherwise it may be that the
669 // conversion succeeds and that errno remains at the value it
670 // was before, whatever that was.
671 char *p;
672 errno = 0;
673 const double d = std::strtod(s.c_str(), &p);
674
675 // We have an error if one of the following conditions is true:
676 // - strtod sets errno != 0
677 // - The original string was empty (we could have checked that
678 // earlier already)
679 // - The string has non-zero length and strtod converted the
680 // first part to something useful, but stopped converting short
681 // of the terminating '\0' character. This happens, for example,
682 // if the given string is "1.234 abc".
683 AssertThrow(!((errno != 0) || (s.size() == 0) ||
684 ((s.size() > 0) && (*p != '\0'))),
685 ExcMessage("Can't convert <" + s + "> to a double."));
686
687 return d;
688 }
689
690
691
692 std::vector<double>
693 string_to_double(const std::vector<std::string> &s)
694 {
695 std::vector<double> tmp(s.size());
696 for (unsigned int i = 0; i < s.size(); ++i)
697 tmp[i] = string_to_double(s[i]);
698 return tmp;
699 }
700
701
702
703 std::vector<std::string>
704 split_string_list(const std::string &s, const std::string &delimiter)
705 {
706 // keep the currently remaining part of the input string in 'tmp' and
707 // keep chopping elements of the list off the front
708 std::string tmp = s;
709
710 // as discussed in the documentation, eat whitespace from the end
711 // of the string
712 while (tmp.size() != 0 && tmp.back() == ' ')
713 tmp.erase(tmp.size() - 1, 1);
714
715 // split the input list until it is empty. since in every iteration
716 // 'tmp' is what's left of the string after the next delimiter,
717 // and since we've stripped trailing space already, 'tmp' will
718 // be empty at one point if 's' ended in a delimiter, even if
719 // there was space after the last delimiter. this matches what's
720 // discussed in the documentation
721 std::vector<std::string> split_list;
722 while (tmp.size() != 0)
723 {
724 std::string name;
725 name = tmp;
726
727 if (name.find(delimiter) != std::string::npos)
728 {
729 name.erase(name.find(delimiter), std::string::npos);
730 tmp.erase(0, tmp.find(delimiter) + delimiter.size());
731 }
732 else
733 tmp = "";
734
735 // strip spaces from this element's front and end
736 while ((name.size() != 0) && (name[0] == ' '))
737 name.erase(0, 1);
738 while (name.size() != 0 && name.back() == ' ')
739 name.erase(name.size() - 1, 1);
740
741 split_list.push_back(name);
742 }
743
744 return split_list;
745 }
746
747
748 std::vector<std::string>
749 split_string_list(const std::string &s, const char delimiter)
750 {
751 std::string d = ",";
752 d[0] = delimiter;
753 return split_string_list(s, d);
754 }
755
756
757 std::vector<std::string>
758 break_text_into_lines(const std::string &original_text,
759 const unsigned int width,
760 const char delimiter)
761 {
762 std::string text = original_text;
763 std::vector<std::string> lines;
764
765 // remove trailing spaces
766 while ((text.size() != 0) && (text.back() == delimiter))
767 text.erase(text.size() - 1, 1);
768
769 // then split the text into lines
770 while (text.size() != 0)
771 {
772 // in each iteration, first remove
773 // leading spaces
774 while ((text.size() != 0) && (text[0] == delimiter))
775 text.erase(0, 1);
776
777 std::size_t pos_newline = text.find_first_of('\n', 0);
778 if (pos_newline != std::string::npos && pos_newline <= width)
779 {
780 std::string line(text, 0, pos_newline);
781 while ((line.size() != 0) && (line.back() == delimiter))
782 line.erase(line.size() - 1, 1);
783 lines.push_back(line);
784 text.erase(0, pos_newline + 1);
785 continue;
786 }
787
788 // if we can fit everything into one
789 // line, then do so. otherwise, we have
790 // to keep breaking
791 if (text.size() < width)
792 {
793 // remove trailing spaces
794 while ((text.size() != 0) && (text.back() == delimiter))
795 text.erase(text.size() - 1, 1);
796 lines.push_back(text);
797 text = "";
798 }
799 else
800 {
801 // starting at position width, find the
802 // location of the previous space, so
803 // that we can break around there
804 int location = std::min<int>(width, text.size() - 1);
805 for (; location > 0; --location)
806 if (text[location] == delimiter)
807 break;
808
809 // if there are no spaces, then try if
810 // there are spaces coming up
811 if (location == 0)
812 for (location = std::min<int>(width, text.size() - 1);
813 location < static_cast<int>(text.size());
814 ++location)
815 if (text[location] == delimiter)
816 break;
817
818 // now take the text up to the found
819 // location and put it into a single
820 // line, and remove it from 'text'
821 std::string line(text, 0, location);
822 while ((line.size() != 0) && (line.back() == delimiter))
823 line.erase(line.size() - 1, 1);
824 lines.push_back(line);
825 text.erase(0, location);
826 }
827 }
828
829 return lines;
830 }
831
832
833
834 bool
835 match_at_string_start(const std::string &name, const std::string &pattern)
836 {
837 if (pattern.size() > name.size())
838 return false;
839
840 for (unsigned int i = 0; i < pattern.size(); ++i)
841 if (pattern[i] != name[i])
842 return false;
843
844 return true;
845 }
846
847
848
849 std::pair<int, unsigned int>
850 get_integer_at_position(const std::string &name, const unsigned int position)
851 {
852 Assert(position < name.size(), ExcInternalError());
853
854 const std::string test_string(name.begin() + position, name.end());
855
856 std::istringstream str(test_string);
857
858 int i;
859 if (str >> i)
860 {
861 // compute the number of
862 // digits of i. assuming it
863 // is less than 8 is likely
864 // ok
865 if (i < 10)
866 return std::make_pair(i, 1U);
867 else if (i < 100)
868 return std::make_pair(i, 2U);
869 else if (i < 1000)
870 return std::make_pair(i, 3U);
871 else if (i < 10000)
872 return std::make_pair(i, 4U);
873 else if (i < 100000)
874 return std::make_pair(i, 5U);
875 else if (i < 1000000)
876 return std::make_pair(i, 6U);
877 else if (i < 10000000)
878 return std::make_pair(i, 7U);
879 else
880 {
881 Assert(false, ExcNotImplemented());
882 return std::make_pair(-1, numbers::invalid_unsigned_int);
883 }
884 }
885 else
886 return std::make_pair(-1, numbers::invalid_unsigned_int);
887 }
888
889
890
891 double
892 generate_normal_random_number(const double a, const double sigma)
893 {
894 // if no noise: return now
895 if (sigma == 0)
896 return a;
897
898 // we would want to use rand(), but that function is not reentrant
899 // in a thread context. one could use rand_r, but this does not
900 // produce reproducible results between threads either (though at
901 // least it is reentrant). these two approaches being
902 // non-workable, use a thread-local random number generator here.
903 // we could use std::mt19937 but doing so results in compiler-dependent
904 // output.
905 static Threads::ThreadLocalStorage<boost::mt19937> random_number_generator;
906 return boost::normal_distribution<>(a,
907 sigma)(random_number_generator.get());
908 }
909
910
911
912 namespace System
913 {
914#ifdef __linux__
915
916 double
918 {
919 std::ifstream cpuinfo;
920 cpuinfo.open("/proc/loadavg");
921
922 AssertThrow(cpuinfo.fail() == false, ExcIO());
923
924 double load;
925 cpuinfo >> load;
926
927 return load;
928 }
929
930#else
931
932 double
934 {
935 return 0.;
936 }
937
938#endif
939
940 const std::string
942 {
944 {
945 case 0:
946 return "disabled";
947 case 128:
948#ifdef __ALTIVEC__
949 return "AltiVec";
950#else
951 return "SSE2";
952#endif
953 case 256:
954 return "AVX";
955 case 512:
956 return "AVX512";
957 default:
958 AssertThrow(false,
960 "Invalid DEAL_II_VECTORIZATION_WIDTH_IN_BITS."));
961 return "ERROR";
962 }
963 }
964
965
966 void
968 {
969 stats.VmPeak = stats.VmSize = stats.VmHWM = stats.VmRSS = 0;
970
971 // parsing /proc/self/stat would be a
972 // lot easier, but it does not contain
973 // VmHWM, so we use /status instead.
974#ifdef __linux__
975 std::ifstream file("/proc/self/status");
976 std::string line;
977 std::string name;
978 while (!file.eof())
979 {
980 file >> name;
981 if (name == "VmPeak:")
982 file >> stats.VmPeak;
983 else if (name == "VmSize:")
984 file >> stats.VmSize;
985 else if (name == "VmHWM:")
986 file >> stats.VmHWM;
987 else if (name == "VmRSS:")
988 {
989 file >> stats.VmRSS;
990 break; // this is always the last entry
991 }
992
993 getline(file, line);
994 }
995#endif
996 }
997
998
999
1000 std::string
1002 {
1003#if defined(DEAL_II_HAVE_UNISTD_H) && defined(DEAL_II_HAVE_GETHOSTNAME)
1004 const unsigned int N = 1024;
1005 char hostname[N];
1006 gethostname(&(hostname[0]), N - 1);
1007#else
1008 std::string hostname("unknown");
1009#endif
1010 return hostname;
1011 }
1012
1013
1014
1015 std::string
1017 {
1018 std::time_t time1 = std::time(nullptr);
1019 std::tm * time = std::localtime(&time1);
1020
1021 std::ostringstream o;
1022 o << time->tm_hour << ":" << (time->tm_min < 10 ? "0" : "")
1023 << time->tm_min << ":" << (time->tm_sec < 10 ? "0" : "")
1024 << time->tm_sec;
1025
1026 return o.str();
1027 }
1028
1029
1030
1031 std::string
1033 {
1034 std::time_t time1 = std::time(nullptr);
1035 std::tm * time = std::localtime(&time1);
1036
1037 std::ostringstream o;
1038 o << time->tm_year + 1900 << "/" << time->tm_mon + 1 << "/"
1039 << time->tm_mday;
1040
1041 return o.str();
1042 }
1043
1044
1045
1046 void
1047 posix_memalign(void **memptr, std::size_t alignment, std::size_t size)
1048 {
1049#ifndef DEAL_II_MSVC
1050 const int ierr = ::posix_memalign(memptr, alignment, size);
1051
1052 AssertThrow(ierr == 0, ExcOutOfMemory(size));
1053 AssertThrow(*memptr != nullptr, ExcOutOfMemory(size));
1054#else
1055 // Windows does not appear to have posix_memalign. just use the
1056 // regular malloc in that case
1057 *memptr = malloc(size);
1058 (void)alignment;
1059 AssertThrow(*memptr != 0, ExcOutOfMemory(size));
1060#endif
1061 }
1062
1063
1064
1065 bool
1067 {
1069 }
1070 } // namespace System
1071
1072
1073#ifdef DEAL_II_WITH_TRILINOS
1074
1075 namespace Trilinos
1076 {
1077 const Epetra_Comm &
1079 {
1080# ifdef DEAL_II_WITH_MPI
1081 static Teuchos::RCP<Epetra_MpiComm> communicator =
1082 Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD), true);
1083# else
1084 static Teuchos::RCP<Epetra_SerialComm> communicator =
1085 Teuchos::rcp(new Epetra_SerialComm(), true);
1086# endif
1087
1088 return *communicator;
1089 }
1090
1091
1092
1093 const Teuchos::RCP<const Teuchos::Comm<int>> &
1095 {
1096# ifdef DEAL_II_WITH_MPI
1097 static auto communicator = Teuchos::RCP<const Teuchos::Comm<int>>(
1098 new Teuchos::MpiComm<int>(MPI_COMM_SELF));
1099# else
1100 static auto communicator =
1101 Teuchos::RCP<const Teuchos::Comm<int>>(new Teuchos::Comm<int>());
1102# endif
1103
1104 return communicator;
1105 }
1106
1107
1108
1109 const Epetra_Comm &
1111 {
1112# ifdef DEAL_II_WITH_MPI
1113 static Teuchos::RCP<Epetra_MpiComm> communicator =
1114 Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_SELF), true);
1115# else
1116 static Teuchos::RCP<Epetra_SerialComm> communicator =
1117 Teuchos::rcp(new Epetra_SerialComm(), true);
1118# endif
1119
1120 return *communicator;
1121 }
1122
1123
1124
1125 Epetra_Comm *
1126 duplicate_communicator(const Epetra_Comm &communicator)
1127 {
1128# ifdef DEAL_II_WITH_MPI
1129
1130 // see if the communicator is in fact a
1131 // parallel MPI communicator; if so,
1132 // return a duplicate of it
1133 const Epetra_MpiComm *mpi_comm =
1134 dynamic_cast<const Epetra_MpiComm *>(&communicator);
1135 if (mpi_comm != nullptr)
1136 return new Epetra_MpiComm(
1137 Utilities::MPI::duplicate_communicator(mpi_comm->GetMpiComm()));
1138# endif
1139
1140 // if we don't support MPI, or if the
1141 // communicator in question was in fact
1142 // not an MPI communicator, return a
1143 // copy of the same object again
1144 Assert(dynamic_cast<const Epetra_SerialComm *>(&communicator) != nullptr,
1146 return new Epetra_SerialComm(
1147 dynamic_cast<const Epetra_SerialComm &>(communicator));
1148 }
1149
1150
1151
1152 void
1153 destroy_communicator(Epetra_Comm &communicator)
1154 {
1155 // save the communicator, reset the map, and delete the communicator if
1156 // this whole thing was created as an MPI communicator
1157# ifdef DEAL_II_WITH_MPI
1158 Epetra_MpiComm *mpi_comm = dynamic_cast<Epetra_MpiComm *>(&communicator);
1159 if (mpi_comm != nullptr)
1160 {
1161 MPI_Comm comm = mpi_comm->GetMpiComm();
1162 *mpi_comm = Epetra_MpiComm(MPI_COMM_SELF);
1163
1165 }
1166# endif
1167 }
1168
1169
1170
1171 unsigned int
1172 get_n_mpi_processes(const Epetra_Comm &mpi_communicator)
1173 {
1174 return mpi_communicator.NumProc();
1175 }
1176
1177
1178 unsigned int
1179 get_this_mpi_process(const Epetra_Comm &mpi_communicator)
1180 {
1181 return static_cast<unsigned int>(mpi_communicator.MyPID());
1182 }
1183
1184
1185
1186 Epetra_Map
1187 duplicate_map(const Epetra_BlockMap &map, const Epetra_Comm &comm)
1188 {
1189 if (map.LinearMap() == true)
1190 {
1191 // each processor stores a
1192 // contiguous range of
1193 // elements in the
1194 // following constructor
1195 // call
1196 return Epetra_Map(map.NumGlobalElements(),
1197 map.NumMyElements(),
1198 map.IndexBase(),
1199 comm);
1200 }
1201 else
1202 {
1203 // the range is not
1204 // contiguous
1205 return Epetra_Map(map.NumGlobalElements(),
1206 map.NumMyElements(),
1207 map.MyGlobalElements(),
1208 0,
1209 comm);
1210 }
1211 }
1212 } // namespace Trilinos
1213
1214#endif
1215
1216#ifndef DOXYGEN
1217 template std::string
1218 to_string<int>(int, unsigned int);
1219 template std::string
1220 to_string<long int>(long int, unsigned int);
1221 template std::string
1222 to_string<long long int>(long long int, unsigned int);
1223 template std::string
1224 to_string<unsigned int>(unsigned int, unsigned int);
1225 template std::string
1226 to_string<unsigned long int>(unsigned long int, unsigned int);
1227 template std::string
1228 to_string<unsigned long long int>(unsigned long long int, unsigned int);
1229 template std::string
1230 to_string<float>(float, unsigned int);
1231 template std::string
1232 to_string<double>(double, unsigned int);
1233 template std::string
1234 to_string<long double>(long double, unsigned int);
1235
1236 template double
1237 truncate_to_n_digits(const double, const unsigned int);
1238 template float
1239 truncate_to_n_digits(const float, const unsigned int);
1240
1241 template std::vector<std::array<std::uint64_t, 1>>
1242 inverse_Hilbert_space_filling_curve<1, double>(
1243 const std::vector<Point<1, double>> &,
1244 const int);
1245 template std::vector<std::array<std::uint64_t, 1>>
1246 inverse_Hilbert_space_filling_curve<1>(
1247 const std::vector<std::array<std::uint64_t, 1>> &,
1248 const int);
1249 template std::vector<std::array<std::uint64_t, 2>>
1250 inverse_Hilbert_space_filling_curve<2, double>(
1251 const std::vector<Point<2, double>> &,
1252 const int);
1253 template std::vector<std::array<std::uint64_t, 2>>
1254 inverse_Hilbert_space_filling_curve<2>(
1255 const std::vector<std::array<std::uint64_t, 2>> &,
1256 const int);
1257 template std::vector<std::array<std::uint64_t, 3>>
1258 inverse_Hilbert_space_filling_curve<3, double>(
1259 const std::vector<Point<3, double>> &,
1260 const int);
1261 template std::vector<std::array<std::uint64_t, 3>>
1262 inverse_Hilbert_space_filling_curve<3>(
1263 const std::vector<std::array<std::uint64_t, 3>> &,
1264 const int);
1265
1266 template std::uint64_t
1267 pack_integers<1>(const std::array<std::uint64_t, 1> &, const int);
1268 template std::uint64_t
1269 pack_integers<2>(const std::array<std::uint64_t, 2> &, const int);
1270 template std::uint64_t
1271 pack_integers<3>(const std::array<std::uint64_t, 3> &, const int);
1272
1273#endif
1274
1275} // namespace Utilities
1276
A class that provides a separate storage location on each thread that accesses the object.
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_PACKAGE_VERSION
Definition: config.h:26
#define DEAL_II_VECTORIZATION_WIDTH_IN_BITS
Definition: config.h:127
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
#define DEAL_II_PACKAGE_NAME
Definition: config.h:24
static ::ExceptionBase & ExcInvalidNumber2StringConversersion(unsigned int arg1, unsigned int arg2)
static ::ExceptionBase & ExcOutOfMemory(std::size_t arg1)
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcInvalidNumber(unsigned int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcCantConvertString(std::string arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Expression fabs(const Expression &x)
Expression ceil(const Expression &x)
Expression floor(const Expression &x)
Expression log10(const Expression &x)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2050
static const char L
static const char N
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void malloc(T *&pointer, const unsigned int n_elements)
Definition: cuda.h:85
void free_communicator(MPI_Comm &mpi_communicator)
Definition: mpi.cc:194
bool job_supports_mpi()
Definition: mpi.cc:1014
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:183
void get_memory_stats(MemoryStats &stats)
Definition: utilities.cc:967
std::string get_hostname()
Definition: utilities.cc:1001
std::string get_time()
Definition: utilities.cc:1016
void posix_memalign(void **memptr, std::size_t alignment, std::size_t size)
Definition: utilities.cc:1047
double get_cpu_load()
Definition: utilities.cc:933
bool job_supports_mpi()
Definition: utilities.cc:1066
const std::string get_current_vectorization_level()
Definition: utilities.cc:941
std::string get_date()
Definition: utilities.cc:1032
void destroy_communicator(Epetra_Comm &communicator)
Definition: utilities.cc:1153
Epetra_Map duplicate_map(const Epetra_BlockMap &map, const Epetra_Comm &comm)
Definition: utilities.cc:1187
unsigned int get_this_mpi_process(const Epetra_Comm &mpi_communicator)
Definition: utilities.cc:1179
unsigned int get_n_mpi_processes(const Epetra_Comm &mpi_communicator)
Definition: utilities.cc:1172
const Epetra_Comm & comm_self()
Definition: utilities.cc:1110
const Teuchos::RCP< const Teuchos::Comm< int > > & tpetra_comm_self()
Definition: utilities.cc:1094
Epetra_Comm * duplicate_communicator(const Epetra_Comm &communicator)
Definition: utilities.cc:1126
const Epetra_Comm & comm_world()
Definition: utilities.cc:1078
std::vector< std::string > split_string_list(const std::string &s, const std::string &delimiter=",")
Definition: utilities.cc:704
Number truncate_to_n_digits(const Number number, const unsigned int n_digits)
Definition: utilities.cc:581
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:558
std::uint64_t pack_integers(const std::array< std::uint64_t, dim > &index, const int bits_per_dim)
Definition: utilities.cc:369
std::pair< int, unsigned int > get_integer_at_position(const std::string &name, const unsigned int position)
Definition: utilities.cc:850
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition: utilities.cc:436
std::string replace_in_string(const std::string &input, const std::string &from, const std::string &to)
Definition: utilities.cc:512
std::vector< unsigned char > decode_base64(const std::string &base64_input)
Definition: utilities.cc:449
std::vector< std::string > break_text_into_lines(const std::string &original_text, const unsigned int width, const char delimiter=' ')
Definition: utilities.cc:758
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
std::vector< std::array< std::uint64_t, dim > > inverse_Hilbert_space_filling_curve(const std::vector< Point< dim, Number > > &points, const int bits_per_dim=64)
Definition: utilities.cc:148
std::string compress(const std::string &input)
Definition: utilities.cc:392
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
bool match_at_string_start(const std::string &name, const std::string &pattern)
Definition: utilities.cc:835
std::string decompress(const std::string &compressed_input)
Definition: utilities.cc:414
unsigned int needed_digits(const unsigned int max_number)
Definition: utilities.cc:568
double string_to_double(const std::string &s)
Definition: utilities.cc:656
std::string dealii_version_string()
Definition: utilities.cc:97
std::string trim(const std::string &input)
Definition: utilities.cc:531
double generate_normal_random_number(const double a, const double sigma)
Definition: utilities.cc:892
int string_to_int(const std::string &s)
Definition: utilities.cc:608
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
Definition: types.h:201
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
unsigned long int VmRSS
Definition: utilities.h:923
unsigned long int VmPeak
Definition: utilities.h:907
unsigned long int VmSize
Definition: utilities.h:912
unsigned long int VmHWM
Definition: utilities.h:917
const MPI_Comm & comm