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