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