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