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