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