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