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