Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
data_out_base.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2023 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 
19 #include <deal.II/base/mpi.h>
23 #include <deal.II/base/utilities.h>
24 
26 
27 #include <algorithm>
28 #include <cmath>
29 #include <cstdint>
30 #include <cstring>
31 #include <ctime>
32 #include <fstream>
33 #include <iomanip>
34 #include <limits>
35 #include <memory>
36 #include <set>
37 #include <sstream>
38 #include <vector>
39 
40 #ifdef DEAL_II_WITH_ZLIB
41 # include <zlib.h>
42 #endif
43 
44 #ifdef DEAL_II_WITH_HDF5
45 # include <hdf5.h>
46 #endif
47 
48 #include <boost/iostreams/copy.hpp>
49 #include <boost/iostreams/device/back_inserter.hpp>
50 #include <boost/iostreams/filtering_stream.hpp>
51 #ifdef DEAL_II_WITH_ZLIB
52 # include <boost/iostreams/filter/zlib.hpp>
53 #endif
54 
55 
56 
58 
59 // we need the following exception from a global function, so can't declare it
60 // in the usual way inside a class
61 namespace
62 {
63  DeclException2(ExcUnexpectedInput,
64  std::string,
65  std::string,
66  << "Unexpected input: expected line\n <" << arg1
67  << ">\nbut got\n <" << arg2 << ">");
68 }
69 
70 
71 namespace
72 {
73 #ifdef DEAL_II_WITH_ZLIB
74  constexpr bool deal_ii_with_zlib = true;
75 #else
76  constexpr bool deal_ii_with_zlib = false;
77 #endif
78 
79 
80 #ifdef DEAL_II_WITH_ZLIB
85  int
86  get_zlib_compression_level(const DataOutBase::CompressionLevel level)
87  {
88  switch (level)
89  {
91  return Z_NO_COMPRESSION;
93  return Z_BEST_SPEED;
95  return Z_BEST_COMPRESSION;
97  return Z_DEFAULT_COMPRESSION;
98  default:
99  Assert(false, ExcNotImplemented());
100  return Z_NO_COMPRESSION;
101  }
102  }
103 
104 # ifdef DEAL_II_WITH_MPI
109  int
110  get_boost_zlib_compression_level(const DataOutBase::CompressionLevel level)
111  {
112  switch (level)
113  {
115  return boost::iostreams::zlib::no_compression;
117  return boost::iostreams::zlib::best_speed;
119  return boost::iostreams::zlib::best_compression;
121  return boost::iostreams::zlib::default_compression;
122  default:
123  Assert(false, ExcNotImplemented());
124  return boost::iostreams::zlib::no_compression;
125  }
126  }
127 # endif
128 #endif
129 
134  template <typename T>
135  std::string
136  compress_array(const std::vector<T> &data,
137  const DataOutBase::CompressionLevel compression_level)
138  {
139 #ifdef DEAL_II_WITH_ZLIB
140  if (data.size() != 0)
141  {
142  const std::size_t uncompressed_size = (data.size() * sizeof(T));
143 
144  // While zlib's compress2 uses unsigned long (which is 64bits
145  // on Linux), the vtu compression header stores the block size
146  // as an std::uint32_t (see below). While we could implement
147  // writing several smaller blocks, we haven't done that. Let's
148  // trigger an error for the user instead:
149  AssertThrow(uncompressed_size <=
152 
153  // allocate a buffer for compressing data and do so
154  auto compressed_data_length = compressBound(uncompressed_size);
155  AssertThrow(compressed_data_length <=
158 
159  std::vector<unsigned char> compressed_data(compressed_data_length);
160 
161  int err = compress2(&compressed_data[0],
162  &compressed_data_length,
163  reinterpret_cast<const Bytef *>(data.data()),
164  uncompressed_size,
165  get_zlib_compression_level(compression_level));
166  (void)err;
167  Assert(err == Z_OK, ExcInternalError());
168 
169  // Discard the unnecessary bytes
170  compressed_data.resize(compressed_data_length);
171 
172  // now encode the compression header
173  const std::uint32_t compression_header[4] = {
174  1, /* number of blocks */
175  static_cast<std::uint32_t>(uncompressed_size), /* size of block */
176  static_cast<std::uint32_t>(
177  uncompressed_size), /* size of last block */
178  static_cast<std::uint32_t>(
179  compressed_data_length)}; /* list of compressed sizes of blocks */
180 
181  const auto *const header_start =
182  reinterpret_cast<const unsigned char *>(&compression_header[0]);
183 
184  return (Utilities::encode_base64(
185  {header_start, header_start + 4 * sizeof(std::uint32_t)}) +
186  Utilities::encode_base64(compressed_data));
187  }
188  else
189  return {};
190 #else
191  (void)data;
192  (void)compression_level;
193  Assert(false,
194  ExcMessage("This function can only be called if cmake found "
195  "a working libz installation."));
196  return {};
197 #endif
198  }
199 
200 
201 
210  template <typename T>
211  std::string
212  vtu_stringize_array(const std::vector<T> &data,
213  const DataOutBase::CompressionLevel compression_level,
214  const int precision)
215  {
216  if (deal_ii_with_zlib &&
217  (compression_level != DataOutBase::CompressionLevel::plain_text))
218  {
219  // compress the data we have in memory
220  return compress_array(data, compression_level);
221  }
222  else
223  {
224  std::ostringstream stream;
225  stream.precision(precision);
226  for (const T &el : data)
227  stream << el << ' ';
228  return stream.str();
229  }
230  }
231 
232 
241  struct ParallelIntermediateHeader
242  {
243  std::uint64_t magic;
244  std::uint64_t version;
245  std::uint64_t compression;
246  std::uint64_t dimension;
247  std::uint64_t space_dimension;
248  std::uint64_t n_ranks;
249  std::uint64_t n_patches;
250  };
251 } // namespace
252 
253 
254 // some declarations of functions and locally used classes
255 namespace DataOutBase
256 {
257  namespace
258  {
264  class SvgCell
265  {
266  public:
267  // Center of the cell (three-dimensional)
269 
274 
279  float depth;
280 
285 
286  // Center of the cell (projected, two-dimensional)
288 
292  bool
293  operator<(const SvgCell &) const;
294  };
295 
296  bool
297  SvgCell::operator<(const SvgCell &e) const
298  {
299  // note the "wrong" order in which we sort the elements
300  return depth > e.depth;
301  }
302 
303 
304 
310  class EpsCell2d
311  {
312  public:
316  Point<2> vertices[4];
317 
322  float color_value;
323 
328  float depth;
329 
333  bool
334  operator<(const EpsCell2d &) const;
335  };
336 
337  bool
338  EpsCell2d::operator<(const EpsCell2d &e) const
339  {
340  // note the "wrong" order in which we sort the elements
341  return depth > e.depth;
342  }
343 
344 
345 
357  template <int dim, int spacedim, typename Number = double>
358  std::unique_ptr<Table<2, Number>>
359  create_global_data_table(const std::vector<Patch<dim, spacedim>> &patches)
360  {
361  // If there is nothing to write, just return
362  if (patches.empty())
363  return std::make_unique<Table<2, Number>>();
364 
365  // unlike in the main function, we don't have here the data_names field,
366  // so we initialize it with the number of data sets in the first patch.
367  // the equivalence of these two definitions is checked in the main
368  // function.
369 
370  // we have to take care, however, whether the points are appended to the
371  // end of the patch.data table
372  const unsigned int n_data_sets = patches[0].points_are_available ?
373  (patches[0].data.n_rows() - spacedim) :
374  patches[0].data.n_rows();
375  const unsigned int n_data_points =
376  std::accumulate(patches.begin(),
377  patches.end(),
378  0U,
379  [](const unsigned int count,
380  const Patch<dim, spacedim> &patch) {
381  return count + patch.data.n_cols();
382  });
383 
384  std::unique_ptr<Table<2, Number>> global_data_table =
385  std::make_unique<Table<2, Number>>(n_data_sets, n_data_points);
386 
387  // loop over all patches
388  unsigned int next_value = 0;
389  for (const auto &patch : patches)
390  {
391  const unsigned int n_subdivisions = patch.n_subdivisions;
392  (void)n_subdivisions;
393 
394  Assert((patch.data.n_rows() == n_data_sets &&
395  !patch.points_are_available) ||
396  (patch.data.n_rows() == n_data_sets + spacedim &&
397  patch.points_are_available),
398  ExcDimensionMismatch(patch.points_are_available ?
399  (n_data_sets + spacedim) :
400  n_data_sets,
401  patch.data.n_rows()));
402  Assert(patch.reference_cell != ReferenceCells::get_hypercube<dim>() ||
403  (n_data_sets == 0) ||
404  (patch.data.n_cols() ==
405  Utilities::fixed_power<dim>(n_subdivisions + 1)),
406  ExcInvalidDatasetSize(patch.data.n_cols(),
407  n_subdivisions + 1));
408 
409  for (unsigned int i = 0; i < patch.data.n_cols(); ++i, ++next_value)
410  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
411  (*global_data_table)[data_set][next_value] =
412  patch.data(data_set, i);
413  }
414  Assert(next_value == n_data_points, ExcInternalError());
415 
416  return global_data_table;
417  }
418  } // namespace
419 
420 
421 
423  : flags(false, true)
424  , node_dim(numbers::invalid_unsigned_int)
425  , num_cells(0)
426  {}
427 
428 
429 
431  : flags(flags)
432  , node_dim(numbers::invalid_unsigned_int)
433  , num_cells(0)
434  {}
435 
436 
437 
438  template <int dim>
439  void
440  DataOutFilter::write_point(const unsigned int index, const Point<dim> &p)
441  {
442  node_dim = dim;
443 
444  Point<3> int_pt;
445  for (unsigned int d = 0; d < dim; ++d)
446  int_pt(d) = p(d);
447 
448  const Map3DPoint::const_iterator it = existing_points.find(int_pt);
449  unsigned int internal_ind;
450 
451  // If the point isn't in the set, or we're not filtering duplicate points,
452  // add it
453  if (it == existing_points.end() || !flags.filter_duplicate_vertices)
454  {
455  internal_ind = existing_points.size();
456  existing_points.insert(std::make_pair(int_pt, internal_ind));
457  }
458  else
459  {
460  internal_ind = it->second;
461  }
462  // Now add the index to the list of filtered points
463  filtered_points[index] = internal_ind;
464  }
465 
466 
467 
468  void
470  const unsigned int pt_index)
471  {
473 
474  // (Re)-initialize counter at any first call to this method.
475  if (cell_index == 0)
476  num_cells = 1;
477  }
478 
479 
480 
481  void
482  DataOutFilter::fill_node_data(std::vector<double> &node_data) const
483  {
484  node_data.resize(existing_points.size() * node_dim);
485 
486  for (const auto &existing_point : existing_points)
487  {
488  for (unsigned int d = 0; d < node_dim; ++d)
489  node_data[node_dim * existing_point.second + d] =
490  existing_point.first(d);
491  }
492  }
493 
494 
495 
496  void
497  DataOutFilter::fill_cell_data(const unsigned int local_node_offset,
498  std::vector<unsigned int> &cell_data) const
499  {
500  cell_data.resize(filtered_cells.size());
501 
502  for (const auto &filtered_cell : filtered_cells)
503  {
504  cell_data[filtered_cell.first] =
505  filtered_cell.second + local_node_offset;
506  }
507  }
508 
509 
510 
511  std::string
512  DataOutFilter::get_data_set_name(const unsigned int set_num) const
513  {
514  return data_set_names.at(set_num);
515  }
516 
517 
518 
519  unsigned int
520  DataOutFilter::get_data_set_dim(const unsigned int set_num) const
521  {
522  return data_set_dims.at(set_num);
523  }
524 
525 
526 
527  const double *
528  DataOutFilter::get_data_set(const unsigned int set_num) const
529  {
530  return data_sets[set_num].data();
531  }
532 
533 
534 
535  unsigned int
537  {
538  return existing_points.size();
539  }
540 
541 
542 
543  unsigned int
545  {
546  return num_cells;
547  }
548 
549 
550 
551  unsigned int
553  {
554  return data_set_names.size();
555  }
556 
557 
558 
559  void
561  {}
562 
563 
564 
565  void
567  {}
568 
569 
570 
571  template <int dim>
572  void
573  DataOutFilter::write_cell(const unsigned int index,
574  const unsigned int start,
575  const std::array<unsigned int, dim> &offsets)
576  {
577  ++num_cells;
578 
579  const unsigned int base_entry =
581 
582  switch (dim)
583  {
584  case 0:
585  {
586  internal_add_cell(base_entry + 0, start);
587  break;
588  }
589 
590  case 1:
591  {
592  const unsigned int d1 = offsets[0];
593 
594  internal_add_cell(base_entry + 0, start);
595  internal_add_cell(base_entry + 1, start + d1);
596  break;
597  }
598 
599  case 2:
600  {
601  const unsigned int d1 = offsets[0];
602  const unsigned int d2 = offsets[1];
603 
604  internal_add_cell(base_entry + 0, start);
605  internal_add_cell(base_entry + 1, start + d1);
606  internal_add_cell(base_entry + 2, start + d2 + d1);
607  internal_add_cell(base_entry + 3, start + d2);
608  break;
609  }
610 
611  case 3:
612  {
613  const unsigned int d1 = offsets[0];
614  const unsigned int d2 = offsets[1];
615  const unsigned int d3 = offsets[2];
616 
617  internal_add_cell(base_entry + 0, start);
618  internal_add_cell(base_entry + 1, start + d1);
619  internal_add_cell(base_entry + 2, start + d2 + d1);
620  internal_add_cell(base_entry + 3, start + d2);
621  internal_add_cell(base_entry + 4, start + d3);
622  internal_add_cell(base_entry + 5, start + d3 + d1);
623  internal_add_cell(base_entry + 6, start + d3 + d2 + d1);
624  internal_add_cell(base_entry + 7, start + d3 + d2);
625  break;
626  }
627 
628  default:
629  Assert(false, ExcNotImplemented());
630  }
631  }
632 
633 
634 
635  void
636  DataOutFilter::write_cell_single(const unsigned int index,
637  const unsigned int start,
638  const unsigned int n_points,
640  {
641  ++num_cells;
642 
643  const unsigned int base_entry = index * n_points;
644 
645  static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
646 
647  for (unsigned int i = 0; i < n_points; ++i)
648  internal_add_cell(base_entry + i,
650  table[i] :
651  i));
652  }
653 
654 
655 
656  void
657  DataOutFilter::write_data_set(const std::string &name,
658  const unsigned int dimension,
659  const unsigned int set_num,
660  const Table<2, double> &data_vectors)
661  {
662  unsigned int new_dim;
663 
664  // HDF5/XDMF output only supports 1d or 3d output, so force rearrangement if
665  // needed
666  if (flags.xdmf_hdf5_output && dimension != 1)
667  new_dim = 3;
668  else
669  new_dim = dimension;
670 
671  // Record the data set name, dimension, and allocate space for it
672  data_set_names.push_back(name);
673  data_set_dims.push_back(new_dim);
674  data_sets.emplace_back(new_dim * existing_points.size());
675 
676  // TODO: averaging, min/max, etc for merged vertices
677  for (unsigned int i = 0; i < filtered_points.size(); ++i)
678  {
679  const unsigned int r = filtered_points[i];
680 
681  for (unsigned int d = 0; d < new_dim; ++d)
682  {
683  if (d < dimension)
684  data_sets.back()[r * new_dim + d] = data_vectors(set_num + d, i);
685  else
686  data_sets.back()[r * new_dim + d] = 0;
687  }
688  }
689  }
690 } // namespace DataOutBase
691 
692 
693 
694 //----------------------------------------------------------------------//
695 // Auxiliary data
696 //----------------------------------------------------------------------//
697 
698 namespace
699 {
700  const char *gmv_cell_type[4] = {"", "line 2", "quad 4", "hex 8"};
701 
702  const char *ucd_cell_type[4] = {"pt", "line", "quad", "hex"};
703 
704  const char *tecplot_cell_type[4] = {"", "lineseg", "quadrilateral", "brick"};
705 
719  template <int dim, int spacedim>
720  std::array<unsigned int, 3>
721  extract_vtk_patch_info(const DataOutBase::Patch<dim, spacedim> &patch,
722  const bool write_higher_order_cells)
723  {
724  std::array<unsigned int, 3> vtk_cell_id = {
725  {/* cell type, tbd: */ numbers::invalid_unsigned_int,
726  /* # of cells, default: just one cell */ 1,
727  /* # of nodes, default: as many nodes as vertices */
728  patch.reference_cell.n_vertices()}};
729 
730  if (write_higher_order_cells)
731  {
732  vtk_cell_id[0] = patch.reference_cell.vtk_lagrange_type();
733  vtk_cell_id[2] = patch.data.n_cols();
734  }
735  else if (patch.data.n_cols() == patch.reference_cell.n_vertices())
736  // One data set per vertex -> a linear cell
737  vtk_cell_id[0] = patch.reference_cell.vtk_linear_type();
738  else if (patch.reference_cell == ReferenceCells::Triangle &&
739  patch.data.n_cols() == 6)
740  {
741  Assert(patch.n_subdivisions == 2, ExcInternalError());
742  vtk_cell_id[0] = patch.reference_cell.vtk_quadratic_type();
743  vtk_cell_id[2] = patch.data.n_cols();
744  }
745  else if (patch.reference_cell == ReferenceCells::Tetrahedron &&
746  patch.data.n_cols() == 10)
747  {
748  Assert(patch.n_subdivisions == 2, ExcInternalError());
749  vtk_cell_id[0] = patch.reference_cell.vtk_quadratic_type();
750  vtk_cell_id[2] = patch.data.n_cols();
751  }
752  else if (patch.reference_cell.is_hyper_cube())
753  {
754  // For hypercubes, we support sub-divided linear cells
755  vtk_cell_id[0] = patch.reference_cell.vtk_linear_type();
756  vtk_cell_id[1] = Utilities::pow(patch.n_subdivisions, dim);
757  }
758  else
759  {
760  Assert(false, ExcNotImplemented());
761  }
762 
763  return vtk_cell_id;
764  }
765 
766  //----------------------------------------------------------------------//
767  // Auxiliary functions
768  //----------------------------------------------------------------------//
769 
770  // For a given patch that corresponds to a hypercube cell, compute the
771  // location of a node interpolating the corner nodes linearly
772  // at the point lattice_location/n_subdivisions where lattice_location
773  // is a dim-dimensional integer vector. If the points are
774  // saved in the patch.data member, return the saved point instead.
775  template <int dim, int spacedim>
776  inline Point<spacedim>
777  get_equispaced_location(
779  const std::initializer_list<unsigned int> &lattice_location,
780  const unsigned int n_subdivisions)
781  {
782  // This function only makes sense when called on hypercube cells
784 
785  Assert(lattice_location.size() == dim, ExcInternalError());
786 
787  const unsigned int xstep = (dim > 0 ? *(lattice_location.begin() + 0) : 0);
788  const unsigned int ystep = (dim > 1 ? *(lattice_location.begin() + 1) : 0);
789  const unsigned int zstep = (dim > 2 ? *(lattice_location.begin() + 2) : 0);
790 
791  // If the patch stores the locations of nodes (rather than of only the
792  // vertices), then obtain the location by direct lookup.
793  if (patch.points_are_available)
794  {
795  Assert(n_subdivisions == patch.n_subdivisions, ExcNotImplemented());
796 
797  unsigned int point_no = 0;
798  switch (dim)
799  {
800  case 3:
801  AssertIndexRange(zstep, n_subdivisions + 1);
802  point_no += (n_subdivisions + 1) * (n_subdivisions + 1) * zstep;
804  case 2:
805  AssertIndexRange(ystep, n_subdivisions + 1);
806  point_no += (n_subdivisions + 1) * ystep;
808  case 1:
809  AssertIndexRange(xstep, n_subdivisions + 1);
810  point_no += xstep;
812  case 0:
813  // break here for dim<=3
814  break;
815 
816  default:
817  Assert(false, ExcNotImplemented());
818  }
819  Point<spacedim> node;
820  for (unsigned int d = 0; d < spacedim; ++d)
821  node[d] = patch.data(patch.data.size(0) - spacedim + d, point_no);
822  return node;
823  }
824  else
825  // The patch does not store node locations, so we have to interpolate
826  // between its vertices:
827  {
828  if (dim == 0)
829  return patch.vertices[0];
830  else
831  {
832  // perform a dim-linear interpolation
833  const double stepsize = 1. / n_subdivisions;
834  const double xfrac = xstep * stepsize;
835 
836  Point<spacedim> node =
837  (patch.vertices[1] * xfrac) + (patch.vertices[0] * (1 - xfrac));
838  if (dim > 1)
839  {
840  const double yfrac = ystep * stepsize;
841  node *= 1 - yfrac;
842  node += ((patch.vertices[3] * xfrac) +
843  (patch.vertices[2] * (1 - xfrac))) *
844  yfrac;
845  if (dim > 2)
846  {
847  const double zfrac = zstep * stepsize;
848  node *= (1 - zfrac);
849  node += (((patch.vertices[5] * xfrac) +
850  (patch.vertices[4] * (1 - xfrac))) *
851  (1 - yfrac) +
852  ((patch.vertices[7] * xfrac) +
853  (patch.vertices[6] * (1 - xfrac))) *
854  yfrac) *
855  zfrac;
856  }
857  }
858  return node;
859  }
860  }
861  }
862 
863  // For a given patch, compute the nodes for arbitrary (non-hypercube) cells.
864  // If the points are saved in the patch.data member, return the saved point
865  // instead.
866  template <int dim, int spacedim>
867  inline Point<spacedim>
868  get_node_location(const DataOutBase::Patch<dim, spacedim> &patch,
869  const unsigned int node_index)
870  {
871  // Due to a historical accident, we are using a different indexing
872  // for pyramids in this file than we do where we create patches.
873  // So translate if necessary.
874  unsigned int point_no_actual = node_index;
876  {
878 
879  static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
880  point_no_actual = table[node_index];
881  }
882 
883  // If the patch stores the locations of nodes (rather than of only the
884  // vertices), then obtain the location by direct lookup.
885  if (patch.points_are_available)
886  {
887  Point<spacedim> node;
888  for (unsigned int d = 0; d < spacedim; ++d)
889  node[d] =
890  patch.data(patch.data.size(0) - spacedim + d, point_no_actual);
891  return node;
892  }
893  else
894  // The patch does not store node locations, so we have to interpolate
895  // between its vertices. This isn't currently implemented for anything
896  // other than one subdivision, but would go here.
897  //
898  // For n_subdivisions==1, the locations are simply those of vertices, so
899  // get the information from there.
900  {
902 
903  return patch.vertices[point_no_actual];
904  }
905  }
906 
907 
908 
914  template <int dim, int spacedim>
915  std::tuple<unsigned int, unsigned int>
916  count_nodes_and_cells(
917  const std::vector<DataOutBase::Patch<dim, spacedim>> &patches)
918  {
919  unsigned int n_nodes = 0;
920  unsigned int n_cells = 0;
921  for (const auto &patch : patches)
922  {
924  ExcMessage(
925  "The reference cell for this patch is set to 'Invalid', "
926  "but that is clearly not a valid choice. Did you forget "
927  "to set the reference cell for the patch?"));
928 
929  if (patch.reference_cell.is_hyper_cube())
930  {
931  n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
932  n_cells += Utilities::fixed_power<dim>(patch.n_subdivisions);
933  }
934  else
935  {
936  Assert(patch.n_subdivisions == 1, ExcNotImplemented());
937  n_nodes += patch.reference_cell.n_vertices();
938  n_cells += 1;
939  }
940  }
941 
942  return std::make_tuple(n_nodes, n_cells);
943  }
944 
945 
946 
952  template <int dim, int spacedim>
953  std::tuple<unsigned int, unsigned int, unsigned int>
954  count_nodes_and_cells_and_points(
955  const std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
956  const bool write_higher_order_cells)
957  {
958  unsigned int n_nodes = 0;
959  unsigned int n_cells = 0;
960  unsigned int n_points_and_n_cells = 0;
961 
962  for (const auto &patch : patches)
963  {
964  if (patch.reference_cell.is_hyper_cube())
965  {
966  n_nodes += Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
967 
968  if (write_higher_order_cells)
969  {
970  // Write all of these nodes as a single higher-order cell. So
971  // add one to the number of cells, and update the number of
972  // points appropriately.
973  n_cells += 1;
974  n_points_and_n_cells +=
975  1 + Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
976  }
977  else
978  {
979  // Write all of these nodes as a collection of d-linear
980  // cells. Add the number of sub-cells to the total number of
981  // cells, and then add one for each cell plus the number of
982  // vertices per cell for each subcell to the number of points.
983  const unsigned int n_subcells =
984  Utilities::fixed_power<dim>(patch.n_subdivisions);
985  n_cells += n_subcells;
986  n_points_and_n_cells +=
987  n_subcells * (1 + GeometryInfo<dim>::vertices_per_cell);
988  }
989  }
990  else
991  {
992  n_nodes += patch.data.n_cols();
993  n_cells += 1;
994  n_points_and_n_cells += patch.data.n_cols() + 1;
995  }
996  }
997 
998  return std::make_tuple(n_nodes, n_cells, n_points_and_n_cells);
999  }
1000 
1006  template <typename FlagsType>
1007  class StreamBase
1008  {
1009  public:
1010  /*
1011  * Constructor. Stores a reference to the output stream for immediate use.
1012  */
1013  StreamBase(std::ostream &stream, const FlagsType &flags)
1014  : selected_component(numbers::invalid_unsigned_int)
1015  , stream(stream)
1016  , flags(flags)
1017  {}
1018 
1023  template <int dim>
1024  void
1025  write_point(const unsigned int, const Point<dim> &)
1026  {
1027  Assert(false,
1028  ExcMessage("The derived class you are using needs to "
1029  "reimplement this function if you want to call "
1030  "it."));
1031  }
1032 
1038  void
1039  flush_points()
1040  {}
1041 
1047  template <int dim>
1048  void
1049  write_cell(const unsigned int /*index*/,
1050  const unsigned int /*start*/,
1051  std::array<unsigned int, dim> & /*offsets*/)
1052  {
1053  Assert(false,
1054  ExcMessage("The derived class you are using needs to "
1055  "reimplement this function if you want to call "
1056  "it."));
1057  }
1058 
1065  void
1066  write_cell_single(const unsigned int index,
1067  const unsigned int start,
1068  const unsigned int n_points,
1070  {
1071  (void)index;
1072  (void)start;
1073  (void)n_points;
1074  (void)reference_cell;
1075 
1076  Assert(false,
1077  ExcMessage("The derived class you are using needs to "
1078  "reimplement this function if you want to call "
1079  "it."));
1080  }
1081 
1088  void
1089  flush_cells()
1090  {}
1091 
1096  template <typename T>
1097  std::ostream &
1098  operator<<(const T &t)
1099  {
1100  stream << t;
1101  return stream;
1102  }
1103 
1110  unsigned int selected_component;
1111 
1112  protected:
1117  std::ostream &stream;
1118 
1122  const FlagsType flags;
1123  };
1124 
1128  class DXStream : public StreamBase<DataOutBase::DXFlags>
1129  {
1130  public:
1131  DXStream(std::ostream &stream, const DataOutBase::DXFlags &flags);
1132 
1133  template <int dim>
1134  void
1135  write_point(const unsigned int index, const Point<dim> &);
1136 
1145  template <int dim>
1146  void
1147  write_cell(const unsigned int index,
1148  const unsigned int start,
1149  const std::array<unsigned int, dim> &offsets);
1150 
1157  template <typename data>
1158  void
1159  write_dataset(const unsigned int index, const std::vector<data> &values);
1160  };
1161 
1165  class GmvStream : public StreamBase<DataOutBase::GmvFlags>
1166  {
1167  public:
1168  GmvStream(std::ostream &stream, const DataOutBase::GmvFlags &flags);
1169 
1170  template <int dim>
1171  void
1172  write_point(const unsigned int index, const Point<dim> &);
1173 
1182  template <int dim>
1183  void
1184  write_cell(const unsigned int index,
1185  const unsigned int start,
1186  const std::array<unsigned int, dim> &offsets);
1187  };
1188 
1192  class TecplotStream : public StreamBase<DataOutBase::TecplotFlags>
1193  {
1194  public:
1195  TecplotStream(std::ostream &stream, const DataOutBase::TecplotFlags &flags);
1196 
1197  template <int dim>
1198  void
1199  write_point(const unsigned int index, const Point<dim> &);
1200 
1209  template <int dim>
1210  void
1211  write_cell(const unsigned int index,
1212  const unsigned int start,
1213  const std::array<unsigned int, dim> &offsets);
1214  };
1215 
1219  class UcdStream : public StreamBase<DataOutBase::UcdFlags>
1220  {
1221  public:
1222  UcdStream(std::ostream &stream, const DataOutBase::UcdFlags &flags);
1223 
1224  template <int dim>
1225  void
1226  write_point(const unsigned int index, const Point<dim> &);
1227 
1238  template <int dim>
1239  void
1240  write_cell(const unsigned int index,
1241  const unsigned int start,
1242  const std::array<unsigned int, dim> &offsets);
1243 
1250  template <typename data>
1251  void
1252  write_dataset(const unsigned int index, const std::vector<data> &values);
1253  };
1254 
1258  class VtkStream : public StreamBase<DataOutBase::VtkFlags>
1259  {
1260  public:
1261  VtkStream(std::ostream &stream, const DataOutBase::VtkFlags &flags);
1262 
1263  template <int dim>
1264  void
1265  write_point(const unsigned int index, const Point<dim> &);
1266 
1275  template <int dim>
1276  void
1277  write_cell(const unsigned int index,
1278  const unsigned int start,
1279  const std::array<unsigned int, dim> &offsets);
1280 
1284  void
1285  write_cell_single(const unsigned int index,
1286  const unsigned int start,
1287  const unsigned int n_points,
1288  const ReferenceCell &reference_cell);
1289 
1297  template <int dim>
1298  void
1299  write_high_order_cell(const unsigned int start,
1300  const std::vector<unsigned> &connectivity);
1301  };
1302 
1303 
1304  //----------------------------------------------------------------------//
1305 
1306  DXStream::DXStream(std::ostream &out, const DataOutBase::DXFlags &f)
1307  : StreamBase<DataOutBase::DXFlags>(out, f)
1308  {}
1309 
1310 
1311  template <int dim>
1312  void
1313  DXStream::write_point(const unsigned int, const Point<dim> &p)
1314  {
1315  if (flags.coordinates_binary)
1316  {
1317  float data[dim];
1318  for (unsigned int d = 0; d < dim; ++d)
1319  data[d] = p(d);
1320  stream.write(reinterpret_cast<const char *>(data), dim * sizeof(*data));
1321  }
1322  else
1323  {
1324  for (unsigned int d = 0; d < dim; ++d)
1325  stream << p(d) << '\t';
1326  stream << '\n';
1327  }
1328  }
1329 
1330 
1331 
1332  // Separate these out to avoid an internal compiler error with intel 17
1334  {
1339  std::array<unsigned int, GeometryInfo<0>::vertices_per_cell>
1340  set_node_numbers(const unsigned int /*start*/,
1341  const std::array<unsigned int, 0> & /*d1*/)
1342  {
1343  Assert(false, ExcInternalError());
1344  return {};
1345  }
1346 
1347 
1348 
1349  std::array<unsigned int, GeometryInfo<1>::vertices_per_cell>
1350  set_node_numbers(const unsigned int start,
1351  const std::array<unsigned int, 1> &offsets)
1352  {
1353  std::array<unsigned int, GeometryInfo<1>::vertices_per_cell> nodes;
1354  nodes[0] = start;
1355  nodes[1] = start + offsets[0];
1356  return nodes;
1357  }
1358 
1359 
1360 
1361  std::array<unsigned int, GeometryInfo<2>::vertices_per_cell>
1362  set_node_numbers(const unsigned int start,
1363  const std::array<unsigned int, 2> &offsets)
1364 
1365  {
1366  const unsigned int d1 = offsets[0];
1367  const unsigned int d2 = offsets[1];
1368 
1369  std::array<unsigned int, GeometryInfo<2>::vertices_per_cell> nodes;
1370  nodes[0] = start;
1371  nodes[1] = start + d1;
1372  nodes[2] = start + d2;
1373  nodes[3] = start + d2 + d1;
1374  return nodes;
1375  }
1376 
1377 
1378 
1379  std::array<unsigned int, GeometryInfo<3>::vertices_per_cell>
1380  set_node_numbers(const unsigned int start,
1381  const std::array<unsigned int, 3> &offsets)
1382  {
1383  const unsigned int d1 = offsets[0];
1384  const unsigned int d2 = offsets[1];
1385  const unsigned int d3 = offsets[2];
1386 
1387  std::array<unsigned int, GeometryInfo<3>::vertices_per_cell> nodes;
1388  nodes[0] = start;
1389  nodes[1] = start + d1;
1390  nodes[2] = start + d2;
1391  nodes[3] = start + d2 + d1;
1392  nodes[4] = start + d3;
1393  nodes[5] = start + d3 + d1;
1394  nodes[6] = start + d3 + d2;
1395  nodes[7] = start + d3 + d2 + d1;
1396  return nodes;
1397  }
1398  } // namespace DataOutBaseImplementation
1399 
1400 
1401 
1402  template <int dim>
1403  void
1404  DXStream::write_cell(const unsigned int,
1405  const unsigned int start,
1406  const std::array<unsigned int, dim> &offsets)
1407  {
1408  const auto nodes =
1409  DataOutBaseImplementation::set_node_numbers(start, offsets);
1410 
1411  if (flags.int_binary)
1412  {
1413  std::array<unsigned int, GeometryInfo<dim>::vertices_per_cell> temp;
1414  for (unsigned int i = 0; i < nodes.size(); ++i)
1415  temp[i] = nodes[GeometryInfo<dim>::dx_to_deal[i]];
1416  stream.write(reinterpret_cast<const char *>(temp.data()),
1417  temp.size() * sizeof(temp[0]));
1418  }
1419  else
1420  {
1421  for (unsigned int i = 0; i < nodes.size() - 1; ++i)
1422  stream << nodes[GeometryInfo<dim>::dx_to_deal[i]] << '\t';
1423  stream << nodes[GeometryInfo<dim>::dx_to_deal[nodes.size() - 1]]
1424  << '\n';
1425  }
1426  }
1427 
1428 
1429 
1430  template <typename data>
1431  inline void
1432  DXStream::write_dataset(const unsigned int, const std::vector<data> &values)
1433  {
1434  if (flags.data_binary)
1435  {
1436  stream.write(reinterpret_cast<const char *>(values.data()),
1437  values.size() * sizeof(data));
1438  }
1439  else
1440  {
1441  for (unsigned int i = 0; i < values.size(); ++i)
1442  stream << '\t' << values[i];
1443  stream << '\n';
1444  }
1445  }
1446 
1447 
1448 
1449  //----------------------------------------------------------------------//
1450 
1451  GmvStream::GmvStream(std::ostream &out, const DataOutBase::GmvFlags &f)
1452  : StreamBase<DataOutBase::GmvFlags>(out, f)
1453  {}
1454 
1455 
1456  template <int dim>
1457  void
1458  GmvStream::write_point(const unsigned int, const Point<dim> &p)
1459  {
1460  Assert(selected_component != numbers::invalid_unsigned_int,
1461  ExcNotInitialized());
1462  stream << p(selected_component) << ' ';
1463  }
1464 
1465 
1466 
1467  template <int dim>
1468  void
1469  GmvStream::write_cell(const unsigned int,
1470  const unsigned int s,
1471  const std::array<unsigned int, dim> &offsets)
1472  {
1473  // Vertices are numbered starting with one.
1474  const unsigned int start = s + 1;
1475  stream << gmv_cell_type[dim] << '\n';
1476 
1477  switch (dim)
1478  {
1479  case 0:
1480  {
1481  stream << start;
1482  break;
1483  }
1484 
1485  case 1:
1486  {
1487  const unsigned int d1 = offsets[0];
1488  stream << start;
1489  stream << '\t' << start + d1;
1490  break;
1491  }
1492 
1493  case 2:
1494  {
1495  const unsigned int d1 = offsets[0];
1496  const unsigned int d2 = offsets[1];
1497  stream << start;
1498  stream << '\t' << start + d1;
1499  stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1500  break;
1501  }
1502 
1503  case 3:
1504  {
1505  const unsigned int d1 = offsets[0];
1506  const unsigned int d2 = offsets[1];
1507  const unsigned int d3 = offsets[2];
1508  stream << start;
1509  stream << '\t' << start + d1;
1510  stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1511  stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
1512  << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
1513  break;
1514  }
1515 
1516  default:
1517  Assert(false, ExcNotImplemented());
1518  }
1519  stream << '\n';
1520  }
1521 
1522 
1523 
1524  TecplotStream::TecplotStream(std::ostream &out,
1525  const DataOutBase::TecplotFlags &f)
1526  : StreamBase<DataOutBase::TecplotFlags>(out, f)
1527  {}
1528 
1529 
1530  template <int dim>
1531  void
1532  TecplotStream::write_point(const unsigned int, const Point<dim> &p)
1533  {
1534  Assert(selected_component != numbers::invalid_unsigned_int,
1535  ExcNotInitialized());
1536  stream << p(selected_component) << '\n';
1537  }
1538 
1539 
1540 
1541  template <int dim>
1542  void
1543  TecplotStream::write_cell(const unsigned int,
1544  const unsigned int s,
1545  const std::array<unsigned int, dim> &offsets)
1546  {
1547  const unsigned int start = s + 1;
1548 
1549  switch (dim)
1550  {
1551  case 0:
1552  {
1553  stream << start;
1554  break;
1555  }
1556 
1557  case 1:
1558  {
1559  const unsigned int d1 = offsets[0];
1560  stream << start;
1561  stream << '\t' << start + d1;
1562  break;
1563  }
1564 
1565  case 2:
1566  {
1567  const unsigned int d1 = offsets[0];
1568  const unsigned int d2 = offsets[1];
1569  stream << start;
1570  stream << '\t' << start + d1;
1571  stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1572  break;
1573  }
1574 
1575  case 3:
1576  {
1577  const unsigned int d1 = offsets[0];
1578  const unsigned int d2 = offsets[1];
1579  const unsigned int d3 = offsets[2];
1580  stream << start;
1581  stream << '\t' << start + d1;
1582  stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1583  stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
1584  << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
1585  break;
1586  }
1587 
1588  default:
1589  Assert(false, ExcNotImplemented());
1590  }
1591  stream << '\n';
1592  }
1593 
1594 
1595 
1596  UcdStream::UcdStream(std::ostream &out, const DataOutBase::UcdFlags &f)
1597  : StreamBase<DataOutBase::UcdFlags>(out, f)
1598  {}
1599 
1600 
1601  template <int dim>
1602  void
1603  UcdStream::write_point(const unsigned int index, const Point<dim> &p)
1604  {
1605  stream << index + 1 << " ";
1606  // write out coordinates
1607  for (unsigned int i = 0; i < dim; ++i)
1608  stream << p(i) << ' ';
1609  // fill with zeroes
1610  for (unsigned int i = dim; i < 3; ++i)
1611  stream << "0 ";
1612  stream << '\n';
1613  }
1614 
1615 
1616 
1617  template <int dim>
1618  void
1619  UcdStream::write_cell(const unsigned int index,
1620  const unsigned int start,
1621  const std::array<unsigned int, dim> &offsets)
1622  {
1623  const auto nodes =
1624  DataOutBaseImplementation::set_node_numbers(start, offsets);
1625 
1626  // Write out all cells and remember that all indices must be shifted by one.
1627  stream << index + 1 << "\t0 " << ucd_cell_type[dim];
1628  for (unsigned int i = 0; i < nodes.size(); ++i)
1629  stream << '\t' << nodes[GeometryInfo<dim>::ucd_to_deal[i]] + 1;
1630  stream << '\n';
1631  }
1632 
1633 
1634 
1635  template <typename data>
1636  inline void
1637  UcdStream::write_dataset(const unsigned int index,
1638  const std::vector<data> &values)
1639  {
1640  stream << index + 1;
1641  for (unsigned int i = 0; i < values.size(); ++i)
1642  stream << '\t' << values[i];
1643  stream << '\n';
1644  }
1645 
1646 
1647 
1648  //----------------------------------------------------------------------//
1649 
1650  VtkStream::VtkStream(std::ostream &out, const DataOutBase::VtkFlags &f)
1651  : StreamBase<DataOutBase::VtkFlags>(out, f)
1652  {}
1653 
1654 
1655  template <int dim>
1656  void
1657  VtkStream::write_point(const unsigned int, const Point<dim> &p)
1658  {
1659  // write out coordinates
1660  stream << p;
1661  // fill with zeroes
1662  for (unsigned int i = dim; i < 3; ++i)
1663  stream << " 0";
1664  stream << '\n';
1665  }
1666 
1667 
1668 
1669  template <int dim>
1670  void
1671  VtkStream::write_cell(const unsigned int,
1672  const unsigned int start,
1673  const std::array<unsigned int, dim> &offsets)
1674  {
1675  stream << GeometryInfo<dim>::vertices_per_cell << '\t';
1676 
1677  switch (dim)
1678  {
1679  case 0:
1680  {
1681  stream << start;
1682  break;
1683  }
1684 
1685  case 1:
1686  {
1687  const unsigned int d1 = offsets[0];
1688  stream << start;
1689  stream << '\t' << start + d1;
1690  break;
1691  }
1692 
1693  case 2:
1694  {
1695  const unsigned int d1 = offsets[0];
1696  const unsigned int d2 = offsets[1];
1697  stream << start;
1698  stream << '\t' << start + d1;
1699  stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1700  break;
1701  }
1702 
1703  case 3:
1704  {
1705  const unsigned int d1 = offsets[0];
1706  const unsigned int d2 = offsets[1];
1707  const unsigned int d3 = offsets[2];
1708  stream << start;
1709  stream << '\t' << start + d1;
1710  stream << '\t' << start + d2 + d1 << '\t' << start + d2;
1711  stream << '\t' << start + d3 << '\t' << start + d3 + d1 << '\t'
1712  << start + d3 + d2 + d1 << '\t' << start + d3 + d2;
1713  break;
1714  }
1715 
1716  default:
1717  Assert(false, ExcNotImplemented());
1718  }
1719  stream << '\n';
1720  }
1721 
1722 
1723 
1724  void
1725  VtkStream::write_cell_single(const unsigned int index,
1726  const unsigned int start,
1727  const unsigned int n_points,
1729  {
1730  (void)index;
1731 
1732  static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
1733 
1734  stream << '\t' << n_points;
1735  for (unsigned int i = 0; i < n_points; ++i)
1736  stream << '\t'
1737  << start +
1738  (reference_cell == ReferenceCells::Pyramid ? table[i] : i);
1739  stream << '\n';
1740  }
1741 
1742  template <int dim>
1743  void
1744  VtkStream::write_high_order_cell(const unsigned int start,
1745  const std::vector<unsigned> &connectivity)
1746  {
1747  stream << connectivity.size();
1748  for (const auto &c : connectivity)
1749  stream << '\t' << start + c;
1750  stream << '\n';
1751  }
1752 } // namespace
1753 
1754 
1755 
1756 namespace DataOutBase
1757 {
1758  const unsigned int Deal_II_IntermediateFlags::format_version = 4;
1759 
1760 
1761  template <int dim, int spacedim>
1762  const unsigned int Patch<dim, spacedim>::space_dim;
1763 
1764 
1765  template <int dim, int spacedim>
1766  const unsigned int Patch<dim, spacedim>::no_neighbor;
1767 
1768 
1769  template <int dim, int spacedim>
1771  : patch_index(no_neighbor)
1772  , n_subdivisions(1)
1773  , points_are_available(false)
1775  // all the other data has a constructor of its own, except for the "neighbors"
1776  // field, which we set to invalid values.
1777  {
1778  for (const unsigned int i : GeometryInfo<dim>::face_indices())
1779  neighbors[i] = no_neighbor;
1780 
1781  AssertIndexRange(dim, spacedim + 1);
1782  Assert(spacedim <= 3, ExcNotImplemented());
1783  }
1784 
1785 
1786 
1787  template <int dim, int spacedim>
1788  bool
1790  {
1791  if (reference_cell != patch.reference_cell)
1792  return false;
1793 
1794  // TODO: make tolerance relative
1795  const double epsilon = 3e-16;
1796  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1797  if (vertices[i].distance(patch.vertices[i]) > epsilon)
1798  return false;
1799 
1800  for (const unsigned int i : GeometryInfo<dim>::face_indices())
1801  if (neighbors[i] != patch.neighbors[i])
1802  return false;
1803 
1804  if (patch_index != patch.patch_index)
1805  return false;
1806 
1807  if (n_subdivisions != patch.n_subdivisions)
1808  return false;
1809 
1810  if (points_are_available != patch.points_are_available)
1811  return false;
1812 
1813  if (data.n_rows() != patch.data.n_rows())
1814  return false;
1815 
1816  if (data.n_cols() != patch.data.n_cols())
1817  return false;
1818 
1819  for (unsigned int i = 0; i < data.n_rows(); ++i)
1820  for (unsigned int j = 0; j < data.n_cols(); ++j)
1821  if (data[i][j] != patch.data[i][j])
1822  return false;
1823 
1824  return true;
1825  }
1826 
1827 
1828 
1829  template <int dim, int spacedim>
1830  std::size_t
1832  {
1833  return (sizeof(vertices) / sizeof(vertices[0]) *
1835  sizeof(neighbors) / sizeof(neighbors[0]) *
1838  MemoryConsumption::memory_consumption(n_subdivisions) +
1840  MemoryConsumption::memory_consumption(points_are_available) +
1841  sizeof(reference_cell));
1842  }
1843 
1844 
1845 
1846  template <int dim, int spacedim>
1847  void
1849  {
1850  std::swap(vertices, other_patch.vertices);
1851  std::swap(neighbors, other_patch.neighbors);
1852  std::swap(patch_index, other_patch.patch_index);
1853  std::swap(n_subdivisions, other_patch.n_subdivisions);
1854  data.swap(other_patch.data);
1855  std::swap(points_are_available, other_patch.points_are_available);
1856  std::swap(reference_cell, other_patch.reference_cell);
1857  }
1858 
1859 
1860 
1861  template <int spacedim>
1862  const unsigned int Patch<0, spacedim>::space_dim;
1863 
1864 
1865  template <int spacedim>
1866  const unsigned int Patch<0, spacedim>::no_neighbor;
1867 
1868 
1869  template <int spacedim>
1870  unsigned int Patch<0, spacedim>::neighbors[1] = {
1872 
1873  template <int spacedim>
1874  const unsigned int Patch<0, spacedim>::n_subdivisions = 1;
1875 
1876  template <int spacedim>
1879 
1880  template <int spacedim>
1882  : patch_index(no_neighbor)
1883  , points_are_available(false)
1884  {
1885  Assert(spacedim <= 3, ExcNotImplemented());
1886  }
1887 
1888 
1889 
1890  template <int spacedim>
1891  bool
1893  {
1894  const unsigned int dim = 0;
1895 
1896  // TODO: make tolerance relative
1897  const double epsilon = 3e-16;
1898  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
1899  if (vertices[i].distance(patch.vertices[i]) > epsilon)
1900  return false;
1901 
1902  if (patch_index != patch.patch_index)
1903  return false;
1904 
1906  return false;
1907 
1908  if (data.n_rows() != patch.data.n_rows())
1909  return false;
1910 
1911  if (data.n_cols() != patch.data.n_cols())
1912  return false;
1913 
1914  for (unsigned int i = 0; i < data.n_rows(); ++i)
1915  for (unsigned int j = 0; j < data.n_cols(); ++j)
1916  if (data[i][j] != patch.data[i][j])
1917  return false;
1918 
1919  return true;
1920  }
1921 
1922 
1923 
1924  template <int spacedim>
1925  std::size_t
1927  {
1928  return (sizeof(vertices) / sizeof(vertices[0]) *
1932  }
1933 
1934 
1935 
1936  template <int spacedim>
1937  void
1939  {
1940  std::swap(vertices, other_patch.vertices);
1941  std::swap(patch_index, other_patch.patch_index);
1942  data.swap(other_patch.data);
1944  }
1945 
1946 
1947 
1948  UcdFlags::UcdFlags(const bool write_preamble)
1949  : write_preamble(write_preamble)
1950  {}
1951 
1952 
1953 
1955  {
1956  space_dimension_labels.emplace_back("x");
1957  space_dimension_labels.emplace_back("y");
1958  space_dimension_labels.emplace_back("z");
1959  }
1960 
1961 
1962 
1963  GnuplotFlags::GnuplotFlags(const std::vector<std::string> &labels)
1964  : space_dimension_labels(labels)
1965  {}
1966 
1967 
1968 
1969  std::size_t
1971  {
1973  }
1974 
1975 
1976 
1977  PovrayFlags::PovrayFlags(const bool smooth,
1978  const bool bicubic_patch,
1979  const bool external_data)
1980  : smooth(smooth)
1981  , bicubic_patch(bicubic_patch)
1982  , external_data(external_data)
1983  {}
1984 
1985 
1986  DataOutFilterFlags::DataOutFilterFlags(const bool filter_duplicate_vertices,
1987  const bool xdmf_hdf5_output)
1988  : filter_duplicate_vertices(filter_duplicate_vertices)
1989  , xdmf_hdf5_output(xdmf_hdf5_output)
1990  {}
1991 
1992 
1993  void
1995  {
1996  prm.declare_entry(
1997  "Filter duplicate vertices",
1998  "false",
1999  Patterns::Bool(),
2000  "Whether to remove duplicate vertex values. deal.II duplicates "
2001  "vertices once for each adjacent cell so that it can output "
2002  "discontinuous quantities for which there may be more than one "
2003  "value for each vertex position. Setting this flag to "
2004  "'true' will merge all of these values by selecting a "
2005  "random one and outputting this as 'the' value for the vertex. "
2006  "As long as the data to be output corresponds to continuous "
2007  "fields, merging vertices has no effect. On the other hand, "
2008  "if the data to be output corresponds to discontinuous fields "
2009  "(either because you are using a discontinuous finite element, "
2010  "or because you are using a DataPostprocessor that yields "
2011  "discontinuous data, or because the data to be output has been "
2012  "produced by entirely different means), then the data in the "
2013  "output file no longer faithfully represents the underlying data "
2014  "because the discontinuous field has been replaced by a "
2015  "continuous one. Note also that the filtering can not occur "
2016  "on processor boundaries. Thus, a filtered discontinuous field "
2017  "looks like a continuous field inside of a subdomain, "
2018  "but like a discontinuous field at the subdomain boundary."
2019  "\n\n"
2020  "In any case, filtering results in drastically smaller output "
2021  "files (smaller by about a factor of 2^dim).");
2022  prm.declare_entry(
2023  "XDMF HDF5 output",
2024  "false",
2025  Patterns::Bool(),
2026  "Whether the data will be used in an XDMF/HDF5 combination.");
2027  }
2028 
2029 
2030 
2031  void
2033  {
2034  filter_duplicate_vertices = prm.get_bool("Filter duplicate vertices");
2035  xdmf_hdf5_output = prm.get_bool("XDMF HDF5 output");
2036  }
2037 
2038 
2039 
2040  DXFlags::DXFlags(const bool write_neighbors,
2041  const bool int_binary,
2042  const bool coordinates_binary,
2043  const bool data_binary)
2044  : write_neighbors(write_neighbors)
2045  , int_binary(int_binary)
2046  , coordinates_binary(coordinates_binary)
2047  , data_binary(data_binary)
2048  , data_double(false)
2049  {}
2050 
2051 
2052  void
2054  {
2055  prm.declare_entry("Write neighbors",
2056  "true",
2057  Patterns::Bool(),
2058  "A boolean field indicating whether neighborship "
2059  "information between cells is to be written to the "
2060  "OpenDX output file");
2061  prm.declare_entry("Integer format",
2062  "ascii",
2063  Patterns::Selection("ascii|32|64"),
2064  "Output format of integer numbers, which is "
2065  "either a text representation (ascii) or binary integer "
2066  "values of 32 or 64 bits length");
2067  prm.declare_entry("Coordinates format",
2068  "ascii",
2069  Patterns::Selection("ascii|32|64"),
2070  "Output format of vertex coordinates, which is "
2071  "either a text representation (ascii) or binary "
2072  "floating point values of 32 or 64 bits length");
2073  prm.declare_entry("Data format",
2074  "ascii",
2075  Patterns::Selection("ascii|32|64"),
2076  "Output format of data values, which is "
2077  "either a text representation (ascii) or binary "
2078  "floating point values of 32 or 64 bits length");
2079  }
2080 
2081 
2082 
2083  void
2085  {
2086  write_neighbors = prm.get_bool("Write neighbors");
2087  // TODO:[GK] Read the new parameters
2088  }
2089 
2090 
2091 
2092  void
2094  {
2095  prm.declare_entry("Write preamble",
2096  "true",
2097  Patterns::Bool(),
2098  "A flag indicating whether a comment should be "
2099  "written to the beginning of the output file "
2100  "indicating date and time of creation as well "
2101  "as the creating program");
2102  }
2103 
2104 
2105 
2106  void
2108  {
2109  write_preamble = prm.get_bool("Write preamble");
2110  }
2111 
2112 
2113 
2114  SvgFlags::SvgFlags(const unsigned int height_vector,
2115  const int azimuth_angle,
2116  const int polar_angle,
2117  const unsigned int line_thickness,
2118  const bool margin,
2119  const bool draw_colorbar)
2120  : height(4000)
2121  , width(0)
2122  , height_vector(height_vector)
2123  , azimuth_angle(azimuth_angle)
2124  , polar_angle(polar_angle)
2125  , line_thickness(line_thickness)
2126  , margin(margin)
2127  , draw_colorbar(draw_colorbar)
2128  {}
2129 
2130 
2131 
2132  void
2134  {
2135  prm.declare_entry("Use smooth triangles",
2136  "false",
2137  Patterns::Bool(),
2138  "A flag indicating whether POVRAY should use smoothed "
2139  "triangles instead of the usual ones");
2140  prm.declare_entry("Use bicubic patches",
2141  "false",
2142  Patterns::Bool(),
2143  "Whether POVRAY should use bicubic patches");
2144  prm.declare_entry("Include external file",
2145  "true",
2146  Patterns::Bool(),
2147  "Whether camera and lighting information should "
2148  "be put into an external file \"data.inc\" or into "
2149  "the POVRAY input file");
2150  }
2151 
2152 
2153 
2154  void
2156  {
2157  smooth = prm.get_bool("Use smooth triangles");
2158  bicubic_patch = prm.get_bool("Use bicubic patches");
2159  external_data = prm.get_bool("Include external file");
2160  }
2161 
2162 
2163 
2164  EpsFlags::EpsFlags(const unsigned int height_vector,
2165  const unsigned int color_vector,
2166  const SizeType size_type,
2167  const unsigned int size,
2168  const double line_width,
2169  const double azimut_angle,
2170  const double turn_angle,
2171  const double z_scaling,
2172  const bool draw_mesh,
2173  const bool draw_cells,
2174  const bool shade_cells,
2175  const ColorFunction color_function)
2176  : height_vector(height_vector)
2177  , color_vector(color_vector)
2178  , size_type(size_type)
2179  , size(size)
2180  , line_width(line_width)
2181  , azimut_angle(azimut_angle)
2182  , turn_angle(turn_angle)
2183  , z_scaling(z_scaling)
2184  , draw_mesh(draw_mesh)
2185  , draw_cells(draw_cells)
2186  , shade_cells(shade_cells)
2187  , color_function(color_function)
2188  {}
2189 
2190 
2191 
2194  const double xmin,
2195  const double xmax)
2196  {
2197  RgbValues rgb_values = {0, 0, 0};
2198 
2199  // A difficult color scale:
2200  // xmin = black (1)
2201  // 3/4*xmin+1/4*xmax = blue (2)
2202  // 1/2*xmin+1/2*xmax = green (3)
2203  // 1/4*xmin+3/4*xmax = red (4)
2204  // xmax = white (5)
2205  // Makes the following color functions:
2206  //
2207  // red green blue
2208  // __
2209  // / /\ / /\ /
2210  // ____/ __/ \/ / \__/
2211 
2212  // { 0 (1) - (3)
2213  // r = { ( 4*x-2*xmin+2*xmax)/(xmax-xmin) (3) - (4)
2214  // { 1 (4) - (5)
2215  //
2216  // { 0 (1) - (2)
2217  // g = { ( 4*x-3*xmin- xmax)/(xmax-xmin) (2) - (3)
2218  // { (-4*x+ xmin+3*xmax)/(xmax-xmin) (3) - (4)
2219  // { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5)
2220  //
2221  // { ( 4*x-4*xmin )/(xmax-xmin) (1) - (2)
2222  // b = { (-4*x+2*xmin+2*xmax)/(xmax-xmin) (2) - (3)
2223  // { 0 (3) - (4)
2224  // { ( 4*x- xmin-3*xmax)/(xmax-xmin) (4) - (5)
2225 
2226  double sum = xmax + xmin;
2227  double sum13 = xmin + 3 * xmax;
2228  double sum22 = 2 * xmin + 2 * xmax;
2229  double sum31 = 3 * xmin + xmax;
2230  double dif = xmax - xmin;
2231  double rezdif = 1.0 / dif;
2232 
2233  int where;
2234 
2235  if (x < (sum31) / 4)
2236  where = 0;
2237  else if (x < (sum22) / 4)
2238  where = 1;
2239  else if (x < (sum13) / 4)
2240  where = 2;
2241  else
2242  where = 3;
2243 
2244  if (dif != 0)
2245  {
2246  switch (where)
2247  {
2248  case 0:
2249  rgb_values.red = 0;
2250  rgb_values.green = 0;
2251  rgb_values.blue = (x - xmin) * 4. * rezdif;
2252  break;
2253  case 1:
2254  rgb_values.red = 0;
2255  rgb_values.green = (4 * x - 3 * xmin - xmax) * rezdif;
2256  rgb_values.blue = (sum22 - 4. * x) * rezdif;
2257  break;
2258  case 2:
2259  rgb_values.red = (4 * x - 2 * sum) * rezdif;
2260  rgb_values.green = (xmin + 3 * xmax - 4 * x) * rezdif;
2261  rgb_values.blue = 0;
2262  break;
2263  case 3:
2264  rgb_values.red = 1;
2265  rgb_values.green = (4 * x - xmin - 3 * xmax) * rezdif;
2266  rgb_values.blue = (4. * x - sum13) * rezdif;
2267  break;
2268  default:
2269  break;
2270  }
2271  }
2272  else // White
2273  rgb_values.red = rgb_values.green = rgb_values.blue = 1;
2274 
2275  return rgb_values;
2276  }
2277 
2278 
2279 
2282  const double xmin,
2283  const double xmax)
2284  {
2285  EpsFlags::RgbValues rgb_values;
2286  rgb_values.red = rgb_values.blue = rgb_values.green =
2287  (x - xmin) / (xmax - xmin);
2288  return rgb_values;
2289  }
2290 
2291 
2292 
2295  const double xmin,
2296  const double xmax)
2297  {
2298  EpsFlags::RgbValues rgb_values;
2299  rgb_values.red = rgb_values.blue = rgb_values.green =
2300  1 - (x - xmin) / (xmax - xmin);
2301  return rgb_values;
2302  }
2303 
2304 
2305 
2306  void
2308  {
2309  prm.declare_entry("Index of vector for height",
2310  "0",
2312  "Number of the input vector that is to be used to "
2313  "generate height information");
2314  prm.declare_entry("Index of vector for color",
2315  "0",
2317  "Number of the input vector that is to be used to "
2318  "generate color information");
2319  prm.declare_entry("Scale to width or height",
2320  "width",
2321  Patterns::Selection("width|height"),
2322  "Whether width or height should be scaled to match "
2323  "the given size");
2324  prm.declare_entry("Size (width or height) in eps units",
2325  "300",
2327  "The size (width or height) to which the eps output "
2328  "file is to be scaled");
2329  prm.declare_entry("Line widths in eps units",
2330  "0.5",
2331  Patterns::Double(),
2332  "The width in which the postscript renderer is to "
2333  "plot lines");
2334  prm.declare_entry("Azimut angle",
2335  "60",
2336  Patterns::Double(0, 180),
2337  "Angle of the viewing position against the vertical "
2338  "axis");
2339  prm.declare_entry("Turn angle",
2340  "30",
2341  Patterns::Double(0, 360),
2342  "Angle of the viewing direction against the y-axis");
2343  prm.declare_entry("Scaling for z-axis",
2344  "1",
2345  Patterns::Double(),
2346  "Scaling for the z-direction relative to the scaling "
2347  "used in x- and y-directions");
2348  prm.declare_entry("Draw mesh lines",
2349  "true",
2350  Patterns::Bool(),
2351  "Whether the mesh lines, or only the surface should be "
2352  "drawn");
2353  prm.declare_entry("Fill interior of cells",
2354  "true",
2355  Patterns::Bool(),
2356  "Whether only the mesh lines, or also the interior of "
2357  "cells should be plotted. If this flag is false, then "
2358  "one can see through the mesh");
2359  prm.declare_entry("Color shading of interior of cells",
2360  "true",
2361  Patterns::Bool(),
2362  "Whether the interior of cells shall be shaded");
2363  prm.declare_entry("Color function",
2364  "default",
2366  "default|grey scale|reverse grey scale"),
2367  "Name of a color function used to colorize mesh lines "
2368  "and/or cell interiors");
2369  }
2370 
2371 
2372 
2373  void
2375  {
2376  height_vector = prm.get_integer("Index of vector for height");
2377  color_vector = prm.get_integer("Index of vector for color");
2378  if (prm.get("Scale to width or height") == "width")
2379  size_type = width;
2380  else
2381  size_type = height;
2382  size = prm.get_integer("Size (width or height) in eps units");
2383  line_width = prm.get_double("Line widths in eps units");
2384  azimut_angle = prm.get_double("Azimut angle");
2385  turn_angle = prm.get_double("Turn angle");
2386  z_scaling = prm.get_double("Scaling for z-axis");
2387  draw_mesh = prm.get_bool("Draw mesh lines");
2388  draw_cells = prm.get_bool("Fill interior of cells");
2389  shade_cells = prm.get_bool("Color shading of interior of cells");
2390  if (prm.get("Color function") == "default")
2392  else if (prm.get("Color function") == "grey scale")
2394  else if (prm.get("Color function") == "reverse grey scale")
2396  else
2397  // we shouldn't get here, since the parameter object should already have
2398  // checked that the given value is valid
2399  Assert(false, ExcInternalError());
2400  }
2401 
2402 
2403  Hdf5Flags::Hdf5Flags(const CompressionLevel compression_level)
2404  : compression_level(compression_level)
2405  {}
2406 
2407 
2408  TecplotFlags::TecplotFlags(const char *zone_name, const double solution_time)
2409  : zone_name(zone_name)
2410  , solution_time(solution_time)
2411  {}
2412 
2413 
2414 
2415  std::size_t
2417  {
2418  return sizeof(*this) + MemoryConsumption::memory_consumption(zone_name);
2419  }
2420 
2421 
2422 
2423  VtkFlags::VtkFlags(const double time,
2424  const unsigned int cycle,
2425  const bool print_date_and_time,
2426  const CompressionLevel compression_level,
2427  const bool write_higher_order_cells,
2428  const std::map<std::string, std::string> &physical_units)
2429  : time(time)
2430  , cycle(cycle)
2431  , print_date_and_time(print_date_and_time)
2432  , compression_level(compression_level)
2433  , write_higher_order_cells(write_higher_order_cells)
2434  , physical_units(physical_units)
2435  {}
2436 
2437 
2438 
2439  OutputFormat
2440  parse_output_format(const std::string &format_name)
2441  {
2442  if (format_name == "none")
2443  return none;
2444 
2445  if (format_name == "dx")
2446  return dx;
2447 
2448  if (format_name == "ucd")
2449  return ucd;
2450 
2451  if (format_name == "gnuplot")
2452  return gnuplot;
2453 
2454  if (format_name == "povray")
2455  return povray;
2456 
2457  if (format_name == "eps")
2458  return eps;
2459 
2460  if (format_name == "gmv")
2461  return gmv;
2462 
2463  if (format_name == "tecplot")
2464  return tecplot;
2465 
2466  if (format_name == "vtk")
2467  return vtk;
2468 
2469  if (format_name == "vtu")
2470  return vtu;
2471 
2472  if (format_name == "deal.II intermediate")
2473  return deal_II_intermediate;
2474 
2475  if (format_name == "hdf5")
2476  return hdf5;
2477 
2478  AssertThrow(false,
2479  ExcMessage("The given file format name is not recognized: <" +
2480  format_name + ">"));
2481 
2482  // return something invalid
2483  return OutputFormat(-1);
2484  }
2485 
2486 
2487 
2488  std::string
2490  {
2491  return "none|dx|ucd|gnuplot|povray|eps|gmv|tecplot|vtk|vtu|hdf5|svg|deal.II intermediate";
2492  }
2493 
2494 
2495 
2496  std::string
2497  default_suffix(const OutputFormat output_format)
2498  {
2499  switch (output_format)
2500  {
2501  case none:
2502  return "";
2503  case dx:
2504  return ".dx";
2505  case ucd:
2506  return ".inp";
2507  case gnuplot:
2508  return ".gnuplot";
2509  case povray:
2510  return ".pov";
2511  case eps:
2512  return ".eps";
2513  case gmv:
2514  return ".gmv";
2515  case tecplot:
2516  return ".dat";
2517  case vtk:
2518  return ".vtk";
2519  case vtu:
2520  return ".vtu";
2521  case deal_II_intermediate:
2522  return ".d2";
2523  case hdf5:
2524  return ".h5";
2525  case svg:
2526  return ".svg";
2527  default:
2528  Assert(false, ExcNotImplemented());
2529  return "";
2530  }
2531  }
2532 
2533 
2534  //----------------------------------------------------------------------//
2535 
2536 
2541  template <int dim, int spacedim>
2542  std::vector<Point<spacedim>>
2543  get_node_positions(const std::vector<Patch<dim, spacedim>> &patches)
2544  {
2545  Assert(dim <= 3, ExcNotImplemented());
2546  static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
2547 
2548  std::vector<Point<spacedim>> node_positions;
2549  for (const auto &patch : patches)
2550  {
2551  // special treatment of non-hypercube cells
2552  if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2553  {
2554  for (unsigned int point_no = 0; point_no < patch.data.n_cols();
2555  ++point_no)
2556  node_positions.emplace_back(get_node_location(
2557  patch,
2559  table[point_no] :
2560  point_no)));
2561  }
2562  else
2563  {
2564  const unsigned int n_subdivisions = patch.n_subdivisions;
2565  const unsigned int n = n_subdivisions + 1;
2566 
2567  switch (dim)
2568  {
2569  case 0:
2570  node_positions.emplace_back(
2571  get_equispaced_location(patch, {}, n_subdivisions));
2572  break;
2573  case 1:
2574  for (unsigned int i1 = 0; i1 < n; ++i1)
2575  node_positions.emplace_back(
2576  get_equispaced_location(patch, {i1}, n_subdivisions));
2577  break;
2578  case 2:
2579  for (unsigned int i2 = 0; i2 < n; ++i2)
2580  for (unsigned int i1 = 0; i1 < n; ++i1)
2581  node_positions.emplace_back(get_equispaced_location(
2582  patch, {i1, i2}, n_subdivisions));
2583  break;
2584  case 3:
2585  for (unsigned int i3 = 0; i3 < n; ++i3)
2586  for (unsigned int i2 = 0; i2 < n; ++i2)
2587  for (unsigned int i1 = 0; i1 < n; ++i1)
2588  node_positions.emplace_back(get_equispaced_location(
2589  patch, {i1, i2, i3}, n_subdivisions));
2590  break;
2591 
2592  default:
2593  Assert(false, ExcInternalError());
2594  }
2595  }
2596  }
2597 
2598  return node_positions;
2599  }
2600 
2601 
2602  template <int dim, int spacedim, typename StreamType>
2603  void
2604  write_nodes(const std::vector<Patch<dim, spacedim>> &patches, StreamType &out)
2605  {
2606  // Obtain the node locations, and then output them via the given stream
2607  // object
2608  const std::vector<Point<spacedim>> node_positions =
2609  get_node_positions(patches);
2610 
2611  int count = 0;
2612  for (const auto &node : node_positions)
2613  out.write_point(count++, node);
2614  out.flush_points();
2615  }
2616 
2617 
2618 
2619  template <int dim, int spacedim, typename StreamType>
2620  void
2621  write_cells(const std::vector<Patch<dim, spacedim>> &patches, StreamType &out)
2622  {
2623  Assert(dim <= 3, ExcNotImplemented());
2624  unsigned int count = 0;
2625  unsigned int first_vertex_of_patch = 0;
2626  for (const auto &patch : patches)
2627  {
2628  // special treatment of simplices since they are not subdivided
2629  if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2630  {
2631  out.write_cell_single(count++,
2632  first_vertex_of_patch,
2633  patch.data.n_cols(),
2634  patch.reference_cell);
2635  first_vertex_of_patch += patch.data.n_cols();
2636  }
2637  else // hypercube cell
2638  {
2639  const unsigned int n_subdivisions = patch.n_subdivisions;
2640  const unsigned int n = n_subdivisions + 1;
2641 
2642  switch (dim)
2643  {
2644  case 0:
2645  {
2646  const unsigned int offset = first_vertex_of_patch;
2647  out.template write_cell<0>(count++, offset, {});
2648  break;
2649  }
2650 
2651  case 1:
2652  {
2653  constexpr unsigned int d1 = 1;
2654 
2655  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2656  {
2657  const unsigned int offset =
2658  first_vertex_of_patch + i1 * d1;
2659  out.template write_cell<1>(count++, offset, {{d1}});
2660  }
2661 
2662  break;
2663  }
2664 
2665  case 2:
2666  {
2667  constexpr unsigned int d1 = 1;
2668  const unsigned int d2 = n;
2670  for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2671  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2672  {
2673  const unsigned int offset =
2674  first_vertex_of_patch + i2 * d2 + i1 * d1;
2675  out.template write_cell<2>(count++,
2676  offset,
2677  {{d1, d2}});
2678  }
2679 
2680  break;
2681  }
2682 
2683  case 3:
2684  {
2685  constexpr unsigned int d1 = 1;
2686  const unsigned int d2 = n;
2687  const unsigned int d3 = n * n;
2688 
2689  for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
2690  for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
2691  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
2692  {
2693  const unsigned int offset = first_vertex_of_patch +
2694  i3 * d3 + i2 * d2 +
2695  i1 * d1;
2696  out.template write_cell<3>(count++,
2697  offset,
2698  {{d1, d2, d3}});
2699  }
2700 
2701  break;
2702  }
2703  default:
2705  }
2706 
2707  // Update the number of the first vertex of this patch
2708  first_vertex_of_patch +=
2709  Utilities::fixed_power<dim>(n_subdivisions + 1);
2710  }
2711  }
2712 
2713  out.flush_cells();
2714  }
2715 
2716 
2717 
2718  template <int dim, int spacedim, typename StreamType>
2719  void
2720  write_high_order_cells(const std::vector<Patch<dim, spacedim>> &patches,
2721  StreamType &out,
2722  const bool legacy_format)
2723  {
2724  Assert(dim <= 3 && dim > 1, ExcNotImplemented());
2725  unsigned int first_vertex_of_patch = 0;
2726  // Array to hold all the node numbers of a cell
2727  std::vector<unsigned> connectivity;
2728 
2729  for (const auto &patch : patches)
2730  {
2731  if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
2732  {
2733  connectivity.resize(patch.data.n_cols());
2734 
2735  for (unsigned int i = 0; i < patch.data.n_cols(); ++i)
2736  connectivity[i] = i;
2737 
2738  out.template write_high_order_cell<dim>(first_vertex_of_patch,
2739  connectivity);
2740 
2741  first_vertex_of_patch += patch.data.n_cols();
2742  }
2743  else
2744  {
2745  const unsigned int n_subdivisions = patch.n_subdivisions;
2746  const unsigned int n = n_subdivisions + 1;
2747 
2748  connectivity.resize(Utilities::fixed_power<dim>(n));
2749 
2750  switch (dim)
2751  {
2752  case 0:
2753  {
2754  Assert(false,
2755  ExcMessage("Point-like cells should not be possible "
2756  "when writing higher-order cells."));
2757  break;
2758  }
2759  case 1:
2760  {
2761  for (unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2762  {
2763  const unsigned int local_index = i1;
2764  const unsigned int connectivity_index =
2765  patch.reference_cell
2766  .template vtk_lexicographic_to_node_index<1>(
2767  {{i1}}, {{n_subdivisions}}, legacy_format);
2768  connectivity[connectivity_index] = local_index;
2769  }
2770 
2771  break;
2772  }
2773  case 2:
2774  {
2775  for (unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2776  for (unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2777  {
2778  const unsigned int local_index = i2 * n + i1;
2779  const unsigned int connectivity_index =
2780  patch.reference_cell
2781  .template vtk_lexicographic_to_node_index<2>(
2782  {{i1, i2}},
2783  {{n_subdivisions, n_subdivisions}},
2784  legacy_format);
2785  connectivity[connectivity_index] = local_index;
2786  }
2787 
2788  break;
2789  }
2790  case 3:
2791  {
2792  for (unsigned int i3 = 0; i3 < n_subdivisions + 1; ++i3)
2793  for (unsigned int i2 = 0; i2 < n_subdivisions + 1; ++i2)
2794  for (unsigned int i1 = 0; i1 < n_subdivisions + 1; ++i1)
2795  {
2796  const unsigned int local_index =
2797  i3 * n * n + i2 * n + i1;
2798  const unsigned int connectivity_index =
2799  patch.reference_cell
2800  .template vtk_lexicographic_to_node_index<3>(
2801  {{i1, i2, i3}},
2802  {{n_subdivisions,
2803  n_subdivisions,
2804  n_subdivisions}},
2805  legacy_format);
2806  connectivity[connectivity_index] = local_index;
2807  }
2808 
2809  break;
2810  }
2811  default:
2812  Assert(false, ExcNotImplemented());
2813  }
2814 
2815  // Having so set up the 'connectivity' data structure,
2816  // output it:
2817  out.template write_high_order_cell<dim>(first_vertex_of_patch,
2818  connectivity);
2819 
2820  // Finally update the number of the first vertex of this patch
2821  first_vertex_of_patch += Utilities::fixed_power<dim>(n);
2822  }
2823  }
2824 
2825  out.flush_cells();
2826  }
2827 
2828 
2829  template <int dim, int spacedim, typename StreamType>
2830  void
2831  write_data(const std::vector<Patch<dim, spacedim>> &patches,
2832  unsigned int n_data_sets,
2833  const bool double_precision,
2834  StreamType &out)
2835  {
2836  Assert(dim <= 3, ExcNotImplemented());
2837  unsigned int count = 0;
2838 
2839  for (const auto &patch : patches)
2840  {
2841  const unsigned int n_subdivisions = patch.n_subdivisions;
2842  const unsigned int n = n_subdivisions + 1;
2843  // Length of loops in all dimensions
2844  Assert((patch.data.n_rows() == n_data_sets &&
2845  !patch.points_are_available) ||
2846  (patch.data.n_rows() == n_data_sets + spacedim &&
2847  patch.points_are_available),
2849  (n_data_sets + spacedim) :
2850  n_data_sets,
2851  patch.data.n_rows()));
2852  Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
2853  ExcInvalidDatasetSize(patch.data.n_cols(), n));
2854 
2855  std::vector<float> floats(n_data_sets);
2856  std::vector<double> doubles(n_data_sets);
2858  // Data is already in lexicographic ordering
2859  for (unsigned int i = 0; i < Utilities::fixed_power<dim>(n);
2860  ++i, ++count)
2861  if (double_precision)
2862  {
2863  for (unsigned int data_set = 0; data_set < n_data_sets;
2864  ++data_set)
2865  doubles[data_set] = patch.data(data_set, i);
2866  out.write_dataset(count, doubles);
2867  }
2868  else
2869  {
2870  for (unsigned int data_set = 0; data_set < n_data_sets;
2871  ++data_set)
2872  floats[data_set] = patch.data(data_set, i);
2873  out.write_dataset(count, floats);
2874  }
2875  }
2876  }
2877 
2878 
2880  namespace
2881  {
2891  svg_project_point(Point<3> point,
2892  Point<3> camera_position,
2893  Point<3> camera_direction,
2894  Point<3> camera_horizontal,
2895  float camera_focus)
2896  {
2897  Point<3> camera_vertical;
2898  camera_vertical[0] = camera_horizontal[1] * camera_direction[2] -
2899  camera_horizontal[2] * camera_direction[1];
2900  camera_vertical[1] = camera_horizontal[2] * camera_direction[0] -
2901  camera_horizontal[0] * camera_direction[2];
2902  camera_vertical[2] = camera_horizontal[0] * camera_direction[1] -
2903  camera_horizontal[1] * camera_direction[0];
2904 
2905  float phi;
2906  phi = camera_focus;
2907  phi /= (point[0] - camera_position[0]) * camera_direction[0] +
2908  (point[1] - camera_position[1]) * camera_direction[1] +
2909  (point[2] - camera_position[2]) * camera_direction[2];
2910 
2911  Point<3> projection;
2912  projection[0] =
2913  camera_position[0] + phi * (point[0] - camera_position[0]);
2914  projection[1] =
2915  camera_position[1] + phi * (point[1] - camera_position[1]);
2916  projection[2] =
2917  camera_position[2] + phi * (point[2] - camera_position[2]);
2918 
2919  Point<2> projection_decomposition;
2920  projection_decomposition[0] = (projection[0] - camera_position[0] -
2921  camera_focus * camera_direction[0]) *
2922  camera_horizontal[0];
2923  projection_decomposition[0] += (projection[1] - camera_position[1] -
2924  camera_focus * camera_direction[1]) *
2925  camera_horizontal[1];
2926  projection_decomposition[0] += (projection[2] - camera_position[2] -
2927  camera_focus * camera_direction[2]) *
2928  camera_horizontal[2];
2929 
2930  projection_decomposition[1] = (projection[0] - camera_position[0] -
2931  camera_focus * camera_direction[0]) *
2932  camera_vertical[0];
2933  projection_decomposition[1] += (projection[1] - camera_position[1] -
2934  camera_focus * camera_direction[1]) *
2935  camera_vertical[1];
2936  projection_decomposition[1] += (projection[2] - camera_position[2] -
2937  camera_focus * camera_direction[2]) *
2938  camera_vertical[2];
2939 
2940  return projection_decomposition;
2941  }
2942 
2943 
2948  Point<6>
2949  svg_get_gradient_parameters(Point<3> points[])
2950  {
2951  Point<3> v_min, v_max, v_inter;
2952 
2953  // Use the Bubblesort algorithm to sort the points with respect to the
2954  // third coordinate
2955  for (int i = 0; i < 2; ++i)
2956  {
2957  for (int j = 0; j < 2 - i; ++j)
2958  {
2959  if (points[j][2] > points[j + 1][2])
2960  {
2961  Point<3> temp = points[j];
2962  points[j] = points[j + 1];
2963  points[j + 1] = temp;
2964  }
2965  }
2966  }
2967 
2968  // save the related three-dimensional vectors v_min, v_inter, and v_max
2969  v_min = points[0];
2970  v_inter = points[1];
2971  v_max = points[2];
2972 
2973  Point<2> A[2];
2974  Point<2> b, gradient;
2975 
2976  // determine the plane offset c
2977  A[0][0] = v_max[0] - v_min[0];
2978  A[0][1] = v_inter[0] - v_min[0];
2979  A[1][0] = v_max[1] - v_min[1];
2980  A[1][1] = v_inter[1] - v_min[1];
2981 
2982  b[0] = -v_min[0];
2983  b[1] = -v_min[1];
2984 
2985  double x, sum;
2986  bool col_change = false;
2987 
2988  if (A[0][0] == 0)
2989  {
2990  col_change = true;
2991 
2992  A[0][0] = A[0][1];
2993  A[0][1] = 0;
2994 
2995  double temp = A[1][0];
2996  A[1][0] = A[1][1];
2997  A[1][1] = temp;
2998  }
2999 
3000  for (unsigned int k = 0; k < 1; ++k)
3001  {
3002  for (unsigned int i = k + 1; i < 2; ++i)
3003  {
3004  x = A[i][k] / A[k][k];
3005 
3006  for (unsigned int j = k + 1; j < 2; ++j)
3007  A[i][j] = A[i][j] - A[k][j] * x;
3008 
3009  b[i] = b[i] - b[k] * x;
3010  }
3011  }
3012 
3013  b[1] = b[1] / A[1][1];
3014 
3015  for (int i = 0; i >= 0; i--)
3016  {
3017  sum = b[i];
3018 
3019  for (unsigned int j = i + 1; j < 2; ++j)
3020  sum = sum - A[i][j] * b[j];
3021 
3022  b[i] = sum / A[i][i];
3023  }
3024 
3025  if (col_change)
3026  {
3027  double temp = b[0];
3028  b[0] = b[1];
3029  b[1] = temp;
3030  }
3031 
3032  double c = b[0] * (v_max[2] - v_min[2]) + b[1] * (v_inter[2] - v_min[2]) +
3033  v_min[2];
3034 
3035  // Determine the first entry of the gradient (phi, cf. documentation)
3036  A[0][0] = v_max[0] - v_min[0];
3037  A[0][1] = v_inter[0] - v_min[0];
3038  A[1][0] = v_max[1] - v_min[1];
3039  A[1][1] = v_inter[1] - v_min[1];
3040 
3041  b[0] = 1.0 - v_min[0];
3042  b[1] = -v_min[1];
3043 
3044  col_change = false;
3046  if (A[0][0] == 0)
3047  {
3048  col_change = true;
3049 
3050  A[0][0] = A[0][1];
3051  A[0][1] = 0;
3052 
3053  double temp = A[1][0];
3054  A[1][0] = A[1][1];
3055  A[1][1] = temp;
3056  }
3057 
3058  for (unsigned int k = 0; k < 1; ++k)
3059  {
3060  for (unsigned int i = k + 1; i < 2; ++i)
3061  {
3062  x = A[i][k] / A[k][k];
3063 
3064  for (unsigned int j = k + 1; j < 2; ++j)
3065  A[i][j] = A[i][j] - A[k][j] * x;
3066 
3067  b[i] = b[i] - b[k] * x;
3068  }
3069  }
3070 
3071  b[1] = b[1] / A[1][1];
3072 
3073  for (int i = 0; i >= 0; i--)
3074  {
3075  sum = b[i];
3076 
3077  for (unsigned int j = i + 1; j < 2; ++j)
3078  sum = sum - A[i][j] * b[j];
3079 
3080  b[i] = sum / A[i][i];
3081  }
3082 
3083  if (col_change)
3084  {
3085  double temp = b[0];
3086  b[0] = b[1];
3087  b[1] = temp;
3088  }
3089 
3090  gradient[0] = b[0] * (v_max[2] - v_min[2]) +
3091  b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3092 
3093  // determine the second entry of the gradient
3094  A[0][0] = v_max[0] - v_min[0];
3095  A[0][1] = v_inter[0] - v_min[0];
3096  A[1][0] = v_max[1] - v_min[1];
3097  A[1][1] = v_inter[1] - v_min[1];
3098 
3099  b[0] = -v_min[0];
3100  b[1] = 1.0 - v_min[1];
3101 
3102  col_change = false;
3103 
3104  if (A[0][0] == 0)
3105  {
3106  col_change = true;
3107 
3108  A[0][0] = A[0][1];
3109  A[0][1] = 0;
3110 
3111  double temp = A[1][0];
3112  A[1][0] = A[1][1];
3113  A[1][1] = temp;
3114  }
3115 
3116  for (unsigned int k = 0; k < 1; ++k)
3117  {
3118  for (unsigned int i = k + 1; i < 2; ++i)
3119  {
3120  x = A[i][k] / A[k][k];
3121 
3122  for (unsigned int j = k + 1; j < 2; ++j)
3123  A[i][j] = A[i][j] - A[k][j] * x;
3124 
3125  b[i] = b[i] - b[k] * x;
3126  }
3127  }
3128 
3129  b[1] = b[1] / A[1][1];
3130 
3131  for (int i = 0; i >= 0; i--)
3132  {
3133  sum = b[i];
3134 
3135  for (unsigned int j = i + 1; j < 2; ++j)
3136  sum = sum - A[i][j] * b[j];
3137 
3138  b[i] = sum / A[i][i];
3139  }
3140 
3141  if (col_change)
3142  {
3143  double temp = b[0];
3144  b[0] = b[1];
3145  b[1] = temp;
3146  }
3147 
3148  gradient[1] = b[0] * (v_max[2] - v_min[2]) +
3149  b[1] * (v_inter[2] - v_min[2]) - c + v_min[2];
3150 
3151  // normalize the gradient
3152  gradient /= gradient.norm();
3153 
3154  const double lambda = -gradient[0] * (v_min[0] - v_max[0]) -
3155  gradient[1] * (v_min[1] - v_max[1]);
3156 
3157  Point<6> gradient_parameters;
3158 
3159  gradient_parameters[0] = v_min[0];
3160  gradient_parameters[1] = v_min[1];
3161 
3162  gradient_parameters[2] = v_min[0] + lambda * gradient[0];
3163  gradient_parameters[3] = v_min[1] + lambda * gradient[1];
3164 
3165  gradient_parameters[4] = v_min[2];
3166  gradient_parameters[5] = v_max[2];
3167 
3168  return gradient_parameters;
3169  }
3170  } // namespace
3171 
3172 
3173 
3174  template <int dim, int spacedim>
3175  void
3177  const std::vector<Patch<dim, spacedim>> &patches,
3178  const std::vector<std::string> &data_names,
3179  const std::vector<
3180  std::tuple<unsigned int,
3181  unsigned int,
3182  std::string,
3184  const UcdFlags &flags,
3185  std::ostream &out)
3186  {
3187  // Note that while in theory dim==0 should be implemented, this is not
3188  // tested, therefore currently not allowed.
3189  AssertThrow(dim > 0, ExcNotImplemented());
3190 
3191  AssertThrow(out.fail() == false, ExcIO());
3192 
3193 #ifndef DEAL_II_WITH_MPI
3194  // verify that there are indeed patches to be written out. most of the
3195  // times, people just forget to call build_patches when there are no
3196  // patches, so a warning is in order. that said, the assertion is disabled
3197  // if we support MPI since then it can happen that on the coarsest mesh, a
3198  // processor simply has no cells it actually owns, and in that case it is
3199  // legit if there are no patches
3200  Assert(patches.size() > 0, ExcNoPatches());
3201 #else
3202  if (patches.empty())
3203  return;
3204 #endif
3205 
3206  const unsigned int n_data_sets = data_names.size();
3207 
3208  UcdStream ucd_out(out, flags);
3209 
3210  // first count the number of cells and cells for later use
3211  unsigned int n_nodes;
3212  unsigned int n_cells;
3213  std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
3214  //---------------------
3215  // preamble
3216  if (flags.write_preamble)
3217  {
3218  out
3219  << "# This file was generated by the deal.II library." << '\n'
3220  << "# Date = " << Utilities::System::get_date() << '\n'
3221  << "# Time = " << Utilities::System::get_time() << '\n'
3222  << "#" << '\n'
3223  << "# For a description of the UCD format see the AVS Developer's guide."
3224  << '\n'
3225  << "#" << '\n';
3226  }
3227 
3228  // start with ucd data
3229  out << n_nodes << ' ' << n_cells << ' ' << n_data_sets << ' ' << 0
3230  << ' ' // no cell data at present
3231  << 0 // no model data
3232  << '\n';
3233 
3234  write_nodes(patches, ucd_out);
3235  out << '\n';
3236 
3237  write_cells(patches, ucd_out);
3238  out << '\n';
3239 
3240  //---------------------------
3241  // now write data
3242  if (n_data_sets != 0)
3243  {
3244  out << n_data_sets << " "; // number of vectors
3245  for (unsigned int i = 0; i < n_data_sets; ++i)
3246  out << 1 << ' '; // number of components;
3247  // only 1 supported presently
3248  out << '\n';
3249 
3250  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3251  out << data_names[data_set]
3252  << ",dimensionless" // no units supported at present
3253  << '\n';
3254 
3255  write_data(patches, n_data_sets, true, ucd_out);
3256  }
3257  // make sure everything now gets to disk
3258  out.flush();
3259 
3260  // assert the stream is still ok
3261  AssertThrow(out.fail() == false, ExcIO());
3262  }
3263 
3264 
3265  template <int dim, int spacedim>
3266  void
3268  const std::vector<Patch<dim, spacedim>> &patches,
3269  const std::vector<std::string> &data_names,
3270  const std::vector<
3271  std::tuple<unsigned int,
3272  unsigned int,
3273  std::string,
3275  const DXFlags &flags,
3276  std::ostream &out)
3277  {
3278  // Point output is currently not implemented.
3279  AssertThrow(dim > 0, ExcNotImplemented());
3280 
3281  AssertThrow(out.fail() == false, ExcIO());
3282 
3283 #ifndef DEAL_II_WITH_MPI
3284  // verify that there are indeed patches to be written out. most of the
3285  // times, people just forget to call build_patches when there are no
3286  // patches, so a warning is in order. that said, the assertion is disabled
3287  // if we support MPI since then it can happen that on the coarsest mesh, a
3288  // processor simply has no cells it actually owns, and in that case it is
3289  // legit if there are no patches
3290  Assert(patches.size() > 0, ExcNoPatches());
3291 #else
3292  if (patches.empty())
3293  return;
3294 #endif
3295  // Stream with special features for dx output
3296  DXStream dx_out(out, flags);
3297 
3298  // Variable counting the offset of binary data.
3299  unsigned int offset = 0;
3300 
3301  const unsigned int n_data_sets = data_names.size();
3302 
3303  // first count the number of cells and cells for later use
3304  unsigned int n_nodes;
3305  unsigned int n_cells;
3306  std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
3307 
3308  // start with vertices order is lexicographical, x varying fastest
3309  out << "object \"vertices\" class array type float rank 1 shape "
3310  << spacedim << " items " << n_nodes;
3311 
3312  if (flags.coordinates_binary)
3313  {
3314  out << " lsb ieee data 0" << '\n';
3315  offset += n_nodes * spacedim * sizeof(float);
3316  }
3317  else
3318  {
3319  out << " data follows" << '\n';
3320  write_nodes(patches, dx_out);
3321  }
3322 
3323  //-----------------------------
3324  // first write the coordinates of all vertices
3325 
3326  //---------------------------------------
3327  // write cells
3328  out << "object \"cells\" class array type int rank 1 shape "
3329  << GeometryInfo<dim>::vertices_per_cell << " items " << n_cells;
3330 
3331  if (flags.int_binary)
3332  {
3333  out << " lsb binary data " << offset << '\n';
3334  offset += n_cells * sizeof(int);
3335  }
3336  else
3337  {
3338  out << " data follows" << '\n';
3339  write_cells(patches, dx_out);
3340  out << '\n';
3341  }
3342 
3343 
3344  out << "attribute \"element type\" string \"";
3345  if (dim == 1)
3346  out << "lines";
3347  if (dim == 2)
3348  out << "quads";
3349  if (dim == 3)
3350  out << "cubes";
3351  out << "\"" << '\n' << "attribute \"ref\" string \"positions\"" << '\n';
3352 
3353  // TODO:[GK] Patches must be of same size!
3354  //---------------------------
3355  // write neighbor information
3356  if (flags.write_neighbors)
3357  {
3358  out << "object \"neighbors\" class array type int rank 1 shape "
3359  << GeometryInfo<dim>::faces_per_cell << " items " << n_cells
3360  << " data follows";
3361 
3362  for (const auto &patch : patches)
3363  {
3364  const unsigned int n = patch.n_subdivisions;
3365  const unsigned int n1 = (dim > 0) ? n : 1;
3366  const unsigned int n2 = (dim > 1) ? n : 1;
3367  const unsigned int n3 = (dim > 2) ? n : 1;
3368  const unsigned int x_minus = (dim > 0) ? 0 : 0;
3369  const unsigned int x_plus = (dim > 0) ? 1 : 0;
3370  const unsigned int y_minus = (dim > 1) ? 2 : 0;
3371  const unsigned int y_plus = (dim > 1) ? 3 : 0;
3372  const unsigned int z_minus = (dim > 2) ? 4 : 0;
3373  const unsigned int z_plus = (dim > 2) ? 5 : 0;
3374  unsigned int cells_per_patch = Utilities::fixed_power<dim>(n);
3375  unsigned int dx = 1;
3376  unsigned int dy = n;
3377  unsigned int dz = n * n;
3378 
3379  const unsigned int patch_start =
3380  patch.patch_index * cells_per_patch;
3381 
3382  for (unsigned int i3 = 0; i3 < n3; ++i3)
3383  for (unsigned int i2 = 0; i2 < n2; ++i2)
3384  for (unsigned int i1 = 0; i1 < n1; ++i1)
3385  {
3386  const unsigned int nx = i1 * dx;
3387  const unsigned int ny = i2 * dy;
3388  const unsigned int nz = i3 * dz;
3389 
3390  // There are no neighbors for dim==0. Note that this case is
3391  // caught by the AssertThrow at the beginning of this
3392  // function anyway. This condition avoids compiler warnings.
3393  if (dim < 1)
3394  continue;
3395 
3396  out << '\n';
3397  // Direction -x Last cell in row of other patch
3398  if (i1 == 0)
3399  {
3400  const unsigned int nn = patch.neighbors[x_minus];
3401  out << '\t';
3402  if (nn != patch.no_neighbor)
3403  out
3404  << (nn * cells_per_patch + ny + nz + dx * (n - 1));
3405  else
3406  out << "-1";
3407  }
3408  else
3409  {
3410  out << '\t' << patch_start + nx - dx + ny + nz;
3411  }
3412  // Direction +x First cell in row of other patch
3413  if (i1 == n - 1)
3414  {
3415  const unsigned int nn = patch.neighbors[x_plus];
3416  out << '\t';
3417  if (nn != patch.no_neighbor)
3418  out << (nn * cells_per_patch + ny + nz);
3419  else
3420  out << "-1";
3421  }
3422  else
3423  {
3424  out << '\t' << patch_start + nx + dx + ny + nz;
3425  }
3426  if (dim < 2)
3427  continue;
3428  // Direction -y
3429  if (i2 == 0)
3430  {
3431  const unsigned int nn = patch.neighbors[y_minus];
3432  out << '\t';
3433  if (nn != patch.no_neighbor)
3434  out
3435  << (nn * cells_per_patch + nx + nz + dy * (n - 1));
3436  else
3437  out << "-1";
3438  }
3439  else
3440  {
3441  out << '\t' << patch_start + nx + ny - dy + nz;
3442  }
3443  // Direction +y
3444  if (i2 == n - 1)
3445  {
3446  const unsigned int nn = patch.neighbors[y_plus];
3447  out << '\t';
3448  if (nn != patch.no_neighbor)
3449  out << (nn * cells_per_patch + nx + nz);
3450  else
3451  out << "-1";
3452  }
3453  else
3454  {
3455  out << '\t' << patch_start + nx + ny + dy + nz;
3456  }
3457  if (dim < 3)
3458  continue;
3459 
3460  // Direction -z
3461  if (i3 == 0)
3462  {
3463  const unsigned int nn = patch.neighbors[z_minus];
3464  out << '\t';
3465  if (nn != patch.no_neighbor)
3466  out
3467  << (nn * cells_per_patch + nx + ny + dz * (n - 1));
3468  else
3469  out << "-1";
3470  }
3471  else
3472  {
3473  out << '\t' << patch_start + nx + ny + nz - dz;
3474  }
3475  // Direction +z
3476  if (i3 == n - 1)
3477  {
3478  const unsigned int nn = patch.neighbors[z_plus];
3479  out << '\t';
3480  if (nn != patch.no_neighbor)
3481  out << (nn * cells_per_patch + nx + ny);
3482  else
3483  out << "-1";
3484  }
3485  else
3486  {
3487  out << '\t' << patch_start + nx + ny + nz + dz;
3488  }
3489  }
3490  out << '\n';
3491  }
3492  }
3493  //---------------------------
3494  // now write data
3495  if (n_data_sets != 0)
3496  {
3497  out << "object \"data\" class array type float rank 1 shape "
3498  << n_data_sets << " items " << n_nodes;
3499 
3500  if (flags.data_binary)
3501  {
3502  out << " lsb ieee data " << offset << '\n';
3503  offset += n_data_sets * n_nodes *
3504  ((flags.data_double) ? sizeof(double) : sizeof(float));
3505  }
3506  else
3507  {
3508  out << " data follows" << '\n';
3509  write_data(patches, n_data_sets, flags.data_double, dx_out);
3510  }
3511 
3512  // loop over all patches
3513  out << "attribute \"dep\" string \"positions\"" << '\n';
3514  }
3515  else
3516  {
3517  out << "object \"data\" class constantarray type float rank 0 items "
3518  << n_nodes << " data follows" << '\n'
3519  << '0' << '\n';
3520  }
3521 
3522  // no model data
3523 
3524  out << "object \"deal data\" class field" << '\n'
3525  << "component \"positions\" value \"vertices\"" << '\n'
3526  << "component \"connections\" value \"cells\"" << '\n'
3527  << "component \"data\" value \"data\"" << '\n';
3528 
3529  if (flags.write_neighbors)
3530  out << "component \"neighbors\" value \"neighbors\"" << '\n';
3531 
3532  {
3533  out << "attribute \"created\" string \"" << Utilities::System::get_date()
3534  << ' ' << Utilities::System::get_time() << '"' << '\n';
3535  }
3536 
3537  out << "end" << '\n';
3538  // Write all binary data now
3539  if (flags.coordinates_binary)
3540  write_nodes(patches, dx_out);
3541  if (flags.int_binary)
3542  write_cells(patches, dx_out);
3543  if (flags.data_binary)
3544  write_data(patches, n_data_sets, flags.data_double, dx_out);
3545 
3546  // make sure everything now gets to disk
3547  out.flush();
3548 
3549  // assert the stream is still ok
3550  AssertThrow(out.fail() == false, ExcIO());
3551  }
3552 
3553 
3554 
3555  template <int dim, int spacedim>
3556  void
3558  const std::vector<Patch<dim, spacedim>> &patches,
3559  const std::vector<std::string> &data_names,
3560  const std::vector<
3561  std::tuple<unsigned int,
3562  unsigned int,
3563  std::string,
3565  const GnuplotFlags &flags,
3566  std::ostream &out)
3567  {
3568  AssertThrow(out.fail() == false, ExcIO());
3569 
3570 #ifndef DEAL_II_WITH_MPI
3571  // verify that there are indeed patches to be written out. most
3572  // of the times, people just forget to call build_patches when there
3573  // are no patches, so a warning is in order. that said, the
3574  // assertion is disabled if we support MPI since then it can
3575  // happen that on the coarsest mesh, a processor simply has no
3576  // cells it actually owns, and in that case it is legit if there
3577  // are no patches
3578  Assert(patches.size() > 0, ExcNoPatches());
3579 #else
3580  if (patches.empty())
3581  return;
3582 #endif
3583 
3584  const unsigned int n_data_sets = data_names.size();
3585 
3586  // write preamble
3587  {
3588  out << "# This file was generated by the deal.II library." << '\n'
3589  << "# Date = " << Utilities::System::get_date() << '\n'
3590  << "# Time = " << Utilities::System::get_time() << '\n'
3591  << "#" << '\n'
3592  << "# For a description of the GNUPLOT format see the GNUPLOT manual."
3593  << '\n'
3594  << "#" << '\n'
3595  << "# ";
3596 
3597  AssertThrow(spacedim <= flags.space_dimension_labels.size(),
3599  for (unsigned int spacedim_n = 0; spacedim_n < spacedim; ++spacedim_n)
3600  {
3601  out << '<' << flags.space_dimension_labels.at(spacedim_n) << "> ";
3602  }
3603 
3604  for (const auto &data_name : data_names)
3605  out << '<' << data_name << "> ";
3606  out << '\n';
3607  }
3608 
3609 
3610  // loop over all patches
3611  for (const auto &patch : patches)
3612  {
3613  const unsigned int n_subdivisions = patch.n_subdivisions;
3614  const unsigned int n_points_per_direction = n_subdivisions + 1;
3615 
3616  Assert((patch.data.n_rows() == n_data_sets &&
3617  !patch.points_are_available) ||
3618  (patch.data.n_rows() == n_data_sets + spacedim &&
3619  patch.points_are_available),
3621  (n_data_sets + spacedim) :
3622  n_data_sets,
3623  patch.data.n_rows()));
3624 
3625  auto output_point_data =
3626  [&out, &patch, n_data_sets](const unsigned int point_index) mutable {
3627  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
3628  out << patch.data(data_set, point_index) << ' ';
3629  };
3630 
3631  switch (dim)
3632  {
3633  case 0:
3634  {
3635  Assert(patch.reference_cell == ReferenceCells::Vertex,
3636  ExcInternalError());
3637  Assert(patch.data.n_cols() == 1,
3638  ExcInvalidDatasetSize(patch.data.n_cols(),
3639  n_subdivisions + 1));
3640 
3641 
3642  // compute coordinates for this patch point
3643  out << get_equispaced_location(patch, {}, n_subdivisions)
3644  << ' ';
3645  output_point_data(0);
3646  out << '\n';
3647  out << '\n';
3648  break;
3649  }
3650 
3651  case 1:
3652  {
3653  Assert(patch.reference_cell == ReferenceCells::Line,
3654  ExcInternalError());
3655  Assert(patch.data.n_cols() ==
3656  Utilities::fixed_power<dim>(n_points_per_direction),
3657  ExcInvalidDatasetSize(patch.data.n_cols(),
3658  n_subdivisions + 1));
3659 
3660  for (unsigned int i1 = 0; i1 < n_points_per_direction; ++i1)
3661  {
3662  // compute coordinates for this patch point
3663  out << get_equispaced_location(patch, {i1}, n_subdivisions)
3664  << ' ';
3665 
3666  output_point_data(i1);
3667  out << '\n';
3668  }
3669  // end of patch
3670  out << '\n';
3671  out << '\n';
3672  break;
3673  }
3674 
3675  case 2:
3676  {
3677  if (patch.reference_cell == ReferenceCells::Quadrilateral)
3678  {
3679  Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3680  n_points_per_direction),
3681  ExcInvalidDatasetSize(patch.data.n_cols(),
3682  n_subdivisions + 1));
3683 
3684  for (unsigned int i2 = 0; i2 < n_points_per_direction; ++i2)
3685  {
3686  for (unsigned int i1 = 0; i1 < n_points_per_direction;
3687  ++i1)
3688  {
3689  // compute coordinates for this patch point
3690  out << get_equispaced_location(patch,
3691  {i1, i2},
3692  n_subdivisions)
3693  << ' ';
3694 
3695  output_point_data(i1 + i2 * n_points_per_direction);
3696  out << '\n';
3697  }
3698  // end of row in patch
3699  out << '\n';
3700  }
3701  }
3702  else if (patch.reference_cell == ReferenceCells::Triangle)
3703  {
3704  Assert(n_subdivisions == 1, ExcNotImplemented());
3705 
3706  Assert(patch.data.n_cols() == 3, ExcInternalError());
3707 
3708  // Gnuplot can only plot surfaces if each facet of the
3709  // surface is a bilinear patch, or a subdivided bilinear
3710  // patch with equally many points along each row of the
3711  // subdivision. This is what the code above for
3712  // quadrilaterals does. We emulate this by repeating the
3713  // third point of a triangle twice so that there are two
3714  // points for that row as well -- i.e., we write a 2x2
3715  // bilinear patch where two of the points are collapsed onto
3716  // one vertex.
3717  //
3718  // This also matches the example here:
3719  // https://stackoverflow.com/questions/42784369/drawing-triangular-mesh-using-gnuplot
3720  out << get_node_location(patch, 0) << ' ';
3721  output_point_data(0);
3722  out << '\n';
3723 
3724  out << get_node_location(patch, 1) << ' ';
3725  output_point_data(1);
3726  out << '\n';
3727  out << '\n'; // end of one row of points
3728 
3729  out << get_node_location(patch, 2) << ' ';
3730  output_point_data(2);
3731  out << '\n';
3732 
3733  out << get_node_location(patch, 2) << ' ';
3734  output_point_data(2);
3735  out << '\n';
3736  out << '\n'; // end of the second row of points
3737  out << '\n'; // end of the entire patch
3738  }
3739  else
3740  // There aren't any other reference cells in 2d than the
3741  // quadrilateral and the triangle. So whatever we got here
3742  // can't be any good
3743  Assert(false, ExcInternalError());
3744  // end of patch
3745  out << '\n';
3746 
3747  break;
3748  }
3749 
3750  case 3:
3751  {
3752  if (patch.reference_cell == ReferenceCells::Hexahedron)
3753  {
3754  Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(
3755  n_points_per_direction),
3756  ExcInvalidDatasetSize(patch.data.n_cols(),
3757  n_subdivisions + 1));
3758 
3759  // for all grid points: draw lines into all positive
3760  // coordinate directions if there is another grid point
3761  // there
3762  for (unsigned int i3 = 0; i3 < n_points_per_direction; ++i3)
3763  for (unsigned int i2 = 0; i2 < n_points_per_direction;
3764  ++i2)
3765  for (unsigned int i1 = 0; i1 < n_points_per_direction;
3766  ++i1)
3767  {
3768  // compute coordinates for this patch point
3769  const Point<spacedim> this_point =
3770  get_equispaced_location(patch,
3771  {i1, i2, i3},
3772  n_subdivisions);
3773  // line into positive x-direction if possible
3774  if (i1 < n_subdivisions)
3775  {
3776  // write point here and its data
3777  out << this_point << ' ';
3778  output_point_data(i1 +
3779  i2 * n_points_per_direction +
3780  i3 * n_points_per_direction *
3781  n_points_per_direction);
3782  out << '\n';
3783 
3784  // write point there and its data
3785  out << get_equispaced_location(patch,
3786  {i1 + 1, i2, i3},
3787  n_subdivisions)
3788  << ' ';
3789 
3790  output_point_data((i1 + 1) +
3791  i2 * n_points_per_direction +
3792  i3 * n_points_per_direction *
3793  n_points_per_direction);
3794  out << '\n';
3795 
3796  // end of line
3797  out << '\n' << '\n';
3798  }
3799 
3800  // line into positive y-direction if possible
3801  if (i2 < n_subdivisions)
3802  {
3803  // write point here and its data
3804  out << this_point << ' ';
3805  output_point_data(i1 +
3806  i2 * n_points_per_direction +
3807  i3 * n_points_per_direction *
3808  n_points_per_direction);
3809  out << '\n';
3810 
3811  // write point there and its data
3812  out << get_equispaced_location(patch,
3813  {i1, i2 + 1, i3},
3814  n_subdivisions)
3815  << ' ';
3816 
3817  output_point_data(
3818  i1 + (i2 + 1) * n_points_per_direction +
3819  i3 * n_points_per_direction *
3820  n_points_per_direction);
3821  out << '\n';
3822 
3823  // end of line
3824  out << '\n' << '\n';
3825  }
3826 
3827  // line into positive z-direction if possible
3828  if (i3 < n_subdivisions)
3829  {
3830  // write point here and its data
3831  out << this_point << ' ';
3832  output_point_data(i1 +
3833  i2 * n_points_per_direction +
3834  i3 * n_points_per_direction *
3835  n_points_per_direction);
3836  out << '\n';
3837 
3838  // write point there and its data
3839  out << get_equispaced_location(patch,
3840  {i1, i2, i3 + 1},
3841  n_subdivisions)
3842  << ' ';
3843 
3844  output_point_data(
3845  i1 + i2 * n_points_per_direction +
3846  (i3 + 1) * n_points_per_direction *
3847  n_points_per_direction);
3848  out << '\n';
3849  // end of line
3850  out << '\n' << '\n';
3851  }
3852  }
3853  }
3854  else if (patch.reference_cell == ReferenceCells::Tetrahedron)
3855  {
3856  Assert(n_subdivisions == 1, ExcNotImplemented());
3857 
3858  // Draw the tetrahedron as a collection of two lines.
3859  for (const unsigned int v : {0, 1, 2, 0, 3, 2})
3860  {
3861  out << get_node_location(patch, v) << ' ';
3862  output_point_data(v);
3863  out << '\n';
3864  }
3865  out << '\n'; // end of first line
3866 
3867  for (const unsigned int v : {3, 1})
3868  {
3869  out << get_node_location(patch, v) << ' ';
3870  output_point_data(v);
3871  out << '\n';
3872  }
3873  out << '\n'; // end of second line
3874  }
3875  else if (patch.reference_cell == ReferenceCells::Pyramid)
3876  {
3877  Assert(n_subdivisions == 1, ExcNotImplemented());
3878 
3879  // Draw the pyramid as a collection of two lines.
3880  for (const unsigned int v : {0, 1, 3, 2, 0, 4, 1})
3881  {
3882  out << get_node_location(patch, v) << ' ';
3883  output_point_data(v);
3884  out << '\n';
3885  }
3886  out << '\n'; // end of first line
3887 
3888  for (const unsigned int v : {2, 4, 3})
3889  {
3890  out << get_node_location(patch, v) << ' ';
3891  output_point_data(v);
3892  out << '\n';
3893  }
3894  out << '\n'; // end of second line
3895  }
3896  else if (patch.reference_cell == ReferenceCells::Wedge)
3897  {
3898  Assert(n_subdivisions == 1, ExcNotImplemented());
3899 
3900  // Draw the wedge as a collection of three
3901  // lines. The first one wraps around the base,
3902  // goes up to the top, and wraps around that. The
3903  // second and third are just individual lines
3904  // going from base to top.
3905  for (const unsigned int v : {0, 1, 2, 0, 3, 4, 5, 3})
3906  {
3907  out << get_node_location(patch, v) << ' ';
3908  output_point_data(v);
3909  out << '\n';
3910  }
3911  out << '\n'; // end of first line
3912 
3913  for (const unsigned int v : {1, 4})
3914  {
3915  out << get_node_location(patch, v) << ' ';
3916  output_point_data(v);
3917  out << '\n';
3918  }
3919  out << '\n'; // end of second line
3920 
3921  for (const unsigned int v : {2, 5})
3922  {
3923  out << get_node_location(patch, v) << ' ';
3924  output_point_data(v);
3925  out << '\n';
3926  }
3927  out << '\n'; // end of second line
3928  }
3929  else
3930  // No other reference cells are currently implemented
3931  Assert(false, ExcNotImplemented());
3932 
3933  break;
3934  }
3935 
3936  default:
3937  Assert(false, ExcNotImplemented());
3938  }
3939  }
3940  // make sure everything now gets to disk
3941  out.flush();
3942 
3943  AssertThrow(out.fail() == false, ExcIO());
3944  }
3945 
3946 
3947  namespace
3948  {
3949  template <int dim, int spacedim>
3950  void
3951  do_write_povray(const std::vector<Patch<dim, spacedim>> &,
3952  const std::vector<std::string> &,
3953  const PovrayFlags &,
3954  std::ostream &)
3955  {
3956  Assert(false,
3957  ExcMessage("Writing files in POVRAY format is only supported "
3958  "for two-dimensional meshes."));
3959  }
3960 
3961 
3962 
3963  void
3964  do_write_povray(const std::vector<Patch<2, 2>> &patches,
3965  const std::vector<std::string> &data_names,
3966  const PovrayFlags &flags,
3967  std::ostream &out)
3968  {
3969  AssertThrow(out.fail() == false, ExcIO());
3970 
3971 #ifndef DEAL_II_WITH_MPI
3972  // verify that there are indeed patches to be written out. most
3973  // of the times, people just forget to call build_patches when there
3974  // are no patches, so a warning is in order. that said, the
3975  // assertion is disabled if we support MPI since then it can
3976  // happen that on the coarsest mesh, a processor simply has no cells it
3977  // actually owns, and in that case it is legit if there are no patches
3978  Assert(patches.size() > 0, ExcNoPatches());
3979 #else
3980  if (patches.empty())
3981  return;
3982 #endif
3983  constexpr int dim = 2;
3984  (void)dim;
3985  constexpr int spacedim = 2;
3986 
3987  const unsigned int n_data_sets = data_names.size();
3988  (void)n_data_sets;
3989 
3990  // write preamble
3991  {
3992  out
3993  << "/* This file was generated by the deal.II library." << '\n'
3994  << " Date = " << Utilities::System::get_date() << '\n'
3995  << " Time = " << Utilities::System::get_time() << '\n'
3996  << '\n'
3997  << " For a description of the POVRAY format see the POVRAY manual."
3998  << '\n'
3999  << "*/ " << '\n';
4000 
4001  // include files
4002  out << "#include \"colors.inc\" " << '\n'
4003  << "#include \"textures.inc\" " << '\n';
4004 
4005 
4006  // use external include file for textures, camera and light
4007  if (flags.external_data)
4008  out << "#include \"data.inc\" " << '\n';
4009  else // all definitions in data file
4010  {
4011  // camera
4012  out << '\n'
4013  << '\n'
4014  << "camera {" << '\n'
4015  << " location <1,4,-7>" << '\n'
4016  << " look_at <0,0,0>" << '\n'
4017  << " angle 30" << '\n'
4018  << "}" << '\n';
4019 
4020  // light
4021  out << '\n'
4022  << "light_source {" << '\n'
4023  << " <1,4,-7>" << '\n'
4024  << " color Grey" << '\n'
4025  << "}" << '\n';
4026  out << '\n'
4027  << "light_source {" << '\n'
4028  << " <0,20,0>" << '\n'
4029  << " color White" << '\n'
4030  << "}" << '\n';
4031  }
4032  }
4033 
4034  // max. and min. height of solution
4035  Assert(patches.size() > 0, ExcNoPatches());
4036  double hmin = patches[0].data(0, 0);
4037  double hmax = patches[0].data(0, 0);
4038 
4039  for (const auto &patch : patches)
4040  {
4041  const unsigned int n_subdivisions = patch.n_subdivisions;
4042 
4043  Assert((patch.data.n_rows() == n_data_sets &&
4044  !patch.points_are_available) ||
4045  (patch.data.n_rows() == n_data_sets + spacedim &&
4046  patch.points_are_available),
4048  (n_data_sets + spacedim) :
4049  n_data_sets,
4050  patch.data.n_rows()));
4051  Assert(patch.data.n_cols() ==
4052  Utilities::fixed_power<dim>(n_subdivisions + 1),
4053  ExcInvalidDatasetSize(patch.data.n_cols(),
4054  n_subdivisions + 1));
4055 
4056  for (unsigned int i = 0; i < n_subdivisions + 1; ++i)
4057  for (unsigned int j = 0; j < n_subdivisions + 1; ++j)
4058  {
4059  const int dl = i * (n_subdivisions + 1) + j;
4060  if (patch.data(0, dl) < hmin)
4061  hmin = patch.data(0, dl);
4062  if (patch.data(0, dl) > hmax)
4063  hmax = patch.data(0, dl);
4064  }
4065  }
4066 
4067  out << "#declare HMIN=" << hmin << ";" << '\n'
4068  << "#declare HMAX=" << hmax << ";" << '\n'
4069  << '\n';
4070 
4071  if (!flags.external_data)
4072  {
4073  // texture with scaled niveau lines 10 lines in the surface
4074  out << "#declare Tex=texture{" << '\n'
4075  << " pigment {" << '\n'
4076  << " gradient y" << '\n'
4077  << " scale y*(HMAX-HMIN)*" << 0.1 << '\n'
4078  << " color_map {" << '\n'
4079  << " [0.00 color Light_Purple] " << '\n'
4080  << " [0.95 color Light_Purple] " << '\n'
4081  << " [1.00 color White] " << '\n'
4082  << "} } }" << '\n'
4083  << '\n';
4084  }
4085 
4086  if (!flags.bicubic_patch)
4087  {
4088  // start of mesh header
4089  out << '\n' << "mesh {" << '\n';
4090  }
4091 
4092  // loop over all patches
4093  for (const auto &patch : patches)
4094  {
4095  const unsigned int n_subdivisions = patch.n_subdivisions;
4096  const unsigned int n = n_subdivisions + 1;
4097  const unsigned int d1 = 1;
4098  const unsigned int d2 = n;
4099 
4100  Assert((patch.data.n_rows() == n_data_sets &&
4101  !patch.points_are_available) ||
4102  (patch.data.n_rows() == n_data_sets + spacedim &&
4103  patch.points_are_available),
4105  (n_data_sets + spacedim) :
4106  n_data_sets,
4107  patch.data.n_rows()));
4108  Assert(patch.data.n_cols() == Utilities::fixed_power<dim>(n),
4109  ExcInvalidDatasetSize(patch.data.n_cols(),
4110  n_subdivisions + 1));
4111 
4112 
4113  std::vector<Point<spacedim>> ver(n * n);
4114 
4115  for (unsigned int i2 = 0; i2 < n; ++i2)
4116  for (unsigned int i1 = 0; i1 < n; ++i1)
4117  {
4118  // compute coordinates for this patch point, storing in ver
4119  ver[i1 * d1 + i2 * d2] =
4120  get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4121  }
4122 
4123 
4124  if (!flags.bicubic_patch)
4125  {
4126  // approximate normal vectors in patch
4127  std::vector<Point<3>> nrml;
4128  // only if smooth triangles are used
4129  if (flags.smooth)
4130  {
4131  nrml.resize(n * n);
4132  // These are difference quotients of the surface
4133  // mapping. We take them symmetric inside the
4134  // patch and one-sided at the edges
4135  Point<3> h1, h2;
4136  // Now compute normals in every point
4137  for (unsigned int i = 0; i < n; ++i)
4138  for (unsigned int j = 0; j < n; ++j)
4139  {
4140  const unsigned int il = (i == 0) ? i : (i - 1);
4141  const unsigned int ir =
4142  (i == n_subdivisions) ? i : (i + 1);
4143  const unsigned int jl = (j == 0) ? j : (j - 1);
4144  const unsigned int jr =
4145  (j == n_subdivisions) ? j : (j + 1);
4146 
4147  h1(0) =
4148  ver[ir * d1 + j * d2](0) - ver[il * d1 + j * d2](0);
4149  h1(1) = patch.data(0, ir * d1 + j * d2) -
4150  patch.data(0, il * d1 + j * d2);
4151  h1(2) =
4152  ver[ir * d1 + j * d2](1) - ver[il * d1 + j * d2](1);
4153 
4154  h2(0) =
4155  ver[i * d1 + jr * d2](0) - ver[i * d1 + jl * d2](0);
4156  h2(1) = patch.data(0, i * d1 + jr * d2) -
4157  patch.data(0, i * d1 + jl * d2);
4158  h2(2) =
4159  ver[i * d1 + jr * d2](1) - ver[i * d1 + jl * d2](1);
4160 
4161  nrml[i * d1 + j * d2](0) =
4162  h1(1) * h2(2) - h1(2) * h2(1);
4163  nrml[i * d1 + j * d2](1) =
4164  h1(2) * h2(0) - h1(0) * h2(2);
4165  nrml[i * d1 + j * d2](2) =
4166  h1(0) * h2(1) - h1(1) * h2(0);
4167 
4168  // normalize Vector
4169  double norm =
4170  std::sqrt(std::pow(nrml[i * d1 + j * d2](0), 2.) +
4171  std::pow(nrml[i * d1 + j * d2](1), 2.) +
4172  std::pow(nrml[i * d1 + j * d2](2), 2.));
4173 
4174  if (nrml[i * d1 + j * d2](1) < 0)
4175  norm *= -1.;
4176 
4177  for (unsigned int k = 0; k < 3; ++k)
4178  nrml[i * d1 + j * d2](k) /= norm;
4179  }
4180  }
4181 
4182  // setting up triangles
4183  for (unsigned int i = 0; i < n_subdivisions; ++i)
4184  for (unsigned int j = 0; j < n_subdivisions; ++j)
4185  {
4186  // down/left vertex of triangle
4187  const int dl = i * d1 + j * d2;
4188  if (flags.smooth)
4189  {
4190  // writing smooth_triangles
4191 
4192  // down/right triangle
4193  out << "smooth_triangle {" << '\n'
4194  << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4195  << "," << ver[dl](1) << ">, <" << nrml[dl](0)
4196  << ", " << nrml[dl](1) << ", " << nrml[dl](2)
4197  << ">," << '\n';
4198  out << " \t<" << ver[dl + d1](0) << ","
4199  << patch.data(0, dl + d1) << "," << ver[dl + d1](1)
4200  << ">, <" << nrml[dl + d1](0) << ", "
4201  << nrml[dl + d1](1) << ", " << nrml[dl + d1](2)
4202  << ">," << '\n';
4203  out << "\t<" << ver[dl + d1 + d2](0) << ","
4204  << patch.data(0, dl + d1 + d2) << ","
4205  << ver[dl + d1 + d2](1) << ">, <"
4206  << nrml[dl + d1 + d2](0) << ", "
4207  << nrml[dl + d1 + d2](1) << ", "
4208  << nrml[dl + d1 + d2](2) << ">}" << '\n';
4209 
4210  // upper/left triangle
4211  out << "smooth_triangle {" << '\n'
4212  << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4213  << "," << ver[dl](1) << ">, <" << nrml[dl](0)
4214  << ", " << nrml[dl](1) << ", " << nrml[dl](2)
4215  << ">," << '\n';
4216  out << "\t<" << ver[dl + d1 + d2](0) << ","
4217  << patch.data(0, dl + d1 + d2) << ","
4218  << ver[dl + d1 + d2](1) << ">, <"
4219  << nrml[dl + d1 + d2](0) << ", "
4220  << nrml[dl + d1 + d2](1) << ", "
4221  << nrml[dl + d1 + d2](2) << ">," << '\n';
4222  out << "\t<" << ver[dl + d2](0) << ","
4223  << patch.data(0, dl + d2) << "," << ver[dl + d2](1)
4224  << ">, <" << nrml[dl + d2](0) << ", "
4225  << nrml[dl + d2](1) << ", " << nrml[dl + d2](2)
4226  << ">}" << '\n';
4227  }
4228  else
4229  {
4230  // writing standard triangles down/right triangle
4231  out << "triangle {" << '\n'
4232  << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4233  << "," << ver[dl](1) << ">," << '\n';
4234  out << "\t<" << ver[dl + d1](0) << ","
4235  << patch.data(0, dl + d1) << "," << ver[dl + d1](1)
4236  << ">," << '\n';
4237  out << "\t<" << ver[dl + d1 + d2](0) << ","
4238  << patch.data(0, dl + d1 + d2) << ","
4239  << ver[dl + d1 + d2](1) << ">}" << '\n';
4240 
4241  // upper/left triangle
4242  out << "triangle {" << '\n'
4243  << "\t<" << ver[dl](0) << "," << patch.data(0, dl)
4244  << "," << ver[dl](1) << ">," << '\n';
4245  out << "\t<" << ver[dl + d1 + d2](0) << ","
4246  << patch.data(0, dl + d1 + d2) << ","
4247  << ver[dl + d1 + d2](1) << ">," << '\n';
4248  out << "\t<" << ver[dl + d2](0) << ","
4249  << patch.data(0, dl + d2) << "," << ver[dl + d2](1)
4250  << ">}" << '\n';
4251  }
4252  }
4253  }
4254  else
4255  {
4256  // writing bicubic_patch
4257  Assert(n_subdivisions == 3,
4258  ExcDimensionMismatch(n_subdivisions, 3));
4259  out << '\n'
4260  << "bicubic_patch {" << '\n'
4261  << " type 0" << '\n'
4262  << " flatness 0" << '\n'
4263  << " u_steps 0" << '\n'
4264  << " v_steps 0" << '\n';
4265  for (int i = 0; i < 16; ++i)
4266  {
4267  out << "\t<" << ver[i](0) << "," << patch.data(0, i) << ","
4268  << ver[i](1) << ">";
4269  if (i != 15)
4270  out << ",";
4271  out << '\n';
4272  }
4273  out << " texture {Tex}" << '\n' << "}" << '\n';
4274  }
4275  }
4276 
4277  if (!flags.bicubic_patch)
4278  {
4279  // the end of the mesh
4280  out << " texture {Tex}" << '\n' << "}" << '\n' << '\n';
4281  }
4282 
4283  // make sure everything now gets to disk
4284  out.flush();
4285 
4286  AssertThrow(out.fail() == false, ExcIO());
4287  }
4288  } // namespace
4289 
4290 
4291 
4292  template <int dim, int spacedim>
4293  void
4295  const std::vector<Patch<dim, spacedim>> &patches,
4296  const std::vector<std::string> &data_names,
4297  const std::vector<
4298  std::tuple<unsigned int,
4299  unsigned int,
4300  std::string,
4302  const PovrayFlags &flags,
4303  std::ostream &out)
4304  {
4305  do_write_povray(patches, data_names, flags, out);
4306  }
4307 
4308 
4309 
4310  template <int dim, int spacedim>
4311  void
4313  const std::vector<Patch<dim, spacedim>> & /*patches*/,
4314  const std::vector<std::string> & /*data_names*/,
4315  const std::vector<
4316  std::tuple<unsigned int,
4317  unsigned int,
4318  std::string,
4320  const EpsFlags & /*flags*/,
4321  std::ostream & /*out*/)
4322  {
4323  // not implemented, see the documentation of the function
4324  AssertThrow(dim == 2, ExcNotImplemented());
4325  }
4326 
4327 
4328  template <int spacedim>
4329  void
4331  const std::vector<Patch<2, spacedim>> &patches,
4332  const std::vector<std::string> & /*data_names*/,
4333  const std::vector<
4334  std::tuple<unsigned int,
4335  unsigned int,
4336  std::string,
4338  const EpsFlags &flags,
4339  std::ostream &out)
4340  {
4341  AssertThrow(out.fail() == false, ExcIO());
4342 
4343 #ifndef DEAL_II_WITH_MPI
4344  // verify that there are indeed patches to be written out. most of the
4345  // times, people just forget to call build_patches when there are no
4346  // patches, so a warning is in order. that said, the assertion is disabled
4347  // if we support MPI since then it can happen that on the coarsest mesh, a
4348  // processor simply has no cells it actually owns, and in that case it is
4349  // legit if there are no patches
4350  Assert(patches.size() > 0, ExcNoPatches());
4351 #else
4352  if (patches.empty())
4353  return;
4354 #endif
4355 
4356  // set up an array of cells to be written later. this array holds the cells
4357  // of all the patches as projected to the plane perpendicular to the line of
4358  // sight.
4359  //
4360  // note that they are kept sorted by the set, where we chose the value of
4361  // the center point of the cell along the line of sight as value for sorting
4362  std::multiset<EpsCell2d> cells;
4363 
4364  // two variables in which we will store the minimum and maximum values of
4365  // the field to be used for colorization
4366  float min_color_value = std::numeric_limits<float>::max();
4367  float max_color_value = std::numeric_limits<float>::min();
4368 
4369  // Array for z-coordinates of points. The elevation determined by a function
4370  // if spacedim=2 or the z-coordinate of the grid point if spacedim=3
4371  double heights[4] = {0, 0, 0, 0};
4372 
4373  // compute the cells for output and enter them into the set above note that
4374  // since dim==2, we have exactly four vertices per patch and per cell
4375  for (const auto &patch : patches)
4376  {
4377  const unsigned int n_subdivisions = patch.n_subdivisions;
4378  const unsigned int n = n_subdivisions + 1;
4379  const unsigned int d1 = 1;
4380  const unsigned int d2 = n;
4381 
4382  for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
4383  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
4384  {
4385  Point<spacedim> points[4];
4386  points[0] =
4387  get_equispaced_location(patch, {i1, i2}, n_subdivisions);
4388  points[1] =
4389  get_equispaced_location(patch, {i1 + 1, i2}, n_subdivisions);
4390  points[2] =
4391  get_equispaced_location(patch, {i1, i2 + 1}, n_subdivisions);
4392  points[3] = get_equispaced_location(patch,
4393  {i1 + 1, i2 + 1},
4394  n_subdivisions);
4395 
4396  switch (spacedim)
4397  {
4398  case 2:
4399  Assert((flags.height_vector < patch.data.n_rows()) ||
4400  patch.data.n_rows() == 0,
4402  0,
4403  patch.data.n_rows()));
4404  heights[0] =
4405  patch.data.n_rows() != 0 ?
4406  patch.data(flags.height_vector, i1 * d1 + i2 * d2) *
4407  flags.z_scaling :
4408  0;
4409  heights[1] = patch.data.n_rows() != 0 ?
4410  patch.data(flags.height_vector,
4411  (i1 + 1) * d1 + i2 * d2) *
4412  flags.z_scaling :
4413  0;
4414  heights[2] = patch.data.n_rows() != 0 ?
4415  patch.data(flags.height_vector,
4416  i1 * d1 + (i2 + 1) * d2) *
4417  flags.z_scaling :
4418  0;
4419  heights[3] = patch.data.n_rows() != 0 ?
4420  patch.data(flags.height_vector,
4421  (i1 + 1) * d1 + (i2 + 1) * d2) *
4422  flags.z_scaling :
4423  0;
4424 
4425  break;
4426  case 3:
4427  // Copy z-coordinates into the height vector
4428  for (unsigned int i = 0; i < 4; ++i)
4429  heights[i] = points[i](2);
4430  break;
4431  default:
4432  Assert(false, ExcNotImplemented());
4433  }
4434 
4435 
4436  // now compute the projection of the bilinear cell given by the
4437  // four vertices and their heights and write them to a proper cell
4438  // object. note that we only need the first two components of the
4439  // projected position for output, but we need the value along the
4440  // line of sight for sorting the cells for back-to- front-output
4441  //
4442  // this computation was first written by Stefan Nauber. please
4443  // no-one ask me why it works that way (or may be not), especially
4444  // not about the angles and the sign of the height field, I don't
4445  // know it.
4446  EpsCell2d eps_cell;
4447  const double pi = numbers::PI;
4448  const double cx =
4449  -std::cos(pi - flags.azimut_angle * 2 * pi / 360.),
4450  cz = -std::cos(flags.turn_angle * 2 * pi / 360.),
4451  sx =
4452  std::sin(pi - flags.azimut_angle * 2 * pi / 360.),
4453  sz = std::sin(flags.turn_angle * 2 * pi / 360.);
4454  for (unsigned int vertex = 0; vertex < 4; ++vertex)
4455  {
4456  const double x = points[vertex](0), y = points[vertex](1),
4457  z = -heights[vertex];
4458 
4459  eps_cell.vertices[vertex](0) = -cz * x + sz * y;
4460  eps_cell.vertices[vertex](1) =
4461  -cx * sz * x - cx * cz * y - sx * z;
4462 
4463  // ( 1 0 0 )
4464  // D1 = ( 0 cx -sx )
4465  // ( 0 sx cx )
4466 
4467  // ( cy 0 sy )
4468  // Dy = ( 0 1 0 )
4469  // (-sy 0 cy )
4470 
4471  // ( cz -sz 0 )
4472  // Dz = ( sz cz 0 )
4473  // ( 0 0 1 )
4474 
4475  // ( cz -sz 0 )( 1 0 0 )(x) (
4476  // cz*x-sz*(cx*y-sx*z)+0*(sx*y+cx*z) )
4477  // Dxz = ( sz cz 0 )( 0 cx -sx )(y) = (
4478  // sz*x+cz*(cx*y-sx*z)+0*(sx*y+cx*z) )
4479  // ( 0 0 1 )( 0 sx cx )(z) ( 0*x+
4480  // *(cx*y-sx*z)+1*(sx*y+cx*z) )
4481  }
4482 
4483  // compute coordinates of center of cell
4484  const Point<spacedim> center_point =
4485  (points[0] + points[1] + points[2] + points[3]) / 4;
4486  const double center_height =
4487  -(heights[0] + heights[1] + heights[2] + heights[3]) / 4;
4488 
4489  // compute the depth into the picture
4490  eps_cell.depth = -sx * sz * center_point(0) -
4491  sx * cz * center_point(1) + cx * center_height;
4492 
4493  if (flags.draw_cells && flags.shade_cells)
4494  {
4495  Assert((flags.color_vector < patch.data.n_rows()) ||
4496  patch.data.n_rows() == 0,
4498  0,
4499  patch.data.n_rows()));
4500  const double color_values[4] = {
4501  patch.data.n_rows() != 0 ?
4502  patch.data(flags.color_vector, i1 * d1 + i2 * d2) :
4503  1,
4504 
4505  patch.data.n_rows() != 0 ?
4506  patch.data(flags.color_vector, (i1 + 1) * d1 + i2 * d2) :
4507  1,
4508 
4509  patch.data.n_rows() != 0 ?
4510  patch.data(flags.color_vector, i1 * d1 + (i2 + 1) * d2) :
4511  1,
4512 
4513  patch.data.n_rows() != 0 ?
4514  patch.data(flags.color_vector,
4515  (i1 + 1) * d1 + (i2 + 1) * d2) :
4516  1};
4517 
4518  // set color value to average of the value at the vertices
4519  eps_cell.color_value = (color_values[0] + color_values[1] +
4520  color_values[3] + color_values[2]) /
4521  4;
4522 
4523  // update bounds of color field
4524  min_color_value =
4525  std::min(min_color_value, eps_cell.color_value);
4526  max_color_value =
4527  std::max(max_color_value, eps_cell.color_value);
4528  }
4529 
4530  // finally add this cell
4531  cells.insert(eps_cell);
4532  }
4533  }
4534 
4535  // find out minimum and maximum x and y coordinates to compute offsets and
4536  // scaling factors
4537  double x_min = cells.begin()->vertices[0](0);
4538  double x_max = x_min;
4539  double y_min = cells.begin()->vertices[0](1);
4540  double y_max = y_min;
4541 
4542  for (const auto &cell : cells)
4543  for (const auto &vertex : cell.vertices)
4544  {
4545  x_min = std::min(x_min, vertex(0));
4546  x_max = std::max(x_max, vertex(0));
4547  y_min = std::min(y_min, vertex(1));
4548  y_max = std::max(y_max, vertex(1));
4549  }
4550 
4551  // scale in x-direction such that in the output 0 <= x <= 300. don't scale
4552  // in y-direction to preserve the shape of the triangulation
4553  const double scale =
4554  (flags.size /
4555  (flags.size_type == EpsFlags::width ? x_max - x_min : y_min - y_max));
4556 
4557  const Point<2> offset(x_min, y_min);
4558 
4559 
4560  // now write preamble
4561  {
4562  out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4563  << "%%Title: deal.II Output" << '\n'
4564  << "%%Creator: the deal.II library" << '\n'
4565  << "%%Creation Date: " << Utilities::System::get_date() << " - "
4566  << Utilities::System::get_time() << '\n'
4567  << "%%BoundingBox: "
4568  // lower left corner
4569  << "0 0 "
4570  // upper right corner
4571  << static_cast<unsigned int>((x_max - x_min) * scale + 0.5) << ' '
4572  << static_cast<unsigned int>((y_max - y_min) * scale + 0.5) << '\n';
4573 
4574  // define some abbreviations to keep the output small:
4575  // m=move turtle to
4576  // l=define a line
4577  // s=set rgb color
4578  // sg=set gray value
4579  // lx=close the line and plot the line
4580  // lf=close the line and fill the interior
4581  out << "/m {moveto} bind def" << '\n'
4582  << "/l {lineto} bind def" << '\n'
4583  << "/s {setrgbcolor} bind def" << '\n'
4584  << "/sg {setgray} bind def" << '\n'
4585  << "/lx {lineto closepath stroke} bind def" << '\n'
4586  << "/lf {lineto closepath fill} bind def" << '\n';
4587 
4588  out << "%%EndProlog" << '\n' << '\n';
4589  // set fine lines
4590  out << flags.line_width << " setlinewidth" << '\n';
4591  }
4592 
4593  // check if min and max values for the color are actually different. If
4594  // that is not the case (such things happen, for example, in the very first
4595  // time step of a time dependent problem, if the initial values are zero),
4596  // all values are equal, and then we can draw everything in an arbitrary
4597  // color. Thus, change one of the two values arbitrarily
4598  if (max_color_value == min_color_value)
4599  max_color_value = min_color_value + 1;
4600 
4601  // now we've got all the information we need. write the cells. note: due to
4602  // the ordering, we traverse the list of cells back-to-front
4603  for (const auto &cell : cells)
4604  {
4605  if (flags.draw_cells)
4606  {
4607  if (flags.shade_cells)
4608  {
4609  const EpsFlags::RgbValues rgb_values =
4610  (*flags.color_function)(cell.color_value,
4611  min_color_value,
4612  max_color_value);
4613 
4614  // write out color
4615  if (rgb_values.is_grey())
4616  out << rgb_values.red << " sg ";
4617  else
4618  out << rgb_values.red << ' ' << rgb_values.green << ' '
4619  << rgb_values.blue << " s ";
4620  }
4621  else
4622  out << "1 sg ";
4623 
4624  out << (cell.vertices[0] - offset) * scale << " m "
4625  << (cell.vertices[1] - offset) * scale << " l "
4626  << (cell.vertices[3] - offset) * scale << " l "
4627  << (cell.vertices[2] - offset) * scale << " lf" << '\n';
4628  }
4629 
4630  if (flags.draw_mesh)
4631  out << "0 sg " // draw lines in black
4632  << (cell.vertices[0] - offset) * scale << " m "
4633  << (cell.vertices[1] - offset) * scale << " l "
4634  << (cell.vertices[3] - offset) * scale << " l "
4635  << (cell.vertices[2] - offset) * scale << " lx" << '\n';
4636  }
4637  out << "showpage" << '\n';
4638 
4639  out.flush();
4640 
4641  AssertThrow(out.fail() == false, ExcIO());
4642  }
4643 
4644 
4645 
4646  template <int dim, int spacedim>
4647  void
4649  const std::vector<Patch<dim, spacedim>> &patches,
4650  const std::vector<std::string> &data_names,
4651  const std::vector<
4652  std::tuple<unsigned int,
4653  unsigned int,
4654  std::string,
4656  const GmvFlags &flags,
4657  std::ostream &out)
4658  {
4659  // The gmv format does not support cells that only consist of a single
4660  // point. It does support the output of point data using the keyword
4661  // 'tracers' instead of 'nodes' and 'cells', but this output format is
4662  // currently not implemented.
4663  AssertThrow(dim > 0, ExcNotImplemented());
4664 
4665  Assert(dim <= 3, ExcNotImplemented());
4666  AssertThrow(out.fail() == false, ExcIO());
4667 
4668 #ifndef DEAL_II_WITH_MPI
4669  // verify that there are indeed patches to be written out. most of the
4670  // times, people just forget to call build_patches when there are no
4671  // patches, so a warning is in order. that said, the assertion is disabled
4672  // if we support MPI since then it can happen that on the coarsest mesh, a
4673  // processor simply has no cells it actually owns, and in that case it is
4674  // legit if there are no patches
4675  Assert(patches.size() > 0, ExcNoPatches());
4676 #else
4677  if (patches.empty())
4678  return;
4679 #endif
4680 
4681  GmvStream gmv_out(out, flags);
4682  const unsigned int n_data_sets = data_names.size();
4683  // check against # of data sets in first patch. checks against all other
4684  // patches are made in write_gmv_reorder_data_vectors
4685  Assert((patches[0].data.n_rows() == n_data_sets &&
4686  !patches[0].points_are_available) ||
4687  (patches[0].data.n_rows() == n_data_sets + spacedim &&
4688  patches[0].points_are_available),
4689  ExcDimensionMismatch(patches[0].points_are_available ?
4690  (n_data_sets + spacedim) :
4691  n_data_sets,
4692  patches[0].data.n_rows()));
4693 
4694  //---------------------
4695  // preamble
4696  out << "gmvinput ascii" << '\n' << '\n';
4697 
4698  // first count the number of cells and cells for later use
4699  unsigned int n_nodes;
4700  unsigned int n_cells;
4701  std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
4702 
4703  // For the format we write here, we need to write all node values relating
4704  // to one variable at a time. We could in principle do this by looping
4705  // over all patches and extracting the values corresponding to the one
4706  // variable we're dealing with right now, and then start the process over
4707  // for the next variable with another loop over all patches.
4708  //
4709  // An easier way is to create a global table that for each variable
4710  // lists all values. This copying of data vectors can be done in the
4711  // background while we're already working on vertices and cells,
4712  // so do this on a separate task and when wanting to write out the
4713  // data, we wait for that task to finish.
4715  create_global_data_table_task = Threads::new_task(
4716  [&patches]() { return create_global_data_table(patches); });
4717 
4718  //-----------------------------
4719  // first make up a list of used vertices along with their coordinates
4720  //
4721  // note that we have to print 3 dimensions
4722  out << "nodes " << n_nodes << '\n';
4723  for (unsigned int d = 0; d < spacedim; ++d)
4724  {
4725  gmv_out.selected_component = d;
4726  write_nodes(patches, gmv_out);
4727  out << '\n';
4728  }
4729  gmv_out.selected_component = numbers::invalid_unsigned_int;
4730 
4731  for (unsigned int d = spacedim; d < 3; ++d)
4732  {
4733  for (unsigned int i = 0; i < n_nodes; ++i)
4734  out << "0 ";
4735  out << '\n';
4736  }
4737 
4738  //-------------------------------
4739  // now for the cells. note that vertices are counted from 1 onwards
4740  out << "cells " << n_cells << '\n';
4741  write_cells(patches, gmv_out);
4742 
4743  //-------------------------------------
4744  // data output.
4745  out << "variable" << '\n';
4746 
4747  // Wait for the reordering to be done and retrieve the reordered data:
4748  const Table<2, double> data_vectors =
4749  std::move(*create_global_data_table_task.return_value());
4750 
4751  // then write data. the '1' means: node data (as opposed to cell data, which
4752  // we do not support explicitly here)
4753  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4754  {
4755  out << data_names[data_set] << " 1" << '\n';
4756  std::copy(data_vectors[data_set].begin(),
4757  data_vectors[data_set].end(),
4758  std::ostream_iterator<double>(out, " "));
4759  out << '\n' << '\n';
4760  }
4761 
4762 
4763 
4764  // end of variable section
4765  out << "endvars" << '\n';
4766 
4767  // end of output
4768  out << "endgmv" << '\n';
4769 
4770  // make sure everything now gets to disk
4771  out.flush();
4772 
4773  // assert the stream is still ok
4774  AssertThrow(out.fail() == false, ExcIO());
4775  }
4776 
4777 
4778 
4779  template <int dim, int spacedim>
4780  void
4782  const std::vector<Patch<dim, spacedim>> &patches,
4783  const std::vector<std::string> &data_names,
4784  const std::vector<
4785  std::tuple<unsigned int,
4786  unsigned int,
4787  std::string,
4789  const TecplotFlags &flags,
4790  std::ostream &out)
4791  {
4792  AssertThrow(out.fail() == false, ExcIO());
4793 
4794  // The FEBLOCK or FEPOINT formats of tecplot only allows full elements (e.g.
4795  // triangles), not single points. Other tecplot format allow point output,
4796  // but they are currently not implemented.
4797  AssertThrow(dim > 0, ExcNotImplemented());
4798 
4799 #ifndef DEAL_II_WITH_MPI
4800  // verify that there are indeed patches to be written out. most of the
4801  // times, people just forget to call build_patches when there are no
4802  // patches, so a warning is in order. that said, the assertion is disabled
4803  // if we support MPI since then it can happen that on the coarsest mesh, a
4804  // processor simply has no cells it actually owns, and in that case it is
4805  // legit if there are no patches
4806  Assert(patches.size() > 0, ExcNoPatches());
4807 #else
4808  if (patches.empty())
4809  return;
4810 #endif
4811 
4812  TecplotStream tecplot_out(out, flags);
4813 
4814  const unsigned int n_data_sets = data_names.size();
4815  // check against # of data sets in first patch. checks against all other
4816  // patches are made in write_gmv_reorder_data_vectors
4817  Assert((patches[0].data.n_rows() == n_data_sets &&
4818  !patches[0].points_are_available) ||
4819  (patches[0].data.n_rows() == n_data_sets + spacedim &&
4820  patches[0].points_are_available),
4821  ExcDimensionMismatch(patches[0].points_are_available ?
4822  (n_data_sets + spacedim) :
4823  n_data_sets,
4824  patches[0].data.n_rows()));
4825 
4826  // first count the number of cells and cells for later use
4827  unsigned int n_nodes;
4828  unsigned int n_cells;
4829  std::tie(n_nodes, n_cells) = count_nodes_and_cells(patches);
4830 
4831  //---------
4832  // preamble
4833  {
4834  out
4835  << "# This file was generated by the deal.II library." << '\n'
4836  << "# Date = " << Utilities::System::get_date() << '\n'
4837  << "# Time = " << Utilities::System::get_time() << '\n'
4838  << "#" << '\n'
4839  << "# For a description of the Tecplot format see the Tecplot documentation."
4840  << '\n'
4841  << "#" << '\n';
4842 
4843 
4844  out << "Variables=";
4845 
4846  switch (spacedim)
4847  {
4848  case 1:
4849  out << "\"x\"";
4850  break;
4851  case 2:
4852  out << "\"x\", \"y\"";
4853  break;
4854  case 3:
4855  out << "\"x\", \"y\", \"z\"";
4856  break;
4857  default:
4858  Assert(false, ExcNotImplemented());
4859  }
4860 
4861  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4862  out << ", \"" << data_names[data_set] << "\"";
4863 
4864  out << '\n';
4865 
4866  out << "zone ";
4867  if (flags.zone_name)
4868  out << "t=\"" << flags.zone_name << "\" ";
4869 
4870  if (flags.solution_time >= 0.0)
4871  out << "strandid=1, solutiontime=" << flags.solution_time << ", ";
4872 
4873  out << "f=feblock, n=" << n_nodes << ", e=" << n_cells
4874  << ", et=" << tecplot_cell_type[dim] << '\n';
4875  }
4876 
4877 
4878  // For the format we write here, we need to write all node values relating
4879  // to one variable at a time. We could in principle do this by looping
4880  // over all patches and extracting the values corresponding to the one
4881  // variable we're dealing with right now, and then start the process over
4882  // for the next variable with another loop over all patches.
4883  //
4884  // An easier way is to create a global table that for each variable
4885  // lists all values. This copying of data vectors can be done in the
4886  // background while we're already working on vertices and cells,
4887  // so do this on a separate task and when wanting to write out the
4888  // data, we wait for that task to finish.
4890  create_global_data_table_task = Threads::new_task(
4891  [&patches]() { return create_global_data_table(patches); });
4892 
4893  //-----------------------------
4894  // first make up a list of used vertices along with their coordinates
4895 
4896 
4897  for (unsigned int d = 0; d < spacedim; ++d)
4898  {
4899  tecplot_out.selected_component = d;
4900  write_nodes(patches, tecplot_out);
4901  out << '\n';
4902  }
4903 
4904 
4905  //-------------------------------------
4906  // data output.
4907  //
4908  // Wait for the reordering to be done and retrieve the reordered data:
4909  const Table<2, double> data_vectors =
4910  std::move(*create_global_data_table_task.return_value());
4911 
4912  // then write data.
4913  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
4914  {
4915  std::copy(data_vectors[data_set].begin(),
4916  data_vectors[data_set].end(),
4917  std::ostream_iterator<double>(out, "\n"));
4918  out << '\n';
4919  }
4920 
4921  write_cells(patches, tecplot_out);
4922 
4923  // make sure everything now gets to disk
4924  out.flush();
4925 
4926  // assert the stream is still ok
4927  AssertThrow(out.fail() == false, ExcIO());
4928  }
4929 
4930 
4931 
4932  template <int dim, int spacedim>
4933  void
4935  const std::vector<Patch<dim, spacedim>> &patches,
4936  const std::vector<std::string> &data_names,
4937  const std::vector<
4938  std::tuple<unsigned int,
4939  unsigned int,
4940  std::string,
4942  &nonscalar_data_ranges,
4943  const VtkFlags &flags,
4944  std::ostream &out)
4945  {
4946  AssertThrow(out.fail() == false, ExcIO());
4947 
4948 #ifndef DEAL_II_WITH_MPI
4949  // verify that there are indeed patches to be written out. most of the
4950  // times, people just forget to call build_patches when there are no
4951  // patches, so a warning is in order. that said, the assertion is disabled
4952  // if we support MPI since then it can happen that on the coarsest mesh, a
4953  // processor simply has no cells it actually owns, and in that case it is
4954  // legit if there are no patches
4955  Assert(patches.size() > 0, ExcNoPatches());
4956 #else
4957  if (patches.empty())
4958  return;
4959 #endif
4960 
4961  VtkStream vtk_out(out, flags);
4962 
4963  const unsigned int n_data_sets = data_names.size();
4964  // check against # of data sets in first patch.
4965  if (patches[0].points_are_available)
4966  {
4967  AssertDimension(n_data_sets + spacedim, patches[0].data.n_rows());
4968  }
4969  else
4970  {
4971  AssertDimension(n_data_sets, patches[0].data.n_rows());
4972  }
4973 
4974  //---------------------
4975  // preamble
4976  {
4977  out << "# vtk DataFile Version 3.0" << '\n'
4978  << "#This file was generated by the deal.II library";
4979  if (flags.print_date_and_time)
4980  {
4981  out << " on " << Utilities::System::get_date() << " at "
4983  }
4984  else
4985  out << '.';
4986  out << '\n' << "ASCII" << '\n';
4987  // now output the data header
4988  out << "DATASET UNSTRUCTURED_GRID\n" << '\n';
4989  }
4990 
4991  // if desired, output time and cycle of the simulation, following the
4992  // instructions at
4993  // http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
4994  {
4995  const unsigned int n_metadata =
4996  ((flags.cycle != std::numeric_limits<unsigned int>::min() ? 1 : 0) +
4997  (flags.time != std::numeric_limits<double>::min() ? 1 : 0));
4998  if (n_metadata > 0)
4999  {
5000  out << "FIELD FieldData " << n_metadata << '\n';
5001 
5003  {
5004  out << "CYCLE 1 1 int\n" << flags.cycle << '\n';
5005  }
5006  if (flags.time != std::numeric_limits<double>::min())
5007  {
5008  out << "TIME 1 1 double\n" << flags.time << '\n';
5009  }
5010  }
5011  }
5012 
5013  // first count the number of cells and cells for later use
5014  unsigned int n_nodes;
5015  unsigned int n_cells;
5016  unsigned int n_points_and_n_cells;
5017  std::tie(n_nodes, n_cells, n_points_and_n_cells) =
5018  count_nodes_and_cells_and_points(patches, flags.write_higher_order_cells);
5019 
5020  // For the format we write here, we need to write all node values relating
5021  // to one variable at a time. We could in principle do this by looping
5022  // over all patches and extracting the values corresponding to the one
5023  // variable we're dealing with right now, and then start the process over
5024  // for the next variable with another loop over all patches.
5025  //
5026  // An easier way is to create a global table that for each variable
5027  // lists all values. This copying of data vectors can be done in the
5028  // background while we're already working on vertices and cells,
5029  // so do this on a separate task and when wanting to write out the
5030  // data, we wait for that task to finish.
5032  create_global_data_table_task = Threads::new_task(
5033  [&patches]() { return create_global_data_table(patches); });
5034 
5035  //-----------------------------
5036  // first make up a list of used vertices along with their coordinates
5037  //
5038  // note that we have to print d=1..3 dimensions
5039  out << "POINTS " << n_nodes << " double" << '\n';
5040  write_nodes(patches, vtk_out);
5041  out << '\n';
5042  //-------------------------------
5043  // now for the cells
5044  out << "CELLS " << n_cells << ' ' << n_points_and_n_cells << '\n';
5045  if (flags.write_higher_order_cells)
5046  write_high_order_cells(patches, vtk_out, /* legacy_format = */ true);
5047  else
5048  write_cells(patches, vtk_out);
5049  out << '\n';
5050  // next output the types of the cells. since all cells are the same, this is
5051  // simple
5052  out << "CELL_TYPES " << n_cells << '\n';
5053 
5054  // need to distinguish between linear cells, simplex cells (linear or
5055  // quadratic), and high order cells
5056  for (const auto &patch : patches)
5057  {
5058  const auto vtk_cell_id =
5059  extract_vtk_patch_info(patch, flags.write_higher_order_cells);
5060 
5061  for (unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5062  out << ' ' << vtk_cell_id[0];
5063  }
5064 
5065  out << '\n';
5066  //-------------------------------------
5067  // data output.
5068 
5069  // Wait for the reordering to be done and retrieve the reordered data:
5070  const Table<2, double> data_vectors =
5071  std::move(*create_global_data_table_task.return_value());
5072 
5073  // then write data. the 'POINT_DATA' means: node data (as opposed to cell
5074  // data, which we do not support explicitly here). all following data sets
5075  // are point data
5076  out << "POINT_DATA " << n_nodes << '\n';
5077 
5078  // when writing, first write out all vector data, then handle the scalar
5079  // data sets that have been left over
5080  std::vector<bool> data_set_written(n_data_sets, false);
5081  for (const auto &nonscalar_data_range : nonscalar_data_ranges)
5082  {
5083  AssertThrow(std::get<3>(nonscalar_data_range) !=
5085  ExcMessage(
5086  "The VTK writer does not currently support outputting "
5087  "tensor data. Use the VTU writer instead."));
5088 
5089  AssertThrow(std::get<1>(nonscalar_data_range) >=
5090  std::get<0>(nonscalar_data_range),
5091  ExcLowerRange(std::get<1>(nonscalar_data_range),
5092  std::get<0>(nonscalar_data_range)));
5093  AssertThrow(std::get<1>(nonscalar_data_range) < n_data_sets,
5094  ExcIndexRange(std::get<1>(nonscalar_data_range),
5095  0,
5096  n_data_sets));
5097  AssertThrow(std::get<1>(nonscalar_data_range) + 1 -
5098  std::get<0>(nonscalar_data_range) <=
5099  3,
5100  ExcMessage(
5101  "Can't declare a vector with more than 3 components "
5102  "in VTK"));
5103 
5104  // mark these components as already written:
5105  for (unsigned int i = std::get<0>(nonscalar_data_range);
5106  i <= std::get<1>(nonscalar_data_range);
5107  ++i)
5108  data_set_written[i] = true;
5109 
5110  // write the header. concatenate all the component names with double
5111  // underscores unless a vector name has been specified
5112  out << "VECTORS ";
5113 
5114  if (!std::get<2>(nonscalar_data_range).empty())
5115  out << std::get<2>(nonscalar_data_range);
5116  else
5117  {
5118  for (unsigned int i = std::get<0>(nonscalar_data_range);
5119  i < std::get<1>(nonscalar_data_range);
5120  ++i)
5121  out << data_names[i] << "__";
5122  out << data_names[std::get<1>(nonscalar_data_range)];
5123  }
5124 
5125  out << " double" << '\n';
5126 
5127  // now write data. pad all vectors to have three components
5128  for (unsigned int n = 0; n < n_nodes; ++n)
5129  {
5130  switch (std::get<1>(nonscalar_data_range) -
5131  std::get<0>(nonscalar_data_range))
5132  {
5133  case 0:
5134  out << data_vectors(std::get<0>(nonscalar_data_range), n)
5135  << " 0 0" << '\n';
5136  break;
5137 
5138  case 1:
5139  out << data_vectors(std::get<0>(nonscalar_data_range), n)
5140  << ' '
5141  << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5142  << " 0" << '\n';
5143  break;
5144  case 2:
5145  out << data_vectors(std::get<0>(nonscalar_data_range), n)
5146  << ' '
5147  << data_vectors(std::get<0>(nonscalar_data_range) + 1, n)
5148  << ' '
5149  << data_vectors(std::get<0>(nonscalar_data_range) + 2, n)
5150  << '\n';
5151  break;
5152 
5153  default:
5154  // VTK doesn't support anything else than vectors with 1, 2,
5155  // or 3 components
5156  Assert(false, ExcInternalError());
5157  }
5158  }
5159  }
5160 
5161  // now do the left over scalar data sets
5162  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
5163  if (data_set_written[data_set] == false)
5164  {
5165  out << "SCALARS " << data_names[data_set] << " double 1" << '\n'
5166  << "LOOKUP_TABLE default" << '\n';
5167  std::copy(data_vectors[data_set].begin(),
5168  data_vectors[data_set].end(),
5169  std::ostream_iterator<double>(out, " "));
5170  out << '\n';
5171  }
5172 
5173  // make sure everything now gets to disk
5174  out.flush();
5175 
5176  // assert the stream is still ok
5177  AssertThrow(out.fail() == false, ExcIO());
5178  }
5179 
5180 
5181  void
5182  write_vtu_header(std::ostream &out, const VtkFlags &flags)
5183  {
5184  AssertThrow(out.fail() == false, ExcIO());
5185  out << "<?xml version=\"1.0\" ?> \n";
5186  out << "<!-- \n";
5187  out << "# vtk DataFile Version 3.0" << '\n'
5188  << "#This file was generated by the deal.II library";
5189  if (flags.print_date_and_time)
5190  {
5191  out << " on " << Utilities::System::get_time() << " at "
5193  }
5194  else
5195  out << '.';
5196  out << "\n-->\n";
5197 
5198  if (flags.write_higher_order_cells)
5199  out << "<VTKFile type=\"UnstructuredGrid\" version=\"2.2\"";
5200  else
5201  out << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\"";
5202  if (deal_ii_with_zlib &&
5204  out << " compressor=\"vtkZLibDataCompressor\"";
5205 #ifdef DEAL_II_WORDS_BIGENDIAN
5206  out << " byte_order=\"BigEndian\"";
5207 #else
5208  out << " byte_order=\"LittleEndian\"";
5209 #endif
5210  out << ">";
5211  out << '\n';
5212  out << "<UnstructuredGrid>";
5213  out << '\n';
5214  }
5215 
5216 
5217 
5218  void
5219  write_vtu_footer(std::ostream &out)
5220  {
5221  AssertThrow(out.fail() == false, ExcIO());
5222  out << " </UnstructuredGrid>\n";
5223  out << "</VTKFile>\n";
5224  }
5225 
5226 
5227 
5228  template <int dim, int spacedim>
5229  void
5231  const std::vector<Patch<dim, spacedim>> &patches,
5232  const std::vector<std::string> &data_names,
5233  const std::vector<
5234  std::tuple<unsigned int,
5235  unsigned int,
5236  std::string,
5238  &nonscalar_data_ranges,
5239  const VtkFlags &flags,
5240  std::ostream &out)
5241  {
5242  write_vtu_header(out, flags);
5243  write_vtu_main(patches, data_names, nonscalar_data_ranges, flags, out);
5244  write_vtu_footer(out);
5245 
5246  out << std::flush;
5247  }
5248 
5249 
5250  template <int dim, int spacedim>
5251  void
5253  const std::vector<Patch<dim, spacedim>> &patches,
5254  const std::vector<std::string> &data_names,
5255  const std::vector<
5256  std::tuple<unsigned int,
5257  unsigned int,
5258  std::string,
5260  &nonscalar_data_ranges,
5261  const VtkFlags &flags,
5262  std::ostream &out)
5263  {
5264  AssertThrow(out.fail() == false, ExcIO());
5265 
5266  // If the user provided physical units, make sure that they don't contain
5267  // quote characters as this would make the VTU file invalid XML and
5268  // probably lead to all sorts of difficult error messages. Other than that,
5269  // trust the user that whatever they provide makes sense somehow.
5270  for (const auto &unit : flags.physical_units)
5271  {
5272  (void)unit;
5273  Assert(
5274  unit.second.find('\"') == std::string::npos,
5275  ExcMessage(
5276  "A physical unit you provided, <" + unit.second +
5277  ">, contained a quotation mark character. This is not allowed."));
5278  }
5279 
5280 #ifndef DEAL_II_WITH_MPI
5281  // verify that there are indeed patches to be written out. most of the
5282  // times, people just forget to call build_patches when there are no
5283  // patches, so a warning is in order. that said, the assertion is disabled
5284  // if we support MPI since then it can happen that on the coarsest mesh, a
5285  // processor simply has no cells it actually owns, and in that case it is
5286  // legit if there are no patches
5287  Assert(patches.size() > 0, ExcNoPatches());
5288 #else
5289  if (patches.empty())
5290  {
5291  // we still need to output a valid vtu file, because other CPUs might
5292  // output data. This is the minimal file that is accepted by paraview
5293  // and visit. if we remove the field definitions, visit is complaining.
5294  out << "<Piece NumberOfPoints=\"0\" NumberOfCells=\"0\" >\n"
5295  << "<Cells>\n"
5296  << "<DataArray type=\"UInt8\" Name=\"types\"></DataArray>\n"
5297  << "</Cells>\n"
5298  << " <PointData Scalars=\"scalars\">\n";
5299  std::vector<bool> data_set_written(data_names.size(), false);
5300  for (const auto &nonscalar_data_range : nonscalar_data_ranges)
5301  {
5302  // mark these components as already written:
5303  for (unsigned int i = std::get<0>(nonscalar_data_range);
5304  i <= std::get<1>(nonscalar_data_range);
5305  ++i)
5306  data_set_written[i] = true;
5307 
5308  // write the header. concatenate all the component names with double
5309  // underscores unless a vector name has been specified
5310  out << " <DataArray type=\"Float32\" Name=\"";
5311 
5312  if (!std::get<2>(nonscalar_data_range).empty())
5313  out << std::get<2>(nonscalar_data_range);
5314  else
5315  {
5316  for (unsigned int i = std::get<0>(nonscalar_data_range);
5317  i < std::get<1>(nonscalar_data_range);
5318  ++i)
5319  out << data_names[i] << "__";
5320  out << data_names[std::get<1>(nonscalar_data_range)];
5321  }
5322 
5323  out << "\" NumberOfComponents=\"3\"></DataArray>\n";
5324  }
5325 
5326  for (unsigned int data_set = 0; data_set < data_names.size();
5327  ++data_set)
5328  if (data_set_written[data_set] == false)
5329  {
5330  out << " <DataArray type=\"Float32\" Name=\""
5331  << data_names[data_set] << "\"></DataArray>\n";
5332  }
5333 
5334  out << " </PointData>\n";
5335  out << "</Piece>\n";
5336 
5337  out << std::flush;
5338 
5339  return;
5340  }
5341 #endif
5342 
5343  // first up: metadata
5344  //
5345  // if desired, output time and cycle of the simulation, following the
5346  // instructions at
5347  // http://www.visitusers.org/index.php?title=Time_and_Cycle_in_VTK_files
5348  {
5349  const unsigned int n_metadata =
5350  ((flags.cycle != std::numeric_limits<unsigned int>::min() ? 1 : 0) +
5351  (flags.time != std::numeric_limits<double>::min() ? 1 : 0));
5352  if (n_metadata > 0)
5353  out << "<FieldData>\n";
5354 
5356  {
5357  out
5358  << "<DataArray type=\"Float32\" Name=\"CYCLE\" NumberOfTuples=\"1\" format=\"ascii\">"
5359  << flags.cycle << "</DataArray>\n";
5360  }
5361  if (flags.time != std::numeric_limits<double>::min())
5362  {
5363  out
5364  << "<DataArray type=\"Float32\" Name=\"TIME\" NumberOfTuples=\"1\" format=\"ascii\">"
5365  << flags.time << "</DataArray>\n";
5366  }
5367 
5368  if (n_metadata > 0)
5369  out << "</FieldData>\n";
5370  }
5371 
5372 
5373  const unsigned int n_data_sets = data_names.size();
5374  // check against # of data sets in first patch. checks against all other
5375  // patches are made in write_gmv_reorder_data_vectors
5376  if (patches[0].points_are_available)
5377  {
5378  AssertDimension(n_data_sets + spacedim, patches[0].data.n_rows());
5379  }
5380  else
5381  {
5382  AssertDimension(n_data_sets, patches[0].data.n_rows());
5383  }
5384 
5385  const char *ascii_or_binary =
5386  (deal_ii_with_zlib &&
5388  "binary" :
5389  "ascii";
5390 
5391 
5392  // first count the number of cells and cells for later use
5393  unsigned int n_nodes;
5394  unsigned int n_cells;
5395  std::tie(n_nodes, n_cells, std::ignore) =
5396  count_nodes_and_cells_and_points(patches, flags.write_higher_order_cells);
5397 
5398  // -----------------
5399  // In the following, let us first set up a number of lambda functions that
5400  // will be used in building the different parts of the VTU file. We will
5401  // later call them in turn on different tasks.
5402  // first make up a list of used vertices along with their coordinates
5403  const auto stringize_vertex_information = [&patches,
5404  &flags,
5405  output_precision =
5406  out.precision(),
5407  ascii_or_binary]() {
5408  std::ostringstream o;
5409  o << " <Points>\n";
5410  o << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\""
5411  << ascii_or_binary << "\">\n";
5412  const std::vector<Point<spacedim>> node_positions =
5413  get_node_positions(patches);
5414 
5415  // VTK/VTU always wants to see three coordinates, even if we are
5416  // in 1d or 2d. So pad node positions with zeros as appropriate.
5417  std::vector<float> node_coordinates_3d;
5418  node_coordinates_3d.reserve(node_positions.size() * 3);
5419  for (const auto &node_position : node_positions)
5420  {
5421  for (unsigned int d = 0; d < 3; ++d)
5422  if (d < spacedim)
5423  node_coordinates_3d.emplace_back(node_position[d]);
5424  else
5425  node_coordinates_3d.emplace_back(0.0f);
5426  }
5427  o << vtu_stringize_array(node_coordinates_3d,
5428  flags.compression_level,
5429  output_precision)
5430  << '\n';
5431  o << " </DataArray>\n";
5432  o << " </Points>\n\n";
5433 
5434  return o.str();
5435  };
5436 
5437 
5438  //-------------------------------
5439  // Now for the cells. The first part of this is how vertices
5440  // build cells.
5441  const auto stringize_cell_to_vertex_information = [&patches,
5442  &flags,
5443  ascii_or_binary,
5444  output_precision =
5445  out.precision()]() {
5446  std::ostringstream o;
5447 
5448  o << " <Cells>\n";
5449  o << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\""
5450  << ascii_or_binary << "\">\n";
5451 
5452  std::vector<int32_t> cells;
5453  Assert(dim <= 3, ExcNotImplemented());
5454 
5455  unsigned int first_vertex_of_patch = 0;
5456 
5457  for (const auto &patch : patches)
5458  {
5459  // First treat a slight oddball case: For triangles and tetrahedra,
5460  // the case with n_subdivisions==2 is treated as if the cell was
5461  // output as a single, quadratic, cell rather than as one would
5462  // expect as 4 sub-cells (for triangles; and the corresponding
5463  // number of sub-cells for tetrahedra). This is courtesy of some
5464  // special-casing in the function extract_vtk_patch_info().
5465  if ((dim >= 2) &&
5466  (patch.reference_cell == ReferenceCells::get_simplex<dim>()) &&
5467  (patch.n_subdivisions == 2))
5468  {
5469  const unsigned int n_points = patch.data.n_cols();
5470  Assert((dim == 2 && n_points == 6) ||
5471  (dim == 3 && n_points == 10),
5472  ExcInternalError());
5473 
5474  if (deal_ii_with_zlib &&
5475  (flags.compression_level !=
5477  {
5478  for (unsigned int i = 0; i < n_points; ++i)
5479  cells.push_back(first_vertex_of_patch + i);
5480  }
5481  else
5482  {
5483  for (unsigned int i = 0; i < n_points; ++i)
5484  o << '\t' << first_vertex_of_patch + i;
5485  o << '\n';
5486  }
5487 
5488  first_vertex_of_patch += n_points;
5489  }
5490  // Then treat all of the other non-hypercube cases since they can
5491  // currently not be subdivided (into sub-cells, or into higher-order
5492  // cells):
5493  else if (patch.reference_cell != ReferenceCells::get_hypercube<dim>())
5494  {
5495  Assert(patch.n_subdivisions == 1, ExcNotImplemented());
5496 
5497  const unsigned int n_points = patch.data.n_cols();
5498 
5499  if (deal_ii_with_zlib &&
5500  (flags.compression_level !=
5502  {
5503  for (unsigned int i = 0; i < n_points; ++i)
5504  cells.push_back(
5505  first_vertex_of_patch +
5507  }
5508  else
5509  {
5510  for (unsigned int i = 0; i < n_points; ++i)
5511  o << '\t'
5512  << (first_vertex_of_patch +
5514  o << '\n';
5515  }
5516 
5517  first_vertex_of_patch += n_points;
5518  }
5519  else // a hypercube cell
5520  {
5521  const unsigned int n_subdivisions = patch.n_subdivisions;
5522  const unsigned int n_points_per_direction = n_subdivisions + 1;
5523 
5524  std::vector<unsigned> local_vertex_order;
5525 
5526  // Output the current state of the local_vertex_order array,
5527  // then clear it:
5528  const auto flush_current_cell = [&flags,
5529  &o,
5530  &cells,
5531  first_vertex_of_patch,
5532  &local_vertex_order]() {
5533  if (deal_ii_with_zlib &&
5534  (flags.compression_level !=
5536  {
5537  for (const auto &c : local_vertex_order)
5538  cells.push_back(first_vertex_of_patch + c);
5539  }
5540  else
5541  {
5542  for (const auto &c : local_vertex_order)
5543  o << '\t' << first_vertex_of_patch + c;
5544  o << '\n';
5545  }
5546 
5547  local_vertex_order.clear();
5548  };
5549 
5550  if (flags.write_higher_order_cells == false)
5551  {
5552  local_vertex_order.reserve(Utilities::fixed_power<dim>(2));
5553 
5554  switch (dim)
5555  {
5556  case 0:
5557  {
5558  local_vertex_order.emplace_back(0);
5559  flush_current_cell();
5560  break;
5561  }
5562 
5563  case 1:
5564  {
5565  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5566  {
5567  const unsigned int starting_offset = i1;
5568  local_vertex_order.emplace_back(starting_offset);
5569  local_vertex_order.emplace_back(starting_offset +
5570  1);
5571  flush_current_cell();
5572  }
5573  break;
5574  }
5575 
5576  case 2:
5577  {
5578  for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5579  for (unsigned int i1 = 0; i1 < n_subdivisions; ++i1)
5580  {
5581  const unsigned int starting_offset =
5582  i2 * n_points_per_direction + i1;
5583  local_vertex_order.emplace_back(
5584  starting_offset);
5585  local_vertex_order.emplace_back(
5586  starting_offset + 1);
5587  local_vertex_order.emplace_back(
5588  starting_offset + n_points_per_direction + 1);
5589  local_vertex_order.emplace_back(
5590  starting_offset + n_points_per_direction);
5591  flush_current_cell();
5592  }
5593  break;
5594  }
5595 
5596  case 3:
5597  {
5598  for (unsigned int i3 = 0; i3 < n_subdivisions; ++i3)
5599  for (unsigned int i2 = 0; i2 < n_subdivisions; ++i2)
5600  for (unsigned int i1 = 0; i1 < n_subdivisions;
5601  ++i1)
5602  {
5603  const unsigned int starting_offset =
5604  i3 * n_points_per_direction *
5605  n_points_per_direction +
5606  i2 * n_points_per_direction + i1;
5607  local_vertex_order.emplace_back(
5608  starting_offset);
5609  local_vertex_order.emplace_back(
5610  starting_offset + 1);
5611  local_vertex_order.emplace_back(
5612  starting_offset + n_points_per_direction +
5613  1);
5614  local_vertex_order.emplace_back(
5615  starting_offset + n_points_per_direction);
5616  local_vertex_order.emplace_back(
5617  starting_offset + n_points_per_direction *
5618  n_points_per_direction);
5619  local_vertex_order.emplace_back(
5620  starting_offset +
5621  n_points_per_direction *
5622  n_points_per_direction +
5623  1);
5624  local_vertex_order.emplace_back(
5625  starting_offset +
5626  n_points_per_direction *
5627  n_points_per_direction +
5628  n_points_per_direction + 1);
5629  local_vertex_order.emplace_back(
5630  starting_offset +
5631  n_points_per_direction *
5632  n_points_per_direction +
5633  n_points_per_direction);
5634  flush_current_cell();
5635  }
5636  break;
5637  }
5638 
5639  default:
5640  Assert(false, ExcNotImplemented());
5641  }
5642  }
5643  else // use higher-order output
5644  {
5645  local_vertex_order.resize(
5646  Utilities::fixed_power<dim>(n_points_per_direction));
5647 
5648  switch (dim)
5649  {
5650  case 0:
5651  {
5652  Assert(false,
5653  ExcMessage(
5654  "Point-like cells should not be possible "
5655  "when writing higher-order cells."));
5656  break;
5657  }
5658  case 1:
5659  {
5660  for (unsigned int i1 = 0; i1 < n_subdivisions + 1;
5661  ++i1)
5662  {
5663  const unsigned int local_index = i1;
5664  const unsigned int connectivity_index =
5665  patch.reference_cell
5666  .template vtk_lexicographic_to_node_index<1>(
5667  {{i1}},
5668  {{n_subdivisions}},
5669  /* use VTU, not VTK: */ false);
5670  local_vertex_order[connectivity_index] =
5671  local_index;
5672  flush_current_cell();
5673  }
5674 
5675  break;
5676  }
5677  case 2:
5678  {
5679  for (unsigned int i2 = 0; i2 < n_subdivisions + 1;
5680  ++i2)
5681  for (unsigned int i1 = 0; i1 < n_subdivisions + 1;
5682  ++i1)
5683  {
5684  const unsigned int local_index =
5685  i2 * n_points_per_direction + i1;
5686  const unsigned int connectivity_index =
5687  patch.reference_cell
5688  .template vtk_lexicographic_to_node_index<
5689  2>({{i1, i2}},
5690  {{n_subdivisions, n_subdivisions}},
5691  /* use VTU, not VTK: */ false);
5692  local_vertex_order[connectivity_index] =
5693  local_index;
5694  }
5695  flush_current_cell();
5696 
5697  break;
5698  }
5699  case 3:
5700  {
5701  for (unsigned int i3 = 0; i3 < n_subdivisions + 1;
5702  ++i3)
5703  for (unsigned int i2 = 0; i2 < n_subdivisions + 1;
5704  ++i2)
5705  for (unsigned int i1 = 0; i1 < n_subdivisions + 1;
5706  ++i1)
5707  {
5708  const unsigned int local_index =
5709  i3 * n_points_per_direction *
5710  n_points_per_direction +
5711  i2 * n_points_per_direction + i1;
5712  const unsigned int connectivity_index =
5713  patch.reference_cell
5714  .template vtk_lexicographic_to_node_index<
5715  3>({{i1, i2, i3}},
5716  {{n_subdivisions,
5717  n_subdivisions,
5718  n_subdivisions}},
5719  /* use VTU, not VTK: */ false);
5720  local_vertex_order[connectivity_index] =
5721  local_index;
5722  }
5723 
5724  flush_current_cell();
5725  break;
5726  }
5727  default:
5728  Assert(false, ExcNotImplemented());
5729  }
5730  }
5731 
5732  // Finally update the number of the first vertex of this
5733  // patch
5734  first_vertex_of_patch +=
5735  Utilities::fixed_power<dim>(patch.n_subdivisions + 1);
5736  }
5737  }
5738 
5739  // Flush the 'cells' object we created herein.
5740  if (deal_ii_with_zlib && (flags.compression_level !=
5742  {
5743  o << vtu_stringize_array(cells,
5744  flags.compression_level,
5745  output_precision)
5746  << '\n';
5747  }
5748  o << " </DataArray>\n";
5749 
5750  return o.str();
5751  };
5752 
5753 
5754  //-------------------------------
5755  // The second part of cell information is the offsets in
5756  // the array built by the previous lambda function that indicate
5757  // individual cells.
5758  //
5759  // Note that this separates XML VTU format from the VTK format; the latter
5760  // puts the number of nodes per cell in front of the connectivity list for
5761  // each cell, whereas the VTU format uses one large list of vertex indices
5762  // and a separate array of offsets.
5763  //
5764  // The third piece to cell information is that we need to
5765  // output the types of the cells.
5766  //
5767  // The following function does both of these pieces.
5768  const auto stringize_cell_offset_and_type_information =
5769  [&patches,
5770  &flags,
5771  ascii_or_binary,
5772  n_cells,
5773  output_precision = out.precision()]() {
5774  std::ostringstream o;
5775 
5776  o << " <DataArray type=\"Int32\" Name=\"offsets\" format=\""
5777  << ascii_or_binary << "\">\n";
5778 
5779  std::vector<int32_t> offsets;
5780  offsets.reserve(n_cells);
5781 
5782  // std::uint8_t might be an alias to unsigned char which is then not
5783  // printed as ascii integers
5784  std::vector<unsigned int> cell_types;
5785  cell_types.reserve(n_cells);
5786 
5787  unsigned int first_vertex_of_patch = 0;
5788 
5789  for (const auto &patch : patches)
5790  {
5791  const auto vtk_cell_id =
5792  extract_vtk_patch_info(patch, flags.write_higher_order_cells);
5793 
5794  for (unsigned int i = 0; i < vtk_cell_id[1]; ++i)
5795  {
5796  cell_types.push_back(vtk_cell_id[0]);
5797  first_vertex_of_patch += vtk_cell_id[2];
5798  offsets.push_back(first_vertex_of_patch);
5799  }
5800  }
5801 
5802  o << vtu_stringize_array(offsets,
5803  flags.compression_level,
5804  output_precision);
5805  o << '\n';
5806  o << " </DataArray>\n";
5807 
5808  o << " <DataArray type=\"UInt8\" Name=\"types\" format=\""
5809  << ascii_or_binary << "\">\n";
5810 
5811  if (deal_ii_with_zlib &&
5813  {
5814  std::vector<uint8_t> cell_types_uint8_t(cell_types.size());
5815  for (unsigned int i = 0; i < cell_types.size(); ++i)
5816  cell_types_uint8_t[i] = static_cast<std::uint8_t>(cell_types[i]);
5817 
5818  o << vtu_stringize_array(cell_types_uint8_t,
5819  flags.compression_level,
5820  output_precision);
5821  }
5822  else
5823  {
5824  o << vtu_stringize_array(cell_types,
5825  flags.compression_level,
5826  output_precision);
5827  }
5828 
5829  o << '\n';
5830  o << " </DataArray>\n";
5831  o << " </Cells>\n";
5832 
5833  return o.str();
5834  };
5835 
5836 
5837  //-------------------------------------
5838  // data output.
5839 
5840  const auto stringize_nonscalar_data_range =
5841  [&flags,
5842  &data_names,
5843  ascii_or_binary,
5844  n_data_sets,
5845  n_nodes,
5846  output_precision = out.precision()](const Table<2, float> &data_vectors,
5847  const auto &range) {
5848  std::ostringstream o;
5849 
5850  const auto first_component = std::get<0>(range);
5851  const auto last_component = std::get<1>(range);
5852  const auto &name = std::get<2>(range);
5853  const bool is_tensor =
5854  (std::get<3>(range) ==
5856  const unsigned int n_components = (is_tensor ? 9 : 3);
5857  AssertThrow(last_component >= first_component,
5858  ExcLowerRange(last_component, first_component));
5859  AssertThrow(last_component < n_data_sets,
5860  ExcIndexRange(last_component, 0, n_data_sets));
5861  if (is_tensor)
5862  {
5863  AssertThrow((last_component + 1 - first_component <= 9),
5864  ExcMessage(
5865  "Can't declare a tensor with more than 9 components "
5866  "in VTK/VTU format."));
5867  }
5868  else
5869  {
5870  AssertThrow((last_component + 1 - first_component <= 3),
5871  ExcMessage(
5872  "Can't declare a vector with more than 3 components "
5873  "in VTK/VTU format."));
5874  }
5875 
5876  // write the header. concatenate all the component names with double
5877  // underscores unless a vector name has been specified
5878  o << " <DataArray type=\"Float32\" Name=\"";
5879 
5880  if (!name.empty())
5881  o << name;
5882  else
5883  {
5884  for (unsigned int i = first_component; i < last_component; ++i)
5885  o << data_names[i] << "__";
5886  o << data_names[last_component];
5887  }
5888 
5889  o << "\" NumberOfComponents=\"" << n_components << "\" format=\""
5890  << ascii_or_binary << "\"";
5891  // If present, also list the physical units for this quantity. Look
5892  // this up for either the name of the whole vector/tensor, or if that
5893  // isn't listed, via its first component.
5894  if (!name.empty())
5895  {
5896  if (flags.physical_units.find(name) != flags.physical_units.end())
5897  o << " units=\"" << flags.physical_units.at(name) << "\"";
5898  }
5899  else
5900  {
5901  if (flags.physical_units.find(data_names[first_component]) !=
5902  flags.physical_units.end())
5903  o << " units=\""
5904  << flags.physical_units.at(data_names[first_component]) << "\"";
5905  }
5906  o << ">\n";
5907 
5908  // now write data. pad all vectors to have three components
5909  std::vector<float> data;
5910  data.reserve(n_nodes * n_components);
5911 
5912  for (unsigned int n = 0; n < n_nodes; ++n)
5913  {
5914  if (!is_tensor)
5915  {
5916  switch (last_component - first_component)
5917  {
5918  case 0:
5919  data.push_back(data_vectors(first_component, n));
5920  data.push_back(0);
5921  data.push_back(0);
5922  break;
5923 
5924  case 1:
5925  data.push_back(data_vectors(first_component, n));
5926  data.push_back(data_vectors(first_component + 1, n));
5927  data.push_back(0);
5928  break;
5929 
5930  case 2:
5931  data.push_back(data_vectors(first_component, n));
5932  data.push_back(data_vectors(first_component + 1, n));
5933  data.push_back(data_vectors(first_component + 2, n));
5934  break;
5935 
5936  default:
5937  // Anything else is not yet implemented
5938  Assert(false, ExcInternalError());
5939  }
5940  }
5941  else
5942  {
5943  Tensor<2, 3> vtk_data;
5944  vtk_data = 0.;
5945 
5946  const unsigned int size = last_component - first_component + 1;
5947  if (size == 1)
5948  // 1d, 1 element
5949  {
5950  vtk_data[0][0] = data_vectors(first_component, n);
5951  }
5952  else if (size == 4)
5953  // 2d, 4 elements
5954  {
5955  for (unsigned int c = 0; c < size; ++c)
5956  {
5957  const auto ind =
5959  vtk_data[ind[0]][ind[1]] =
5960  data_vectors(first_component + c, n);
5961  }
5962  }
5963  else if (size == 9)
5964  // 3d 9 elements
5965  {
5966  for (unsigned int c = 0; c < size; ++c)
5967  {
5968  const auto ind =
5970  vtk_data[ind[0]][ind[1]] =
5971  data_vectors(first_component + c, n);
5972  }
5973  }
5974  else
5975  {
5976  Assert(false, ExcInternalError());
5977  }
5978 
5979  // now put the tensor into data
5980  // note we pad with zeros because VTK format always wants to
5981  // see a 3x3 tensor, regardless of dimension
5982  for (unsigned int i = 0; i < 3; ++i)
5983  for (unsigned int j = 0; j < 3; ++j)
5984  data.push_back(vtk_data[i][j]);
5985  }
5986  } // loop over nodes
5987 
5988  o << vtu_stringize_array(data,
5989  flags.compression_level,
5990  output_precision);
5991  o << '\n';
5992  o << " </DataArray>\n";
5993 
5994  return o.str();
5995  };
5996 
5997  const auto stringize_scalar_data_set =
5998  [&flags,
5999  &data_names,
6000  ascii_or_binary,
6001  output_precision = out.precision()](const Table<2, float> &data_vectors,
6002  const unsigned int data_set) {
6003  std::ostringstream o;
6004 
6005  o << " <DataArray type=\"Float32\" Name=\"" << data_names[data_set]
6006  << "\" format=\"" << ascii_or_binary << "\"";
6007  // If present, also list the physical units for this quantity.
6008  if (flags.physical_units.find(data_names[data_set]) !=
6009  flags.physical_units.end())
6010  o << " units=\"" << flags.physical_units.at(data_names[data_set])
6011  << "\"";
6012 
6013  o << ">\n";
6014 
6015  const std::vector<float> data(data_vectors[data_set].begin(),
6016  data_vectors[data_set].end());
6017  o << vtu_stringize_array(data,
6018  flags.compression_level,
6019  output_precision);
6020  o << '\n';
6021  o << " </DataArray>\n";
6022 
6023  return o.str();
6024  };
6025 
6026 
6027  // For the format we write here, we need to write all node values relating
6028  // to one variable at a time. We could in principle do this by looping
6029  // over all patches and extracting the values corresponding to the one
6030  // variable we're dealing with right now, and then start the process over
6031  // for the next variable with another loop over all patches.
6032  //
6033  // An easier way is to create a global table that for each variable
6034  // lists all values. This copying of data vectors can be done in the
6035  // background while we're already working on vertices and cells,
6036  // so do this on a separate task and when wanting to write out the
6037  // data, we wait for that task to finish.
6039  create_global_data_table_task = Threads::new_task([&patches]() {
6040  return create_global_data_table<dim, spacedim, float>(patches);
6041  });
6042 
6043  // -----------------------------
6044  // Now finally get around to actually doing anything. Let's start with
6045  // running the first three tasks generating the vertex and cell information:
6047  mesh_tasks += Threads::new_task(stringize_vertex_information);
6048  mesh_tasks += Threads::new_task(stringize_cell_to_vertex_information);
6049  mesh_tasks += Threads::new_task(stringize_cell_offset_and_type_information);
6050 
6051  // For what follows, we have to have the reordered data available. So wait
6052  // for that task to conclude and get the resulting data table:
6053  const Table<2, float> data_vectors =
6054  std::move(*create_global_data_table_task.return_value());
6055 
6056  // Then create the strings for the actual values of the solution vectors,
6057  // again on separate tasks:
6059  // When writing, first write out all vector and tensor data
6060  std::vector<bool> data_set_handled(n_data_sets, false);
6061  for (const auto &range : nonscalar_data_ranges)
6062  {
6063  // Mark these components as already handled:
6064  const auto first_component = std::get<0>(range);
6065  const auto last_component = std::get<1>(range);
6066  for (unsigned int i = first_component; i <= last_component; ++i)
6067  data_set_handled[i] = true;
6068 
6069  data_tasks += Threads::new_task([&, range]() {
6070  return stringize_nonscalar_data_range(data_vectors, range);
6071  });
6072  }
6073 
6074  // Now do the left over scalar data sets
6075  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6076  if (data_set_handled[data_set] == false)
6077  {
6078  data_tasks += Threads::new_task([&, data_set]() {
6079  return stringize_scalar_data_set(data_vectors, data_set);
6080  });
6081  }
6082 
6083  // Alright, all tasks are now running. Wait for their conclusion and output
6084  // all of the data they have produced:
6085  out << "<Piece NumberOfPoints=\"" << n_nodes << "\" NumberOfCells=\""
6086  << n_cells << "\" >\n";
6087  for (const auto &s : mesh_tasks.return_values())
6088  out << s;
6089  out << " <PointData Scalars=\"scalars\">\n";
6090  for (const auto &s : data_tasks.return_values())
6091  out << s;
6092  out << " </PointData>\n";
6093  out << " </Piece>\n";
6094 
6095  // make sure everything now gets to disk
6096  out.flush();
6097 
6098  // assert the stream is still ok
6099  AssertThrow(out.fail() == false, ExcIO());
6100  }
6101 
6102 
6103 
6104  void
6106  std::ostream &out,
6107  const std::vector<std::string> &piece_names,
6108  const std::vector<std::string> &data_names,
6109  const std::vector<
6110  std::tuple<unsigned int,
6111  unsigned int,
6112  std::string,
6114  &nonscalar_data_ranges,
6115  const VtkFlags &flags)
6116  {
6117  AssertThrow(out.fail() == false, ExcIO());
6118 
6119  // If the user provided physical units, make sure that they don't contain
6120  // quote characters as this would make the VTU file invalid XML and
6121  // probably lead to all sorts of difficult error messages. Other than that,
6122  // trust the user that whatever they provide makes sense somehow.
6123  for (const auto &unit : flags.physical_units)
6124  {
6125  (void)unit;
6126  Assert(
6127  unit.second.find('\"') == std::string::npos,
6128  ExcMessage(
6129  "A physical unit you provided, <" + unit.second +
6130  ">, contained a quotation mark character. This is not allowed."));
6131  }
6132 
6133  const unsigned int n_data_sets = data_names.size();
6134 
6135  out << "<?xml version=\"1.0\"?>\n";
6136 
6137  out << "<!--\n";
6138  out << "#This file was generated by the deal.II library"
6139  << " on " << Utilities::System::get_date() << " at "
6140  << Utilities::System::get_time() << "\n-->\n";
6141 
6142  out
6143  << "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
6144  out << " <PUnstructuredGrid GhostLevel=\"0\">\n";
6145  out << " <PPointData Scalars=\"scalars\">\n";
6146 
6147  // We need to output in the same order as the write_vtu function does:
6148  std::vector<bool> data_set_written(n_data_sets, false);
6149  for (const auto &nonscalar_data_range : nonscalar_data_ranges)
6150  {
6151  const auto first_component = std::get<0>(nonscalar_data_range);
6152  const auto last_component = std::get<1>(nonscalar_data_range);
6153  const bool is_tensor =
6154  (std::get<3>(nonscalar_data_range) ==
6156  const unsigned int n_components = (is_tensor ? 9 : 3);
6157  AssertThrow(last_component >= first_component,
6158  ExcLowerRange(last_component, first_component));
6159  AssertThrow(last_component < n_data_sets,
6160  ExcIndexRange(last_component, 0, n_data_sets));
6161  if (is_tensor)
6162  {
6163  AssertThrow((last_component + 1 - first_component <= 9),
6164  ExcMessage(
6165  "Can't declare a tensor with more than 9 components "
6166  "in VTK"));
6167  }
6168  else
6169  {
6170  Assert((last_component + 1 - first_component <= 3),
6171  ExcMessage(
6172  "Can't declare a vector with more than 3 components "
6173  "in VTK"));
6174  }
6175 
6176  // mark these components as already written:
6177  for (unsigned int i = std::get<0>(nonscalar_data_range);
6178  i <= std::get<1>(nonscalar_data_range);
6179  ++i)
6180  data_set_written[i] = true;
6181 
6182  // write the header. concatenate all the component names with double
6183  // underscores unless a vector name has been specified
6184  out << " <PDataArray type=\"Float32\" Name=\"";
6185 
6186  const std::string &name = std::get<2>(nonscalar_data_range);
6187  if (!name.empty())
6188  out << name;
6189  else
6190  {
6191  for (unsigned int i = std::get<0>(nonscalar_data_range);
6192  i < std::get<1>(nonscalar_data_range);
6193  ++i)
6194  out << data_names[i] << "__";
6195  out << data_names[std::get<1>(nonscalar_data_range)];
6196  }
6197 
6198  out << "\" NumberOfComponents=\"" << n_components
6199  << "\" format=\"ascii\"";
6200  // If present, also list the physical units for this quantity. Look this
6201  // up for either the name of the whole vector/tensor, or if that isn't
6202  // listed, via its first component.
6203  if (!name.empty())
6204  {
6205  if (flags.physical_units.find(name) != flags.physical_units.end())
6206  out << " units=\"" << flags.physical_units.at(name) << "\"";
6207  }
6208  else
6209  {
6210  if (flags.physical_units.find(
6211  data_names[std::get<1>(nonscalar_data_range)]) !=
6212  flags.physical_units.end())
6213  out << " units=\""
6214  << flags.physical_units.at(
6215  data_names[std::get<1>(nonscalar_data_range)])
6216  << "\"";
6217  }
6218 
6219  out << "/>\n";
6220  }
6221 
6222  // Now for the scalar fields
6223  for (unsigned int data_set = 0; data_set < n_data_sets; ++data_set)
6224  if (data_set_written[data_set] == false)
6225  {
6226  out << " <PDataArray type=\"Float32\" Name=\""
6227  << data_names[data_set] << "\" format=\"ascii\"";
6228 
6229  if (flags.physical_units.find(data_names[data_set]) !=
6230  flags.physical_units.end())
6231  out << " units=\"" << flags.physical_units.at(data_names[data_set])
6232  << "\"";
6233 
6234  out << "/>\n";
6235  }
6236 
6237  out << " </PPointData>\n";
6238 
6239  out << " <PPoints>\n";
6240  out << " <PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n";
6241  out << " </PPoints>\n";
6242 
6243  for (const auto &piece_name : piece_names)
6244  out << " <Piece Source=\"" << piece_name << "\"/>\n";
6245 
6246  out << " </PUnstructuredGrid>\n";
6247  out << "</VTKFile>\n";
6248 
6249  out.flush();
6250 
6251  // assert the stream is still ok
6252  AssertThrow(out.fail() == false, ExcIO());
6253  }
6254 
6255 
6256 
6257  void
6259  std::ostream &out,
6260  const std::vector<std::pair<double, std::string>> &times_and_names)
6261  {
6262  AssertThrow(out.fail() == false, ExcIO());
6263 
6264  out << "<?xml version=\"1.0\"?>\n";
6265 
6266  out << "<!--\n";
6267  out << "#This file was generated by the deal.II library"
6268  << " on " << Utilities::System::get_date() << " at "
6269  << Utilities::System::get_time() << "\n-->\n";
6270 
6271  out
6272  << "<VTKFile type=\"Collection\" version=\"0.1\" ByteOrder=\"LittleEndian\">\n";
6273  out << " <Collection>\n";
6274 
6275  std::streamsize ss = out.precision();
6276  out.precision(12);
6277 
6278  for (const auto &time_and_name : times_and_names)
6279  out << " <DataSet timestep=\"" << time_and_name.first
6280  << "\" group=\"\" part=\"0\" file=\"" << time_and_name.second
6281  << "\"/>\n";
6282 
6283  out << " </Collection>\n";
6284  out << "</VTKFile>\n";
6285 
6286  out.flush();
6287  out.precision(ss);
6288 
6289  AssertThrow(out.fail() == false, ExcIO());
6290  }
6291 
6292 
6293 
6294  void
6295  write_visit_record(std::ostream &out,
6296  const std::vector<std::string> &piece_names)
6297  {
6298  out << "!NBLOCKS " << piece_names.size() << '\n';
6299  for (const auto &piece_name : piece_names)
6300  out << piece_name << '\n';
6301 
6302  out << std::flush;
6303  }
6304 
6305 
6306 
6307  void
6308  write_visit_record(std::ostream &out,
6309  const std::vector<std::vector<std::string>> &piece_names)
6310  {
6311  AssertThrow(out.fail() == false, ExcIO());
6312 
6313  if (piece_names.empty())
6314  return;
6315 
6316  const double nblocks = piece_names[0].size();
6317  Assert(nblocks > 0,
6318  ExcMessage("piece_names should be a vector of nonempty vectors."));
6319 
6320  out << "!NBLOCKS " << nblocks << '\n';
6321  for (const auto &domain : piece_names)
6322  {
6323  Assert(domain.size() == nblocks,
6324  ExcMessage(
6325  "piece_names should be a vector of equal sized vectors."));
6326  for (const auto &subdomain : domain)
6327  out << subdomain << '\n';
6328  }
6329 
6330  out << std::flush;
6331  }
6332 
6333 
6334 
6335  void
6337  std::ostream &out,
6338  const std::vector<std::pair<double, std::vector<std::string>>>
6339  &times_and_piece_names)
6340  {
6341  AssertThrow(out.fail() == false, ExcIO());
6342 
6343  if (times_and_piece_names.empty())
6344  return;
6345 
6346  const double nblocks = times_and_piece_names[0].second.size();
6347  Assert(
6348