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