Reference documentation for deal.II version Git 4abc4a1666 2020-07-04 19:58:34 +0200
\(\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 - 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 #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>
29 #include <deal.II/base/utilities.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 
82 namespace 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  bio::gzip_params(boost::iostreams::gzip::default_compression)));
403  out.push(origin);
404  bio::copy(out, compressed);
405 
406  return compressed.str();
407 #else
408  return input;
409 #endif
410  }
411 
412 
413 
414  std::string
415  decompress(const std::string &compressed_input)
416  {
417 #ifdef DEAL_II_WITH_ZLIB
418  namespace bio = boost::iostreams;
419 
420  std::stringstream compressed(compressed_input);
421  std::stringstream decompressed;
422 
423  bio::filtering_streambuf<bio::input> out;
424  out.push(bio::gzip_decompressor());
425  out.push(compressed);
426  bio::copy(out, decompressed);
427 
428  return decompressed.str();
429 #else
430  return compressed_input;
431 #endif
432  }
433 
434 
435 
436  std::string
437  encode_base64(const std::vector<unsigned char> &binary_input)
438  {
439  using namespace boost::archive::iterators;
440  using It = base64_from_binary<
441  transform_width<std::vector<unsigned char>::const_iterator, 6, 8>>;
442  auto base64 = std::string(It(binary_input.begin()), It(binary_input.end()));
443  // Add padding.
444  return base64.append((3 - binary_input.size() % 3) % 3, '=');
445  }
446 
447 
448 
449  std::vector<unsigned char>
450  decode_base64(const std::string &base64_input)
451  {
452  using namespace boost::archive::iterators;
453  using It =
454  transform_width<binary_from_base64<std::string::const_iterator>, 8, 6>;
455  auto binary = std::vector<unsigned char>(It(base64_input.begin()),
456  It(base64_input.end()));
457  // Remove padding.
458  auto length = base64_input.size();
459  if (binary.size() > 2 && base64_input[length - 1] == '=' &&
460  base64_input[length - 2] == '=')
461  {
462  binary.erase(binary.end() - 2, binary.end());
463  }
464  else if (binary.size() > 1 && base64_input[length - 1] == '=')
465  {
466  binary.erase(binary.end() - 1, binary.end());
467  }
468  return binary;
469  }
470 
471 
472 
473  std::string
474  int_to_string(const unsigned int value, const unsigned int digits)
475  {
476  return to_string(value, digits);
477  }
478 
479 
480 
481  template <typename number>
482  std::string
483  to_string(const number value, const unsigned int digits)
484  {
485  // For integer data types, use the standard std::to_string()
486  // function. On the other hand, that function is defined in terms
487  // of std::sprintf, which does not use the usual std::iostream
488  // interface and tries to render floating point numbers in awkward
489  // ways (see
490  // https://en.cppreference.com/w/cpp/string/basic_string/to_string). So
491  // resort to boost::lexical_cast for all other types (in
492  // particular for floating point types.
493  std::string lc_string = (std::is_integral<number>::value ?
494  std::to_string(value) :
495  boost::lexical_cast<std::string>(value));
496 
497  if ((digits != numbers::invalid_unsigned_int) &&
498  (lc_string.size() < digits))
499  {
500  // We have to add the padding zeroes in front of the number
501  const unsigned int padding_position = (lc_string[0] == '-') ? 1 : 0;
502 
503  const std::string padding(digits - lc_string.size(), '0');
504  lc_string.insert(padding_position, padding);
505  }
506 
507  return lc_string;
508  }
509 
510 
511 
512  std::string
513  replace_in_string(const std::string &input,
514  const std::string &from,
515  const std::string &to)
516  {
517  if (from.empty())
518  return input;
519 
520  std::string out = input;
521  std::string::size_type pos = out.find(from);
522 
523  while (pos != std::string::npos)
524  {
525  out.replace(pos, from.size(), to);
526  pos = out.find(from, pos + to.size());
527  }
528  return out;
529  }
530 
531  std::string
532  trim(const std::string &input)
533  {
534  std::string::size_type left = 0;
535  std::string::size_type right = input.size() > 0 ? input.size() - 1 : 0;
536 
537  for (; left < input.size(); ++left)
538  {
539  if (!std::isspace(input[left]))
540  {
541  break;
542  }
543  }
544 
545  for (; right >= left; --right)
546  {
547  if (!std::isspace(input[right]))
548  {
549  break;
550  }
551  }
552 
553  return std::string(input, left, right - left + 1);
554  }
555 
556 
557 
558  std::string
559  dim_string(const int dim, const int spacedim)
560  {
561  if (dim == spacedim)
562  return int_to_string(dim);
563  else
564  return int_to_string(dim) + "," + int_to_string(spacedim);
565  }
566 
567 
568  unsigned int
569  needed_digits(const unsigned int max_number)
570  {
571  if (max_number > 0)
572  return static_cast<int>(
573  std::ceil(std::log10(std::fabs(max_number + 0.1))));
574 
575  return 1;
576  }
577 
578 
579 
580  template <typename Number>
581  Number
582  truncate_to_n_digits(const Number number, const unsigned int n_digits)
583  {
584  AssertThrow(n_digits >= 1, ExcMessage("invalid parameter."));
585 
586  if (!(std::fabs(number) > std::numeric_limits<Number>::min()))
587  return number;
588 
589  const int order =
590  static_cast<int>(std::floor(std::log10(std::fabs(number))));
591 
592  const int shift = -order + static_cast<int>(n_digits) - 1;
593 
594  Assert(shift <= static_cast<int>(std::floor(
596  ExcMessage(
597  "Overflow. Use a smaller value for n_digits and/or make sure "
598  "that the absolute value of 'number' does not become too small."));
599 
600  const Number factor = std::pow(10.0, static_cast<Number>(shift));
601 
602  const Number number_cutoff = std::trunc(number * factor) / factor;
603 
604  return number_cutoff;
605  }
606 
607 
608  int
609  string_to_int(const std::string &s_)
610  {
611  // trim whitespace on either side of the text if necessary
612  std::string s = s_;
613  while ((s.size() > 0) && (s[0] == ' '))
614  s.erase(s.begin());
615  while ((s.size() > 0) && (s[s.size() - 1] == ' '))
616  s.erase(s.end() - 1);
617 
618  // Now convert and see whether we succeed. Note that strtol only
619  // touches errno if an error occurred, so if we want to check
620  // whether an error happened, we need to make sure that errno==0
621  // before calling strtol since otherwise it may be that the
622  // conversion succeeds and that errno remains at the value it
623  // was before, whatever that was.
624  char *p;
625  errno = 0;
626  const int i = std::strtol(s.c_str(), &p, 10);
627 
628  // We have an error if one of the following conditions is true:
629  // - strtol sets errno != 0
630  // - The original string was empty (we could have checked that
631  // earlier already)
632  // - The string has non-zero length and strtol converted the
633  // first part to something useful, but stopped converting short
634  // of the terminating '\0' character. This happens, for example,
635  // if the given string is "1234 abc".
636  AssertThrow(!((errno != 0) || (s.size() == 0) ||
637  ((s.size() > 0) && (*p != '\0'))),
638  ExcMessage("Can't convert <" + s + "> to an integer."));
639 
640  return i;
641  }
642 
643 
644 
645  std::vector<int>
646  string_to_int(const std::vector<std::string> &s)
647  {
648  std::vector<int> tmp(s.size());
649  for (unsigned int i = 0; i < s.size(); ++i)
650  tmp[i] = string_to_int(s[i]);
651  return tmp;
652  }
653 
654 
655 
656  double
657  string_to_double(const std::string &s_)
658  {
659  // trim whitespace on either side of the text if necessary
660  std::string s = s_;
661  while ((s.size() > 0) && (s[0] == ' '))
662  s.erase(s.begin());
663  while ((s.size() > 0) && (s[s.size() - 1] == ' '))
664  s.erase(s.end() - 1);
665 
666  // Now convert and see whether we succeed. Note that strtol only
667  // touches errno if an error occurred, so if we want to check
668  // whether an error happened, we need to make sure that errno==0
669  // before calling strtol since otherwise it may be that the
670  // conversion succeeds and that errno remains at the value it
671  // was before, whatever that was.
672  char *p;
673  errno = 0;
674  const double d = std::strtod(s.c_str(), &p);
675 
676  // We have an error if one of the following conditions is true:
677  // - strtod sets errno != 0
678  // - The original string was empty (we could have checked that
679  // earlier already)
680  // - The string has non-zero length and strtod converted the
681  // first part to something useful, but stopped converting short
682  // of the terminating '\0' character. This happens, for example,
683  // if the given string is "1.234 abc".
684  AssertThrow(!((errno != 0) || (s.size() == 0) ||
685  ((s.size() > 0) && (*p != '\0'))),
686  ExcMessage("Can't convert <" + s + "> to a double."));
687 
688  return d;
689  }
690 
691 
692 
693  std::vector<double>
694  string_to_double(const std::vector<std::string> &s)
695  {
696  std::vector<double> tmp(s.size());
697  for (unsigned int i = 0; i < s.size(); ++i)
698  tmp[i] = string_to_double(s[i]);
699  return tmp;
700  }
701 
702 
703 
704  std::vector<std::string>
705  split_string_list(const std::string &s, const std::string &delimiter)
706  {
707  // keep the currently remaining part of the input string in 'tmp' and
708  // keep chopping elements of the list off the front
709  std::string tmp = s;
710 
711  // as discussed in the documentation, eat whitespace from the end
712  // of the string
713  while (tmp.length() != 0 && tmp[tmp.length() - 1] == ' ')
714  tmp.erase(tmp.length() - 1, 1);
715 
716  // split the input list until it is empty. since in every iteration
717  // 'tmp' is what's left of the string after the next delimiter,
718  // and since we've stripped trailing space already, 'tmp' will
719  // be empty at one point if 's' ended in a delimiter, even if
720  // there was space after the last delimiter. this matches what's
721  // discussed in the documentation
722  std::vector<std::string> split_list;
723  while (tmp.length() != 0)
724  {
725  std::string name;
726  name = tmp;
727 
728  if (name.find(delimiter) != std::string::npos)
729  {
730  name.erase(name.find(delimiter), std::string::npos);
731  tmp.erase(0, tmp.find(delimiter) + delimiter.size());
732  }
733  else
734  tmp = "";
735 
736  // strip spaces from this element's front and end
737  while ((name.length() != 0) && (name[0] == ' '))
738  name.erase(0, 1);
739  while (name.length() != 0 && name[name.length() - 1] == ' ')
740  name.erase(name.length() - 1, 1);
741 
742  split_list.push_back(name);
743  }
744 
745  return split_list;
746  }
747 
748 
749  std::vector<std::string>
750  split_string_list(const std::string &s, const char delimiter)
751  {
752  std::string d = ",";
753  d[0] = delimiter;
754  return split_string_list(s, d);
755  }
756 
757 
758  std::vector<std::string>
759  break_text_into_lines(const std::string &original_text,
760  const unsigned int width,
761  const char delimiter)
762  {
763  std::string text = original_text;
764  std::vector<std::string> lines;
765 
766  // remove trailing spaces
767  while ((text.length() != 0) && (text[text.length() - 1] == delimiter))
768  text.erase(text.length() - 1, 1);
769 
770  // then split the text into lines
771  while (text.length() != 0)
772  {
773  // in each iteration, first remove
774  // leading spaces
775  while ((text.length() != 0) && (text[0] == delimiter))
776  text.erase(0, 1);
777 
778  std::size_t pos_newline = text.find_first_of('\n', 0);
779  if (pos_newline != std::string::npos && pos_newline <= width)
780  {
781  std::string line(text, 0, pos_newline);
782  while ((line.length() != 0) &&
783  (line[line.length() - 1] == delimiter))
784  line.erase(line.length() - 1, 1);
785  lines.push_back(line);
786  text.erase(0, pos_newline + 1);
787  continue;
788  }
789 
790  // if we can fit everything into one
791  // line, then do so. otherwise, we have
792  // to keep breaking
793  if (text.length() < width)
794  {
795  // remove trailing spaces
796  while ((text.length() != 0) &&
797  (text[text.length() - 1] == delimiter))
798  text.erase(text.length() - 1, 1);
799  lines.push_back(text);
800  text = "";
801  }
802  else
803  {
804  // starting at position width, find the
805  // location of the previous space, so
806  // that we can break around there
807  int location = std::min<int>(width, text.length() - 1);
808  for (; location > 0; --location)
809  if (text[location] == delimiter)
810  break;
811 
812  // if there are no spaces, then try if
813  // there are spaces coming up
814  if (location == 0)
815  for (location = std::min<int>(width, text.length() - 1);
816  location < static_cast<int>(text.length());
817  ++location)
818  if (text[location] == delimiter)
819  break;
820 
821  // now take the text up to the found
822  // location and put it into a single
823  // line, and remove it from 'text'
824  std::string line(text, 0, location);
825  while ((line.length() != 0) &&
826  (line[line.length() - 1] == delimiter))
827  line.erase(line.length() - 1, 1);
828  lines.push_back(line);
829  text.erase(0, location);
830  }
831  }
832 
833  return lines;
834  }
835 
836 
837 
838  bool
839  match_at_string_start(const std::string &name, const std::string &pattern)
840  {
841  if (pattern.size() > name.size())
842  return false;
843 
844  for (unsigned int i = 0; i < pattern.size(); ++i)
845  if (pattern[i] != name[i])
846  return false;
847 
848  return true;
849  }
850 
851 
852 
853  std::pair<int, unsigned int>
854  get_integer_at_position(const std::string &name, const unsigned int position)
855  {
856  Assert(position < name.size(), ExcInternalError());
857 
858  const std::string test_string(name.begin() + position, name.end());
859 
860  std::istringstream str(test_string);
861 
862  int i;
863  if (str >> i)
864  {
865  // compute the number of
866  // digits of i. assuming it
867  // is less than 8 is likely
868  // ok
869  if (i < 10)
870  return std::make_pair(i, 1U);
871  else if (i < 100)
872  return std::make_pair(i, 2U);
873  else if (i < 1000)
874  return std::make_pair(i, 3U);
875  else if (i < 10000)
876  return std::make_pair(i, 4U);
877  else if (i < 100000)
878  return std::make_pair(i, 5U);
879  else if (i < 1000000)
880  return std::make_pair(i, 6U);
881  else if (i < 10000000)
882  return std::make_pair(i, 7U);
883  else
884  {
885  Assert(false, ExcNotImplemented());
886  return std::make_pair(-1, numbers::invalid_unsigned_int);
887  }
888  }
889  else
890  return std::make_pair(-1, numbers::invalid_unsigned_int);
891  }
892 
893 
894 
895  double
896  generate_normal_random_number(const double a, const double sigma)
897  {
898  // if no noise: return now
899  if (sigma == 0)
900  return a;
901 
902  // we would want to use rand(), but that function is not reentrant
903  // in a thread context. one could use rand_r, but this does not
904  // produce reproducible results between threads either (though at
905  // least it is reentrant). these two approaches being
906  // non-workable, use a thread-local random number generator here.
907  // we could use std::mt19937 but doing so results in compiler-dependent
908  // output.
909  static Threads::ThreadLocalStorage<boost::mt19937> random_number_generator;
910  return boost::normal_distribution<>(a,
911  sigma)(random_number_generator.get());
912  }
913 
914 
915 
916  namespace System
917  {
918 #if defined(__linux__)
919 
920  double
921  get_cpu_load()
922  {
923  std::ifstream cpuinfo;
924  cpuinfo.open("/proc/loadavg");
925 
926  AssertThrow(cpuinfo, ExcIO());
927 
928  double load;
929  cpuinfo >> load;
930 
931  return load;
932  }
933 
934 #else
935 
936  double
938  {
939  return 0.;
940  }
941 
942 #endif
943 
944  const std::string
946  {
948  {
949  case 0:
950  return "disabled";
951  case 128:
952 #ifdef __ALTIVEC__
953  return "AltiVec";
954 #else
955  return "SSE2";
956 #endif
957  case 256:
958  return "AVX";
959  case 512:
960  return "AVX512";
961  default:
962  AssertThrow(false,
964  "Invalid DEAL_II_VECTORIZATION_WIDTH_IN_BITS."));
965  return "ERROR";
966  }
967  }
968 
969 
970  void
972  {
973  stats.VmPeak = stats.VmSize = stats.VmHWM = stats.VmRSS = 0;
974 
975  // parsing /proc/self/stat would be a
976  // lot easier, but it does not contain
977  // VmHWM, so we use /status instead.
978 #if defined(__linux__)
979  std::ifstream file("/proc/self/status");
980  std::string line;
981  std::string name;
982  while (!file.eof())
983  {
984  file >> name;
985  if (name == "VmPeak:")
986  file >> stats.VmPeak;
987  else if (name == "VmSize:")
988  file >> stats.VmSize;
989  else if (name == "VmHWM:")
990  file >> stats.VmHWM;
991  else if (name == "VmRSS:")
992  {
993  file >> stats.VmRSS;
994  break; // this is always the last entry
995  }
996 
997  getline(file, line);
998  }
999 #endif
1000  }
1001 
1002 
1003 
1004  std::string
1006  {
1007 #if defined(DEAL_II_HAVE_UNISTD_H) && defined(DEAL_II_HAVE_GETHOSTNAME)
1008  const unsigned int N = 1024;
1009  char hostname[N];
1010  gethostname(&(hostname[0]), N - 1);
1011 #else
1012  std::string hostname("unknown");
1013 #endif
1014  return hostname;
1015  }
1016 
1017 
1018 
1019  std::string
1021  {
1022  std::time_t time1 = std::time(nullptr);
1023  std::tm * time = std::localtime(&time1);
1024 
1025  std::ostringstream o;
1026  o << time->tm_hour << ":" << (time->tm_min < 10 ? "0" : "")
1027  << time->tm_min << ":" << (time->tm_sec < 10 ? "0" : "")
1028  << time->tm_sec;
1029 
1030  return o.str();
1031  }
1032 
1033 
1034 
1035  std::string
1037  {
1038  std::time_t time1 = std::time(nullptr);
1039  std::tm * time = std::localtime(&time1);
1040 
1041  std::ostringstream o;
1042  o << time->tm_year + 1900 << "/" << time->tm_mon + 1 << "/"
1043  << time->tm_mday;
1044 
1045  return o.str();
1046  }
1047 
1048 
1049 
1050  void
1051  posix_memalign(void **memptr, std::size_t alignment, std::size_t size)
1052  {
1053 #ifndef DEAL_II_MSVC
1054  const int ierr = ::posix_memalign(memptr, alignment, size);
1055 
1056  AssertThrow(ierr == 0, ExcOutOfMemory(size));
1057  AssertThrow(*memptr != nullptr, ExcOutOfMemory(size));
1058 #else
1059  // Windows does not appear to have posix_memalign. just use the
1060  // regular malloc in that case
1061  *memptr = malloc(size);
1062  (void)alignment;
1063  AssertThrow(*memptr != 0, ExcOutOfMemory(size));
1064 #endif
1065  }
1066 
1067 
1068 
1069  bool
1071  {
1073  }
1074  } // namespace System
1075 
1076 
1077 #ifdef DEAL_II_WITH_TRILINOS
1078 
1079  namespace Trilinos
1080  {
1081  const Epetra_Comm &
1083  {
1084 # ifdef DEAL_II_WITH_MPI
1085  static Teuchos::RCP<Epetra_MpiComm> communicator =
1086  Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD), true);
1087 # else
1088  static Teuchos::RCP<Epetra_SerialComm> communicator =
1089  Teuchos::rcp(new Epetra_SerialComm(), true);
1090 # endif
1091 
1092  return *communicator;
1093  }
1094 
1095 
1096 
1097  const Teuchos::RCP<const Teuchos::Comm<int>> &
1099  {
1100 # ifdef DEAL_II_WITH_MPI
1101  static auto communicator = Teuchos::RCP<const Teuchos::Comm<int>>(
1102  new Teuchos::MpiComm<int>(MPI_COMM_SELF));
1103 # else
1104  static auto communicator =
1105  Teuchos::RCP<const Teuchos::Comm<int>>(new Teuchos::Comm<int>());
1106 # endif
1107 
1108  return communicator;
1109  }
1110 
1111 
1112 
1113  const Epetra_Comm &
1115  {
1116 # ifdef DEAL_II_WITH_MPI
1117  static Teuchos::RCP<Epetra_MpiComm> communicator =
1118  Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_SELF), true);
1119 # else
1120  static Teuchos::RCP<Epetra_SerialComm> communicator =
1121  Teuchos::rcp(new Epetra_SerialComm(), true);
1122 # endif
1123 
1124  return *communicator;
1125  }
1126 
1127 
1128 
1129  Epetra_Comm *
1130  duplicate_communicator(const Epetra_Comm &communicator)
1131  {
1132 # ifdef DEAL_II_WITH_MPI
1133 
1134  // see if the communicator is in fact a
1135  // parallel MPI communicator; if so,
1136  // return a duplicate of it
1137  const Epetra_MpiComm *mpi_comm =
1138  dynamic_cast<const Epetra_MpiComm *>(&communicator);
1139  if (mpi_comm != nullptr)
1140  return new Epetra_MpiComm(
1141  Utilities::MPI::duplicate_communicator(mpi_comm->GetMpiComm()));
1142 # endif
1143 
1144  // if we don't support MPI, or if the
1145  // communicator in question was in fact
1146  // not an MPI communicator, return a
1147  // copy of the same object again
1148  Assert(dynamic_cast<const Epetra_SerialComm *>(&communicator) != nullptr,
1149  ExcInternalError());
1150  return new Epetra_SerialComm(
1151  dynamic_cast<const Epetra_SerialComm &>(communicator));
1152  }
1153 
1154 
1155 
1156  void
1157  destroy_communicator(Epetra_Comm &communicator)
1158  {
1159  // save the communicator, reset the map, and delete the communicator if
1160  // this whole thing was created as an MPI communicator
1161 # ifdef DEAL_II_WITH_MPI
1162  Epetra_MpiComm *mpi_comm = dynamic_cast<Epetra_MpiComm *>(&communicator);
1163  if (mpi_comm != nullptr)
1164  {
1165  MPI_Comm comm = mpi_comm->GetMpiComm();
1166  *mpi_comm = Epetra_MpiComm(MPI_COMM_SELF);
1167  const int ierr = MPI_Comm_free(&comm);
1168  AssertThrowMPI(ierr);
1169  }
1170 # endif
1171  }
1172 
1173 
1174 
1175  unsigned int
1176  get_n_mpi_processes(const Epetra_Comm &mpi_communicator)
1177  {
1178  return mpi_communicator.NumProc();
1179  }
1180 
1181 
1182  unsigned int
1183  get_this_mpi_process(const Epetra_Comm &mpi_communicator)
1184  {
1185  return static_cast<unsigned int>(mpi_communicator.MyPID());
1186  }
1187 
1188 
1189 
1190  Epetra_Map
1191  duplicate_map(const Epetra_BlockMap &map, const Epetra_Comm &comm)
1192  {
1193  if (map.LinearMap() == true)
1194  {
1195  // each processor stores a
1196  // contiguous range of
1197  // elements in the
1198  // following constructor
1199  // call
1200  return Epetra_Map(map.NumGlobalElements(),
1201  map.NumMyElements(),
1202  map.IndexBase(),
1203  comm);
1204  }
1205  else
1206  {
1207  // the range is not
1208  // contiguous
1209  return Epetra_Map(map.NumGlobalElements(),
1210  map.NumMyElements(),
1211  map.MyGlobalElements(),
1212  0,
1213  comm);
1214  }
1215  }
1216  } // namespace Trilinos
1217 
1218 #endif
1219 
1220  template std::string
1221  to_string<int>(int, unsigned int);
1222  template std::string
1223  to_string<long int>(long int, unsigned int);
1224  template std::string
1225  to_string<long long int>(long long int, unsigned int);
1226  template std::string
1227  to_string<unsigned int>(unsigned int, unsigned int);
1228  template std::string
1229  to_string<unsigned long int>(unsigned long int, unsigned int);
1230  template std::string
1231  to_string<unsigned long long int>(unsigned long long int, unsigned int);
1232  template std::string
1233  to_string<float>(float, unsigned int);
1234  template std::string
1235  to_string<double>(double, unsigned int);
1236  template std::string
1237  to_string<long double>(long double, unsigned int);
1238 
1239  template double
1240  truncate_to_n_digits(const double, const unsigned int);
1241  template float
1242  truncate_to_n_digits(const float, const unsigned int);
1243 
1244  template std::vector<std::array<std::uint64_t, 1>>
1246  const std::vector<Point<1, double>> &,
1247  const int);
1248  template std::vector<std::array<std::uint64_t, 1>>
1250  const std::vector<std::array<std::uint64_t, 1>> &,
1251  const int);
1252  template std::vector<std::array<std::uint64_t, 2>>
1254  const std::vector<Point<2, double>> &,
1255  const int);
1256  template std::vector<std::array<std::uint64_t, 2>>
1258  const std::vector<std::array<std::uint64_t, 2>> &,
1259  const int);
1260  template std::vector<std::array<std::uint64_t, 3>>
1262  const std::vector<Point<3, double>> &,
1263  const int);
1264  template std::vector<std::array<std::uint64_t, 3>>
1266  const std::vector<std::array<std::uint64_t, 3>> &,
1267  const int);
1268 
1269  template std::uint64_t
1270  pack_integers<1>(const std::array<std::uint64_t, 1> &, const int);
1271  template std::uint64_t
1272  pack_integers<2>(const std::array<std::uint64_t, 2> &, const int);
1273  template std::uint64_t
1274  pack_integers<3>(const std::array<std::uint64_t, 3> &, const int);
1275 
1276 
1277 } // namespace Utilities
1278 
template std::vector< std::array< std::uint64_t, 2 > > inverse_Hilbert_space_filling_curve< 2, double >(const std::vector< Point< 2, double >> &, const int)
std::vector< std::string > split_string_list(const std::string &s, const std::string &delimiter=",")
Definition: utilities.cc:705
std::string encode_base64(const std::vector< unsigned char > &binary_input)
Definition: utilities.cc:437
void posix_memalign(void **memptr, std::size_t alignment, std::size_t size)
Definition: utilities.cc:1051
static const unsigned int invalid_unsigned_int
Definition: types.h:191
#define DEAL_II_VECTORIZATION_WIDTH_IN_BITS
Definition: config.h:119
Number truncate_to_n_digits(const Number number, const unsigned int n_digits)
Definition: utilities.cc:582
std::uint64_t pack_integers(const std::array< std::uint64_t, dim > &index, const int bits_per_dim)
Definition: utilities.cc:369
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
static const char L
std::vector< unsigned char > decode_base64(const std::string &base64_input)
Definition: utilities.cc:450
A class that provides a separate storage location on each thread that accesses the object...
types::global_dof_index size_type
Definition: cuda_kernels.h:45
Expression log10(const Expression &x)
static ::ExceptionBase & ExcIO()
std::string dealii_version_string()
Definition: utilities.cc:97
template std::string to_string< float >(float, unsigned int)
unsigned int get_n_mpi_processes(const Epetra_Comm &mpi_communicator)
Definition: utilities.cc:1176
std::string trim(const std::string &input)
Definition: utilities.cc:532
template std::vector< std::array< std::uint64_t, 1 > > inverse_Hilbert_space_filling_curve< 1, double >(const std::vector< Point< 1, double >> &, const int)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1628
std::pair< int, unsigned int > get_integer_at_position(const std::string &name, const unsigned int position)
Definition: utilities.cc:854
std::vector< std::string > break_text_into_lines(const std::string &original_text, const unsigned int width, const char delimiter=' ')
Definition: utilities.cc:759
template std::string to_string< long long int >(long long int, unsigned int)
template std::vector< std::array< std::uint64_t, 3 > > inverse_Hilbert_space_filling_curve< 3 >(const std::vector< std::array< std::uint64_t, 3 >> &, const int)
template std::uint64_t pack_integers< 2 >(const std::array< std::uint64_t, 2 > &, const int)
template std::string to_string< long double >(long double, unsigned int)
static const char U
template std::string to_string< int >(int, unsigned int)
template std::uint64_t pack_integers< 3 >(const std::array< std::uint64_t, 3 > &, const int)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
std::string decompress(const std::string &compressed_input)
Definition: utilities.cc:415
const Epetra_Comm & comm_self()
Definition: utilities.cc:1114
template int min< int >(const int &, const MPI_Comm &)
Definition: point.h:110
std::string get_date()
Definition: utilities.cc:1036
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:483
void get_memory_stats(MemoryStats &stats)
Definition: utilities.cc:971
#define DEAL_II_PACKAGE_NAME
Definition: config.h:24
void malloc(T *&pointer, const unsigned int n_elements)
Definition: cuda.h:84
double string_to_double(const std::string &s)
Definition: utilities.cc:657
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:376
template std::vector< std::array< std::uint64_t, 2 > > inverse_Hilbert_space_filling_curve< 2 >(const std::vector< std::array< std::uint64_t, 2 >> &, const int)
double get_cpu_load()
Definition: utilities.cc:937
const Teuchos::RCP< const Teuchos::Comm< int > > & tpetra_comm_self()
Definition: utilities.cc:1098
template std::string to_string< double >(double, unsigned int)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_PACKAGE_VERSION
Definition: config.h:26
std::string compress(const std::string &input)
Definition: utilities.cc:392
Expression ceil(const Expression &x)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
template std::vector< std::array< std::uint64_t, 1 > > inverse_Hilbert_space_filling_curve< 1 >(const std::vector< std::array< std::uint64_t, 1 >> &, const int)
double generate_normal_random_number(const double a, const double sigma)
Definition: utilities.cc:896
#define Assert(cond, exc)
Definition: exceptions.h:1403
static ::ExceptionBase & ExcOutOfMemory(std::size_t arg1)
template std::string to_string< long int >(long int, unsigned int)
unsigned long int VmSize
Definition: utilities.h:834
static ::ExceptionBase & ExcCantConvertString(std::string arg1)
static ::ExceptionBase & ExcInvalidNumber(unsigned int arg1)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
bool match_at_string_start(const std::string &name, const std::string &pattern)
Definition: utilities.cc:839
std::string to_string(const T &t)
Definition: patterns.h:2341
template std::string to_string< unsigned long int >(unsigned long int, unsigned int)
template std::uint64_t pack_integers< 1 >(const std::array< std::uint64_t, 1 > &, const int)
Expression fabs(const Expression &x)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
std::string replace_in_string(const std::string &input, const std::string &from, const std::string &to)
Definition: utilities.cc:513
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:559
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
void destroy_communicator(Epetra_Comm &communicator)
Definition: utilities.cc:1157
Definition: cuda.h:31
template std::string to_string< unsigned long long int >(unsigned long long int, unsigned int)
*braid_SplitCommworld & comm
std::string get_hostname()
Definition: utilities.cc:1005
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1684
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:139
unsigned long int VmHWM
Definition: utilities.h:839
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:413
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
T min(const T &t, const MPI_Comm &mpi_communicator)
const Epetra_Comm & comm_world()
Definition: utilities.cc:1082
static const char N
int string_to_int(const std::string &s)
Definition: utilities.cc:609
bool job_supports_mpi()
Definition: utilities.cc:1070
static ::ExceptionBase & ExcNotImplemented()
template std::string to_string< unsigned int >(unsigned int, unsigned int)
unsigned long int VmRSS
Definition: utilities.h:845
unsigned int get_this_mpi_process(const Epetra_Comm &mpi_communicator)
Definition: utilities.cc:1183
const std::string get_current_vectorization_level()
Definition: utilities.cc:945
Epetra_Map duplicate_map(const Epetra_BlockMap &map, const Epetra_Comm &comm)
Definition: utilities.cc:1191
std::string get_time()
Definition: utilities.cc:1020
bool job_supports_mpi()
Definition: mpi.cc:1030
template std::vector< std::array< std::uint64_t, 3 > > inverse_Hilbert_space_filling_curve< 3, double >(const std::vector< Point< 3, double >> &, const int)
void copy(const T *begin, const T *end, U *dest)
T max(const T &t, const MPI_Comm &mpi_communicator)
unsigned long int VmPeak
Definition: utilities.h:829
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:819
unsigned int needed_digits(const unsigned int max_number)
Definition: utilities.cc:569
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidNumber2StringConversersion(unsigned int arg1, unsigned int arg2)
Expression floor(const Expression &x)