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