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