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