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