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