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