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