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