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