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