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