Reference documentation for deal.II version Git cd5f69ed86 2019-12-07 22:13:35 -0500
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
grid_out.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/exceptions.h>
17 #include <deal.II/base/geometry_info.h>
18 #include <deal.II/base/parameter_handler.h>
19 #include <deal.II/base/point.h>
20 #include <deal.II/base/qprojector.h>
21 #include <deal.II/base/quadrature.h>
22 
23 #include <deal.II/fe/mapping.h>
24 
25 #include <deal.II/grid/grid_out.h>
26 #include <deal.II/grid/tria.h>
27 #include <deal.II/grid/tria_accessor.h>
28 #include <deal.II/grid/tria_iterator.h>
29 
30 #include <deal.II/numerics/data_out.h>
31 
32 #include <algorithm>
33 #include <cmath>
34 #include <cstring>
35 #include <ctime>
36 #include <fstream>
37 #include <iomanip>
38 #include <list>
39 #include <set>
40 
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 
45 namespace GridOutFlags
46 {
47  DX::DX(const bool write_cells,
48  const bool write_faces,
49  const bool write_diameter,
50  const bool write_measure,
51  const bool write_all_faces)
52  : write_cells(write_cells)
53  , write_faces(write_faces)
54  , write_diameter(write_diameter)
55  , write_measure(write_measure)
56  , write_all_faces(write_all_faces)
57  {}
58 
59  void
61  {
62  param.declare_entry("Write cells",
63  "true",
65  "Write the mesh connectivity as DX grid cells");
66  param.declare_entry("Write faces",
67  "false",
69  "Write faces of cells. These may be boundary faces "
70  "or all faces between mesh cells, according to "
71  "\"Write all faces\"");
72  param.declare_entry("Write diameter",
73  "false",
75  "If cells are written, additionally write their"
76  " diameter as data for visualization");
77  param.declare_entry("Write measure",
78  "false",
80  "Write the volume of each cell as data");
81  param.declare_entry("Write all faces",
82  "true",
84  "Write all faces, not only boundary");
85  }
86 
87  void
89  {
90  write_cells = param.get_bool("Write cells");
91  write_faces = param.get_bool("Write faces");
92  write_diameter = param.get_bool("Write diameter");
93  write_measure = param.get_bool("Write measure");
94  write_all_faces = param.get_bool("Write all faces");
95  }
96 
97 
98  Msh::Msh(const bool write_faces, const bool write_lines)
99  : write_faces(write_faces)
100  , write_lines(write_lines)
101  {}
102 
103  void
105  {
106  param.declare_entry("Write faces", "false", Patterns::Bool());
107  param.declare_entry("Write lines", "false", Patterns::Bool());
108  }
109 
110 
111  void
113  {
114  write_faces = param.get_bool("Write faces");
115  write_lines = param.get_bool("Write lines");
116  }
117 
118 
119  Ucd::Ucd(const bool write_preamble,
120  const bool write_faces,
121  const bool write_lines)
122  : write_preamble(write_preamble)
123  , write_faces(write_faces)
124  , write_lines(write_lines)
125  {}
126 
127 
128 
129  void
131  {
132  param.declare_entry("Write preamble", "true", Patterns::Bool());
133  param.declare_entry("Write faces", "false", Patterns::Bool());
134  param.declare_entry("Write lines", "false", Patterns::Bool());
135  }
136 
137 
138  void
140  {
141  write_preamble = param.get_bool("Write preamble");
142  write_faces = param.get_bool("Write faces");
143  write_lines = param.get_bool("Write lines");
144  }
145 
146 
147 
148  Gnuplot::Gnuplot(const bool write_cell_numbers,
149  const unsigned int n_extra_curved_line_points,
150  const bool curved_inner_cells,
151  const bool write_additional_boundary_lines)
152  : write_cell_numbers(write_cell_numbers)
153  , n_extra_curved_line_points(n_extra_curved_line_points)
154  , n_boundary_face_points(this->n_extra_curved_line_points)
155  , curved_inner_cells(curved_inner_cells)
156  , write_additional_boundary_lines(write_additional_boundary_lines)
157  {}
158 
159 
160  // TODO we can get rid of these extra constructors and assignment operators
161  // once we remove the reference member variable.
162  Gnuplot::Gnuplot(const Gnuplot &flags)
163  : Gnuplot(flags.write_cell_numbers,
165  flags.curved_inner_cells,
167  {}
168 
169 
170 
171  Gnuplot &
173  {
178 
179  return *this;
180  }
181 
182 
183 
184  void
186  {
187  param.declare_entry("Cell number", "false", Patterns::Bool());
188  param.declare_entry("Boundary points", "2", Patterns::Integer());
189  }
190 
191 
192  void
194  {
195  write_cell_numbers = param.get_bool("Cell number");
196  n_boundary_face_points = param.get_integer("Boundary points");
197  }
198 
199 
201  const unsigned int size,
202  const double line_width,
203  const bool color_lines_on_user_flag,
204  const unsigned int n_boundary_face_points,
205  const bool color_lines_level)
206  : size_type(size_type)
207  , size(size)
208  , line_width(line_width)
209  , color_lines_on_user_flag(color_lines_on_user_flag)
210  , n_boundary_face_points(n_boundary_face_points)
211  , color_lines_level(color_lines_level)
212  {}
213 
214 
215  void
217  {
218  param.declare_entry("Size by",
219  "width",
220  Patterns::Selection("width|height"),
221  "Depending on this parameter, either the "
222  "width or height "
223  "of the eps is scaled to \"Size\"");
224  param.declare_entry("Size",
225  "300",
227  "Size of the output in points");
228  param.declare_entry("Line width",
229  "0.5",
231  "Width of the lines drawn in points");
232  param.declare_entry("Color by flag",
233  "false",
234  Patterns::Bool(),
235  "Draw lines with user flag set in different color");
236  param.declare_entry("Boundary points",
237  "2",
239  "Number of points on boundary edges. "
240  "Increase this beyond 2 to see curved boundaries.");
241  param.declare_entry("Color by level",
242  "false",
243  Patterns::Bool(),
244  "Draw different colors according to grid level.");
245  }
246 
247 
248  void
250  {
251  if (param.get("Size by") == std::string("width"))
252  size_type = width;
253  else if (param.get("Size by") == std::string("height"))
254  size_type = height;
255  size = param.get_integer("Size");
256  line_width = param.get_double("Line width");
257  color_lines_on_user_flag = param.get_bool("Color by flag");
258  n_boundary_face_points = param.get_integer("Boundary points");
259  color_lines_level = param.get_bool("Color by level");
260  }
261 
262 
263 
265  const unsigned int size,
266  const double line_width,
267  const bool color_lines_on_user_flag,
268  const unsigned int n_boundary_face_points)
269  : EpsFlagsBase(size_type,
270  size,
271  line_width,
272  color_lines_on_user_flag,
273  n_boundary_face_points)
274  {}
275 
276 
277  void
279  {}
280 
281 
282  void
284  {
286  }
287 
288 
289 
290  Eps<2>::Eps(const SizeType size_type,
291  const unsigned int size,
292  const double line_width,
293  const bool color_lines_on_user_flag,
294  const unsigned int n_boundary_face_points,
295  const bool write_cell_numbers,
296  const bool write_cell_number_level,
297  const bool write_vertex_numbers,
298  const bool color_lines_level)
299  : EpsFlagsBase(size_type,
300  size,
301  line_width,
302  color_lines_on_user_flag,
303  n_boundary_face_points,
304  color_lines_level)
305  , write_cell_numbers(write_cell_numbers)
306  , write_cell_number_level(write_cell_number_level)
307  , write_vertex_numbers(write_vertex_numbers)
308  {}
309 
310 
311  void
313  {
314  param.declare_entry("Cell number",
315  "false",
316  Patterns::Bool(),
317  "(2D only) Write cell numbers"
318  " into the centers of cells");
319  param.declare_entry("Level number",
320  "false",
321  Patterns::Bool(),
322  "(2D only) if \"Cell number\" is true, write "
323  "numbers in the form level.number");
324  param.declare_entry("Vertex number",
325  "false",
326  Patterns::Bool(),
327  "Write numbers for each vertex");
328  }
329 
330 
331  void
333  {
335  write_cell_numbers = param.get_bool("Cell number");
336  write_cell_number_level = param.get_bool("Level number");
337  write_vertex_numbers = param.get_bool("Vertex number");
338  }
339 
340 
341 
342  Eps<3>::Eps(const SizeType size_type,
343  const unsigned int size,
344  const double line_width,
345  const bool color_lines_on_user_flag,
346  const unsigned int n_boundary_face_points,
347  const double azimut_angle,
348  const double turn_angle)
349  : EpsFlagsBase(size_type,
350  size,
351  line_width,
352  color_lines_on_user_flag,
353  n_boundary_face_points)
354  , azimut_angle(azimut_angle)
355  , turn_angle(turn_angle)
356  {}
357 
358 
359  void
361  {
362  param.declare_entry("Azimuth",
363  "30",
365  "Azimuth of the viw point, that is, the angle "
366  "in the plane from the x-axis.");
367  param.declare_entry("Elevation",
368  "30",
370  "Elevation of the view point above the xy-plane.");
371  }
372 
373 
374  void
376  {
378  azimut_angle = 90 - param.get_double("Elevation");
379  turn_angle = param.get_double("Azimuth");
380  }
381 
382 
383 
385  : draw_boundary(true)
386  , color_by(material_id)
387  , level_depth(true)
388  , n_boundary_face_points(0)
389  , scaling(1., 1.)
390  , fill_style(20)
391  , line_style(0)
392  , line_thickness(1)
393  , boundary_style(0)
394  , boundary_thickness(3)
395  {}
396 
397 
398  void
400  {
401  param.declare_entry("Boundary", "true", Patterns::Bool());
402  param.declare_entry("Level color", "false", Patterns::Bool());
403  param.declare_entry("Level depth", "true", Patterns::Bool());
404  // TODO: Unify this number with other output formats
405  param.declare_entry("Boundary points", "0", Patterns::Integer());
406  param.declare_entry("Fill style", "20", Patterns::Integer());
407  param.declare_entry("Line style", "0", Patterns::Integer());
408  param.declare_entry("Line width", "1", Patterns::Integer());
409  param.declare_entry("Boundary style", "0", Patterns::Integer());
410  param.declare_entry("Boundary width", "3", Patterns::Integer());
411  }
412 
413 
414  void
416  {
417  draw_boundary = param.get_bool("Boundary");
418  level_depth = param.get_bool("Level depth");
419  n_boundary_face_points = param.get_integer("Boundary points");
420  fill_style = param.get_integer("Fill style");
421  line_style = param.get_integer("Line style");
422  line_thickness = param.get_integer("Line width");
423  boundary_style = param.get_integer("Boundary style");
424  boundary_thickness = param.get_integer("Boundary width");
425  }
426 
427  Svg::Svg(const unsigned int line_thickness,
428  const unsigned int boundary_line_thickness,
429  bool margin,
430  const Background background,
431  const int azimuth_angle,
432  const int polar_angle,
433  const Coloring coloring,
434  const bool convert_level_number_to_height,
435  const bool label_level_number,
436  const bool label_cell_index,
437  const bool label_material_id,
438  const bool label_subdomain_id,
439  const bool draw_colorbar,
440  const bool draw_legend)
441  : height(1000)
442  , width(0)
443  , line_thickness(line_thickness)
444  , boundary_line_thickness(boundary_line_thickness)
445  , margin(margin)
446  , background(background)
447  , azimuth_angle(azimuth_angle)
448  , polar_angle(polar_angle)
449  , coloring(coloring)
450  , convert_level_number_to_height(convert_level_number_to_height)
451  , level_height_factor(0.3f)
452  , cell_font_scaling(1.f)
453  , label_level_number(label_level_number)
454  , label_cell_index(label_cell_index)
455  , label_material_id(label_material_id)
456  , label_subdomain_id(label_subdomain_id)
457  , label_level_subdomain_id(false)
458  , draw_colorbar(draw_colorbar)
459  , draw_legend(draw_legend)
460  {}
461 
463  : draw_bounding_box(false) // box
464  {}
465 
466  void
468  {
469  param.declare_entry("Draw bounding box", "false", Patterns::Bool());
470  }
471 
472  void
474  {
475  draw_bounding_box = param.get_bool("Draw bounding box");
476  }
477 } // end namespace GridOutFlags
478 
479 
480 
482  : default_format(none)
483 {}
484 
485 
486 void
488 {
489  dx_flags = flags;
490 }
491 
492 
493 
494 void
496 {
497  msh_flags = flags;
498 }
499 
500 
501 void
503 {
504  ucd_flags = flags;
505 }
506 
507 
508 
509 void
511 {
512  gnuplot_flags = flags;
513 }
514 
515 
516 
517 void
519 {
520  eps_flags_1 = flags;
521 }
522 
523 
524 
525 void
527 {
528  eps_flags_2 = flags;
529 }
530 
531 
532 
533 void
535 {
536  eps_flags_3 = flags;
537 }
538 
539 
540 
541 void
543 {
544  xfig_flags = flags;
545 }
546 
547 
548 void
550 {
551  svg_flags = flags;
552 }
553 
554 
555 void
557 {
558  mathgl_flags = flags;
559 }
560 
561 void
563 {
564  vtk_flags = flags;
565 }
566 
567 void
569 {
570  vtu_flags = flags;
571 }
572 
573 std::string
575 {
576  switch (output_format)
577  {
578  case none:
579  return "";
580  case dx:
581  return ".dx";
582  case gnuplot:
583  return ".gnuplot";
584  case ucd:
585  return ".inp";
586  case eps:
587  return ".eps";
588  case xfig:
589  return ".fig";
590  case msh:
591  return ".msh";
592  case svg:
593  return ".svg";
594  case mathgl:
595  return ".mathgl";
596  case vtk:
597  return ".vtk";
598  case vtu:
599  return ".vtu";
600  default:
601  Assert(false, ExcNotImplemented());
602  return "";
603  }
604 }
605 
606 
607 
608 std::string
610 {
612 }
613 
614 
615 
617 GridOut::parse_output_format(const std::string &format_name)
618 {
619  if (format_name == "none" || format_name == "false")
620  return none;
621 
622  if (format_name == "dx")
623  return dx;
624 
625  if (format_name == "ucd")
626  return ucd;
627 
628  if (format_name == "gnuplot")
629  return gnuplot;
630 
631  if (format_name == "eps")
632  return eps;
633 
634  if (format_name == "xfig")
635  return xfig;
636 
637  if (format_name == "msh")
638  return msh;
639 
640  if (format_name == "svg")
641  return svg;
642 
643  if (format_name == "mathgl")
644  return mathgl;
645 
646  if (format_name == "vtk")
647  return vtk;
648 
649  if (format_name == "vtu")
650  return vtu;
651 
652  AssertThrow(false, ExcInvalidState());
653  // return something weird
654  return OutputFormat(-1);
655 }
656 
657 
658 
659 std::string
661 {
662  return "none|dx|gnuplot|eps|ucd|xfig|msh|svg|mathgl|vtk|vtu";
663 }
664 
665 
666 void
668 {
669  param.declare_entry("Format",
670  "none",
672 
673  param.enter_subsection("DX");
675  param.leave_subsection();
676 
677  param.enter_subsection("Msh");
679  param.leave_subsection();
680 
681  param.enter_subsection("Ucd");
683  param.leave_subsection();
684 
685  param.enter_subsection("Gnuplot");
687  param.leave_subsection();
688 
689  param.enter_subsection("Eps");
694  param.leave_subsection();
695 
696  param.enter_subsection("XFig");
698  param.leave_subsection();
699 
700  param.enter_subsection("MathGL");
702  param.leave_subsection();
703 
704  param.enter_subsection("Vtk");
706  param.leave_subsection();
707 
708  param.enter_subsection("Vtu");
710  param.leave_subsection();
711 }
712 
713 
714 
715 void
717 {
718  default_format = parse_output_format(param.get("Format"));
719 
720  param.enter_subsection("DX");
721  dx_flags.parse_parameters(param);
722  param.leave_subsection();
723 
724  param.enter_subsection("Msh");
726  param.leave_subsection();
727 
728  param.enter_subsection("Ucd");
730  param.leave_subsection();
731 
732  param.enter_subsection("Gnuplot");
734  param.leave_subsection();
735 
736  param.enter_subsection("Eps");
740  param.leave_subsection();
741 
742  param.enter_subsection("XFig");
744  param.leave_subsection();
745 
746  param.enter_subsection("MathGL");
748  param.leave_subsection();
749 
750  param.enter_subsection("Vtk");
752  param.leave_subsection();
753 
754  param.enter_subsection("Vtu");
756  param.leave_subsection();
757 }
758 
759 
760 
761 std::size_t
763 {
764  return (sizeof(dx_flags) + sizeof(msh_flags) + sizeof(ucd_flags) +
765  sizeof(gnuplot_flags) + sizeof(eps_flags_1) + sizeof(eps_flags_2) +
766  sizeof(eps_flags_3) + sizeof(xfig_flags) + sizeof(svg_flags) +
767  sizeof(mathgl_flags) + sizeof(vtk_flags) + sizeof(vtu_flags));
768 }
769 
770 
771 
772 template <>
773 void
774 GridOut::write_dx(const Triangulation<1> &, std::ostream &) const
775 {
776  Assert(false, ExcNotImplemented());
777 }
778 
779 template <>
780 void
781 GridOut::write_dx(const Triangulation<1, 2> &, std::ostream &) const
782 {
783  Assert(false, ExcNotImplemented());
784 }
785 
786 template <>
787 void
788 GridOut::write_dx(const Triangulation<1, 3> &, std::ostream &) const
789 {
790  Assert(false, ExcNotImplemented());
791 }
792 
793 
794 
795 template <int dim, int spacedim>
796 void
798  std::ostream & out) const
799 {
800  // TODO:[GK] allow for boundary faces only
802  AssertThrow(out, ExcIO());
803  // Copied and adapted from write_ucd
804  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
805  const std::vector<bool> & vertex_used = tria.get_used_vertices();
806 
807  const unsigned int n_vertices = tria.n_used_vertices();
808 
809  // vertices are implicitly numbered from 0 to
810  // n_vertices-1. we have to renumber the
811  // vertices, because otherwise we would end
812  // up with wrong results, if there are unused
813  // vertices
814  std::vector<unsigned int> renumber(vertices.size());
815  // fill this vector with new vertex numbers
816  // ranging from 0 to n_vertices-1
817  unsigned int new_number = 0;
818  for (unsigned int i = 0; i < vertices.size(); ++i)
819  if (vertex_used[i])
820  renumber[i] = new_number++;
821  Assert(new_number == n_vertices, ExcInternalError());
822 
825  tria.end();
826 
827 
828  // write the vertices
829  out << "object \"vertices\" class array type float rank 1 shape " << dim
830  << " items " << n_vertices << " data follows" << '\n';
831 
832  for (unsigned int i = 0; i < vertices.size(); ++i)
833  if (vertex_used[i])
834  out << '\t' << vertices[i] << '\n';
835 
836  // write cells or faces
837  const bool write_cells = dx_flags.write_cells;
838  const bool write_faces = (dim > 1) ? dx_flags.write_faces : false;
839 
840  const unsigned int n_cells = tria.n_active_cells();
841  const unsigned int n_faces =
843 
844  const unsigned int n_vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
845  const unsigned int n_vertices_per_face = GeometryInfo<dim>::vertices_per_face;
846 
847  if (write_cells)
848  {
849  out << "object \"cells\" class array type int rank 1 shape "
850  << n_vertices_per_cell << " items " << n_cells << " data follows"
851  << '\n';
852 
853  for (cell = tria.begin_active(); cell != endc; ++cell)
854  {
855  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
856  ++v)
857  out
858  << '\t'
859  << renumber[cell->vertex_index(GeometryInfo<dim>::dx_to_deal[v])];
860  out << '\n';
861  }
862  out << "attribute \"element type\" string \"";
863  if (dim == 1)
864  out << "lines";
865  if (dim == 2)
866  out << "quads";
867  if (dim == 3)
868  out << "cubes";
869  out << "\"" << '\n'
870  << "attribute \"ref\" string \"positions\"" << '\n'
871  << '\n';
872 
873  // Additional cell information
874 
875  out << "object \"material\" class array type int rank 0 items " << n_cells
876  << " data follows" << '\n';
877  for (cell = tria.begin_active(); cell != endc; ++cell)
878  out << ' ' << cell->material_id();
879  out << '\n' << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
880 
881  out << "object \"level\" class array type int rank 0 items " << n_cells
882  << " data follows" << '\n';
883  for (cell = tria.begin_active(); cell != endc; ++cell)
884  out << ' ' << cell->level();
885  out << '\n' << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
886 
888  {
889  out << "object \"measure\" class array type float rank 0 items "
890  << n_cells << " data follows" << '\n';
891  for (cell = tria.begin_active(); cell != endc; ++cell)
892  out << '\t' << cell->measure();
893  out << '\n'
894  << "attribute \"dep\" string \"connections\"" << '\n'
895  << '\n';
896  }
897 
899  {
900  out << "object \"diameter\" class array type float rank 0 items "
901  << n_cells << " data follows" << '\n';
902  for (cell = tria.begin_active(); cell != endc; ++cell)
903  out << '\t' << cell->diameter();
904  out << '\n'
905  << "attribute \"dep\" string \"connections\"" << '\n'
906  << '\n';
907  }
908  }
909 
910  if (write_faces)
911  {
912  out << "object \"faces\" class array type int rank 1 shape "
913  << n_vertices_per_face << " items " << n_faces << " data follows"
914  << '\n';
915 
916  for (cell = tria.begin_active(); cell != endc; ++cell)
917  {
918  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
919  {
921  cell->face(f);
922 
923  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
924  ++v)
925  out << '\t'
926  << renumber[face->vertex_index(
928  out << '\n';
929  }
930  }
931  out << "attribute \"element type\" string \"";
932  if (dim == 2)
933  out << "lines";
934  if (dim == 3)
935  out << "quads";
936  out << "\"" << '\n'
937  << "attribute \"ref\" string \"positions\"" << '\n'
938  << '\n';
939 
940 
941  // Additional face information
942 
943  out << "object \"boundary\" class array type int rank 0 items " << n_faces
944  << " data follows" << '\n';
945  for (cell = tria.begin_active(); cell != endc; ++cell)
946  {
947  // Little trick to get -1 for the interior
948  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
949  {
950  out << ' '
951  << static_cast<std::make_signed<types::boundary_id>::type>(
952  cell->face(f)->boundary_id());
953  }
954  out << '\n';
955  }
956  out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
957 
959  {
960  out << "object \"face measure\" class array type float rank 0 items "
961  << n_faces << " data follows" << '\n';
962  for (cell = tria.begin_active(); cell != endc; ++cell)
963  {
964  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
965  ++f)
966  out << ' ' << cell->face(f)->measure();
967  out << '\n';
968  }
969  out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
970  }
971 
973  {
974  out << "object \"face diameter\" class array type float rank 0 items "
975  << n_faces << " data follows" << '\n';
976  for (cell = tria.begin_active(); cell != endc; ++cell)
977  {
978  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell;
979  ++f)
980  out << ' ' << cell->face(f)->diameter();
981  out << '\n';
982  }
983  out << "attribute \"dep\" string \"connections\"" << '\n' << '\n';
984  }
985  }
986 
987 
988  // Write additional face information
989 
990  if (write_faces)
991  {
992  }
993  else
994  {}
995 
996  // The wrapper
997  out << "object \"deal data\" class field" << '\n'
998  << "component \"positions\" value \"vertices\"" << '\n'
999  << "component \"connections\" value \"cells\"" << '\n';
1000 
1001  if (write_cells)
1002  {
1003  out << "object \"cell data\" class field" << '\n'
1004  << "component \"positions\" value \"vertices\"" << '\n'
1005  << "component \"connections\" value \"cells\"" << '\n';
1006  out << "component \"material\" value \"material\"" << '\n';
1007  out << "component \"level\" value \"level\"" << '\n';
1008  if (dx_flags.write_measure)
1009  out << "component \"measure\" value \"measure\"" << '\n';
1011  out << "component \"diameter\" value \"diameter\"" << '\n';
1012  }
1013 
1014  if (write_faces)
1015  {
1016  out << "object \"face data\" class field" << '\n'
1017  << "component \"positions\" value \"vertices\"" << '\n'
1018  << "component \"connections\" value \"faces\"" << '\n';
1019  out << "component \"boundary\" value \"boundary\"" << '\n';
1020  if (dx_flags.write_measure)
1021  out << "component \"measure\" value \"face measure\"" << '\n';
1023  out << "component \"diameter\" value \"face diameter\"" << '\n';
1024  }
1025 
1026  out << '\n' << "object \"grid data\" class group" << '\n';
1027  if (write_cells)
1028  out << "member \"cells\" value \"cell data\"" << '\n';
1029  if (write_faces)
1030  out << "member \"faces\" value \"face data\"" << '\n';
1031  out << "end" << '\n';
1032 
1033  // make sure everything now gets to
1034  // disk
1035  out.flush();
1036 
1037  AssertThrow(out, ExcIO());
1038 }
1039 
1040 
1041 
1042 template <int dim, int spacedim>
1043 void
1045  std::ostream & out) const
1046 {
1047  AssertThrow(out, ExcIO());
1048 
1049  // get the positions of the
1050  // vertices and whether they are
1051  // used.
1052  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1053  const std::vector<bool> & vertex_used = tria.get_used_vertices();
1054 
1055  const unsigned int n_vertices = tria.n_used_vertices();
1056 
1058  tria.begin_active();
1060  tria.end();
1061 
1062  // Write Header
1063  // The file format is:
1064  /*
1065 
1066 
1067  @f$NOD
1068  number-of-nodes
1069  node-number x-coord y-coord z-coord
1070  ...
1071  @f$ENDNOD
1072  @f$ELM
1073  number-of-elements
1074  elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
1075  ...
1076  @f$ENDELM
1077  */
1078  out << "@f$NOD" << '\n' << n_vertices << '\n';
1079 
1080  // actually write the vertices.
1081  // note that we shall number them
1082  // with first index 1 instead of 0
1083  for (unsigned int i = 0; i < vertices.size(); ++i)
1084  if (vertex_used[i])
1085  {
1086  out << i + 1 // vertex index
1087  << " " << vertices[i];
1088  for (unsigned int d = spacedim + 1; d <= 3; ++d)
1089  out << " 0"; // fill with zeroes
1090  out << '\n';
1091  }
1092 
1093  // Write cells preamble
1094  out << "@f$ENDNOD" << '\n'
1095  << "@f$ELM" << '\n'
1096  << tria.n_active_cells() +
1097  ((msh_flags.write_faces ? n_boundary_faces(tria) : 0) +
1098  (msh_flags.write_lines ? n_boundary_lines(tria) : 0))
1099  << '\n';
1100 
1101  /*
1102  elm-type
1103  defines the geometrical type of the n-th element:
1104  1
1105  Line (2 nodes).
1106  2
1107  Triangle (3 nodes).
1108  3
1109  Quadrangle (4 nodes).
1110  4
1111  Tetrahedron (4 nodes).
1112  5
1113  Hexahedron (8 nodes).
1114  6
1115  Prism (6 nodes).
1116  7
1117  Pyramid (5 nodes).
1118  8
1119  Second order line (3 nodes: 2 associated with the vertices and 1 with the
1120  edge).
1121  9
1122  Second order triangle (6 nodes: 3 associated with the vertices and 3 with
1123  the edges). 10 Second order quadrangle (9 nodes: 4 associated with the
1124  vertices, 4 with the edges and 1 with the face). 11 Second order tetrahedron
1125  (10 nodes: 4 associated with the vertices and 6 with the edges). 12 Second
1126  order hexahedron (27 nodes: 8 associated with the vertices, 12 with the
1127  edges, 6 with the faces and 1 with the volume). 13 Second order prism (18
1128  nodes: 6 associated with the vertices, 9 with the edges and 3 with the
1129  quadrangular faces). 14 Second order pyramid (14 nodes: 5 associated with
1130  the vertices, 8 with the edges and 1 with the quadrangular face). 15 Point
1131  (1 node).
1132  */
1133  unsigned int elm_type;
1134  switch (dim)
1135  {
1136  case 1:
1137  elm_type = 1;
1138  break;
1139  case 2:
1140  elm_type = 3;
1141  break;
1142  case 3:
1143  elm_type = 5;
1144  break;
1145  default:
1146  Assert(false, ExcNotImplemented());
1147  }
1148 
1149  // write cells. Enumerate cells
1150  // consecutively, starting with 1
1151  for (cell = tria.begin_active(); cell != endc; ++cell)
1152  {
1153  out << cell->active_cell_index() + 1 << ' ' << elm_type << ' '
1154  << static_cast<unsigned int>(cell->material_id()) << ' '
1155  << cell->subdomain_id() << ' ' << GeometryInfo<dim>::vertices_per_cell
1156  << ' ';
1157 
1158  // Vertex numbering follows UCD conventions.
1159 
1160  for (unsigned int vertex = 0;
1161  vertex < GeometryInfo<dim>::vertices_per_cell;
1162  ++vertex)
1163  out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) + 1
1164  << ' ';
1165  out << '\n';
1166  }
1167 
1168  // write faces and lines with non-zero boundary indicator
1169  unsigned int next_element_index = tria.n_active_cells() + 1;
1170  if (msh_flags.write_faces)
1171  {
1172  next_element_index = write_msh_faces(tria, next_element_index, out);
1173  }
1174  if (msh_flags.write_lines)
1175  {
1176  next_element_index = write_msh_lines(tria, next_element_index, out);
1177  }
1178 
1179  out << "@f$ENDELM\n";
1180 
1181  // make sure everything now gets to
1182  // disk
1183  out.flush();
1184 
1185  AssertThrow(out, ExcIO());
1186 }
1187 
1188 
1189 template <int dim, int spacedim>
1190 void
1192  std::ostream & out) const
1193 {
1194  AssertThrow(out, ExcIO());
1195 
1196  // get the positions of the
1197  // vertices and whether they are
1198  // used.
1199  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1200  const std::vector<bool> & vertex_used = tria.get_used_vertices();
1201 
1202  const unsigned int n_vertices = tria.n_used_vertices();
1203 
1205  tria.begin_active();
1207  tria.end();
1208 
1209  // write preamble
1211  {
1212  // block this to have local
1213  // variables destroyed after
1214  // use
1215  std::time_t time1 = std::time(nullptr);
1216  std::tm * time = std::localtime(&time1);
1217  out
1218  << "# This file was generated by the deal.II library." << '\n'
1219  << "# Date = " << time->tm_year + 1900 << "/" << time->tm_mon + 1
1220  << "/" << time->tm_mday << '\n'
1221  << "# Time = " << time->tm_hour << ":" << std::setw(2) << time->tm_min
1222  << ":" << std::setw(2) << time->tm_sec << '\n'
1223  << "#" << '\n'
1224  << "# For a description of the UCD format see the AVS Developer's guide."
1225  << '\n'
1226  << "#" << '\n';
1227  }
1228 
1229  // start with ucd data
1230  out << n_vertices << ' '
1231  << tria.n_active_cells() +
1232  ((ucd_flags.write_faces ? n_boundary_faces(tria) : 0) +
1233  (ucd_flags.write_lines ? n_boundary_lines(tria) : 0))
1234  << " 0 0 0" // no data
1235  << '\n';
1236 
1237  // actually write the vertices.
1238  // note that we shall number them
1239  // with first index 1 instead of 0
1240  for (unsigned int i = 0; i < vertices.size(); ++i)
1241  if (vertex_used[i])
1242  {
1243  out << i + 1 // vertex index
1244  << " " << vertices[i];
1245  for (unsigned int d = spacedim + 1; d <= 3; ++d)
1246  out << " 0"; // fill with zeroes
1247  out << '\n';
1248  }
1249 
1250  // write cells. Enumerate cells
1251  // consecutively, starting with 1
1252  for (cell = tria.begin_active(); cell != endc; ++cell)
1253  {
1254  out << cell->active_cell_index() + 1 << ' '
1255  << static_cast<unsigned int>(cell->material_id()) << ' ';
1256  switch (dim)
1257  {
1258  case 1:
1259  out << "line ";
1260  break;
1261  case 2:
1262  out << "quad ";
1263  break;
1264  case 3:
1265  out << "hex ";
1266  break;
1267  default:
1268  Assert(false, ExcNotImplemented());
1269  }
1270 
1271  // it follows a list of the
1272  // vertices of each cell. in 1d
1273  // this is simply a list of the
1274  // two vertices, in 2d its counter
1275  // clockwise, as usual in this
1276  // library. in 3d, the same applies
1277  // (special thanks to AVS for
1278  // numbering their vertices in a
1279  // way compatible to deal.II!)
1280  //
1281  // technical reference:
1282  // AVS Developer's Guide, Release 4,
1283  // May, 1992, p. E6
1284  //
1285  // note: vertex numbers are 1-base
1286  for (unsigned int vertex = 0;
1287  vertex < GeometryInfo<dim>::vertices_per_cell;
1288  ++vertex)
1289  out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) + 1
1290  << ' ';
1291  out << '\n';
1292  }
1293 
1294  // write faces and lines with non-zero boundary indicator
1295  unsigned int next_element_index = tria.n_active_cells() + 1;
1296  if (ucd_flags.write_faces)
1297  {
1298  next_element_index = write_ucd_faces(tria, next_element_index, out);
1299  }
1300  if (ucd_flags.write_lines)
1301  {
1302  next_element_index = write_ucd_lines(tria, next_element_index, out);
1303  }
1304 
1305  // make sure everything now gets to
1306  // disk
1307  out.flush();
1308 
1309  AssertThrow(out, ExcIO());
1310 }
1311 
1312 
1313 
1314 template <int dim, int spacedim>
1315 void
1317  std::ostream &,
1318  const Mapping<dim, spacedim> *) const
1319 {
1320  Assert(false, ExcNotImplemented());
1321 }
1322 
1323 
1324 // TODO:[GK] Obey parameters
1325 template <>
1326 void
1328  std::ostream & out,
1329  const Mapping<2> * /*mapping*/) const
1330 {
1331  const int dim = 2;
1332  const int spacedim = 2;
1333 
1334  const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
1335 
1336  // The following text was copied
1337  // from an existing XFig file.
1338  out << "#FIG 3.2\nLandscape\nCenter\nInches" << std::endl
1339  << "A4\n100.00\nSingle"
1340  << std::endl
1341  // Background is transparent
1342  << "-3" << std::endl
1343  << "# generated by deal.II GridOut class" << std::endl
1344  << "# reduce first number to scale up image" << std::endl
1345  << "1200 2" << std::endl;
1346  // Write custom palette
1347  // grey
1348  unsigned int colno = 32;
1349  out << "0 " << colno++ << " #ff0000" << std::endl;
1350  out << "0 " << colno++ << " #ff8000" << std::endl;
1351  out << "0 " << colno++ << " #ffd000" << std::endl;
1352  out << "0 " << colno++ << " #ffff00" << std::endl;
1353  out << "0 " << colno++ << " #c0ff00" << std::endl;
1354  out << "0 " << colno++ << " #80ff00" << std::endl;
1355  out << "0 " << colno++ << " #00f000" << std::endl;
1356  out << "0 " << colno++ << " #00f0c0" << std::endl;
1357  out << "0 " << colno++ << " #00f0ff" << std::endl;
1358  out << "0 " << colno++ << " #00c0ff" << std::endl;
1359  out << "0 " << colno++ << " #0080ff" << std::endl;
1360  out << "0 " << colno++ << " #0040ff" << std::endl;
1361  out << "0 " << colno++ << " #0000c0" << std::endl;
1362  out << "0 " << colno++ << " #5000ff" << std::endl;
1363  out << "0 " << colno++ << " #8000ff" << std::endl;
1364  out << "0 " << colno++ << " #b000ff" << std::endl;
1365  out << "0 " << colno++ << " #ff00ff" << std::endl;
1366  out << "0 " << colno++ << " #ff80ff" << std::endl;
1367  // grey
1368  for (unsigned int i = 0; i < 8; ++i)
1369  out << "0 " << colno++ << " #" << std::hex << 32 * i + 31 << 32 * i + 31
1370  << 32 * i + 31 << std::dec << std::endl;
1371  // green
1372  for (unsigned int i = 1; i < 16; ++i)
1373  out << "0 " << colno++ << " #00" << std::hex << 16 * i + 15 << std::dec
1374  << "00" << std::endl;
1375  // yellow
1376  for (unsigned int i = 1; i < 16; ++i)
1377  out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << 16 * i + 15
1378  << std::dec << "00" << std::endl;
1379  // red
1380  for (unsigned int i = 1; i < 16; ++i)
1381  out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << std::dec
1382  << "0000" << std::endl;
1383  // purple
1384  for (unsigned int i = 1; i < 16; ++i)
1385  out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << "00"
1386  << 16 * i + 15 << std::dec << std::endl;
1387  // blue
1388  for (unsigned int i = 1; i < 16; ++i)
1389  out << "0 " << colno++ << " #0000" << std::hex << 16 * i + 15 << std::dec
1390  << std::endl;
1391  // cyan
1392  for (unsigned int i = 1; i < 16; ++i)
1393  out << "0 " << colno++ << " #00" << std::hex << 16 * i + 15 << 16 * i + 15
1394  << std::dec << std::endl;
1395 
1396  // We write all cells and cells on
1397  // coarser levels are behind cells
1398  // on finer levels. Level 0
1399  // corresponds to a depth of 900,
1400  // each level subtracting 1
1403 
1404  for (; cell != end; ++cell)
1405  {
1406  // If depth is not encoded, write finest level only
1407  if (!xfig_flags.level_depth && !cell->active())
1408  continue;
1409  // Code for polygon
1410  out << "2 3 " << xfig_flags.line_style << ' '
1412  // with black line
1413  << " 0 ";
1414  // Fill color
1415  switch (xfig_flags.color_by)
1416  {
1417  // TODO[GK]: Simplify after deprecation period is over
1419  out << static_cast<unsigned int>(cell->material_id()) + 32;
1420  break;
1422  out << cell->level() + 8;
1423  break;
1425  out << cell->subdomain_id() + 32;
1426  break;
1428  out << cell->level_subdomain_id() + 32;
1429  break;
1430  default:
1431  Assert(false, ExcInternalError());
1432  }
1433 
1434  // Depth, unused, fill
1435  out << ' '
1436  << (xfig_flags.level_depth ? (900 - cell->level()) :
1437  (900 + cell->material_id()))
1438  << " 0 " << xfig_flags.fill_style
1439  << " 0.0 "
1440  // some style parameters
1441  << " 0 0 -1 0 0 "
1442  // number of points
1443  << nv + 1 << std::endl;
1444 
1445  // For each point, write scaled
1446  // and shifted coordinates
1447  // multiplied by 1200
1448  // (dots/inch)
1449  for (unsigned int k = 0; k <= nv; ++k)
1450  {
1451  const Point<dim> &p =
1452  cell->vertex(GeometryInfo<dim>::ucd_to_deal[k % nv]);
1453  for (unsigned int d = 0; d < static_cast<unsigned int>(dim); ++d)
1454  {
1455  int val = static_cast<int>(1200 * xfig_flags.scaling(d) *
1456  (p(d) - xfig_flags.offset(d)));
1457  out << '\t' << ((d == 0) ? val : -val);
1458  }
1459  out << std::endl;
1460  }
1461  // Now write boundary edges
1462  static const unsigned int face_reorder[4] = {2, 1, 3, 0};
1464  for (const unsigned int f : face_reorder)
1465  {
1466  Triangulation<dim, spacedim>::face_iterator face = cell->face(f);
1467  const types::boundary_id bi = face->boundary_id();
1469  {
1470  // Code for polyline
1471  out << "2 1 "
1472  // with line style and thickness
1473  << xfig_flags.boundary_style << ' '
1474  << xfig_flags.boundary_thickness << ' ' << 1 + bi;
1475  // Fill color
1476  out << " -1 ";
1477  // Depth 100 less than cells
1478  out << (xfig_flags.level_depth ? (800 - cell->level()) :
1479  800 + bi)
1480  // unused, no fill
1481  << " 0 -1 0.0 "
1482  // some style parameters
1483  << " 0 0 -1 0 0 "
1484  // number of points
1485  << GeometryInfo<dim>::vertices_per_face << std::endl;
1486 
1487  // For each point, write scaled
1488  // and shifted coordinates
1489  // multiplied by 1200
1490  // (dots/inch)
1491 
1492  for (unsigned int k = 0;
1493  k < GeometryInfo<dim>::vertices_per_face;
1494  ++k)
1495  {
1496  const Point<dim> &p = face->vertex(k % nv);
1497  for (unsigned int d = 0; d < static_cast<unsigned int>(dim);
1498  ++d)
1499  {
1500  int val =
1501  static_cast<int>(1200 * xfig_flags.scaling(d) *
1502  (p(d) - xfig_flags.offset(d)));
1503  out << '\t' << ((d == 0) ? val : -val);
1504  }
1505  out << std::endl;
1506  }
1507  }
1508  }
1509  }
1510 
1511  // make sure everything now gets to
1512  // disk
1513  out.flush();
1514 
1515  AssertThrow(out, ExcIO());
1516 }
1517 
1518 
1519 
1520 namespace
1521 {
1532  Point<2>
1533  svg_project_point(const Point<3> & point,
1534  const Point<3> & camera_position,
1535  const Tensor<1, 3> &camera_direction,
1536  const Tensor<1, 3> &camera_horizontal,
1537  const float camera_focus)
1538  {
1539  const Tensor<1, 3> camera_vertical =
1540  cross_product_3d(camera_horizontal, camera_direction);
1541 
1542  const float phi =
1543  camera_focus / ((point - camera_position) * camera_direction);
1544 
1545  const Point<3> projection =
1546  camera_position + phi * (point - camera_position);
1547 
1548  return {(projection - camera_position - camera_focus * camera_direction) *
1549  camera_horizontal,
1550  (projection - camera_position - camera_focus * camera_direction) *
1551  camera_vertical};
1552  }
1553 } // namespace
1554 
1555 
1556 
1557 template <int dim, int spacedim>
1558 void
1560  std::ostream & /*out*/) const
1561 {
1562  Assert(false, ExcNotImplemented());
1563 }
1564 
1565 
1566 void
1567 GridOut::write_svg(const Triangulation<2, 2> &tria, std::ostream &out) const
1568 {
1569  unsigned int n = 0;
1570 
1571  unsigned int min_level, max_level;
1572 
1573  // Svg files require an underlying drawing grid. The size of this
1574  // grid is provided in the parameters height and width. Each of them
1575  // may be zero, such that it is computed from the other. Obviously,
1576  // both of them zero does not produce reasonable output.
1577  unsigned int height = svg_flags.height;
1578  unsigned int width = svg_flags.width;
1579  Assert(height != 0 || width != 0,
1580  ExcMessage("You have to set at least one of width and height"));
1581 
1582  unsigned int margin_in_percent = 0;
1584  margin_in_percent = 8;
1585 
1586  // initial font size for cell labels
1587  unsigned int cell_label_font_size;
1588 
1589  // get date and time
1590  // time_t time_stamp;
1591  // tm *now;
1592  // time_stamp = time(0);
1593  // now = localtime(&time_stamp);
1594 
1595  float camera_focus;
1596 
1597  Point<3> point;
1598  Point<2> projection_decomposition;
1599 
1600  float x_max_perspective, x_min_perspective;
1601  float y_max_perspective, y_min_perspective;
1602 
1603  float x_dimension_perspective, y_dimension_perspective;
1604 
1605 
1606  // auxiliary variables for the bounding box and the range of cell levels
1607  double x_min = tria.begin()->vertex(0)[0];
1608  double x_max = x_min;
1609  double y_min = tria.begin()->vertex(0)[1];
1610  double y_max = y_min;
1611 
1612  double x_dimension, y_dimension;
1613 
1614  min_level = max_level = tria.begin()->level();
1615 
1616  // auxiliary set for the materials being used
1617  std::set<unsigned int> materials;
1618 
1619  // auxiliary set for the levels being used
1620  std::set<unsigned int> levels;
1621 
1622  // auxiliary set for the subdomains being used
1623  std::set<unsigned int> subdomains;
1624 
1625  // auxiliary set for the level subdomains being used
1626  std::set<int> level_subdomains;
1627 
1628  // We use an active cell iterator to determine the
1629  // bounding box of the given triangulation and check
1630  // the cells for material id, level number, subdomain id
1631  // (, and level subdomain id).
1632  for (const auto &cell : tria.cell_iterators())
1633  {
1634  for (unsigned int vertex_index = 0; vertex_index < 4; vertex_index++)
1635  {
1636  if (cell->vertex(vertex_index)[0] < x_min)
1637  x_min = cell->vertex(vertex_index)[0];
1638  if (cell->vertex(vertex_index)[0] > x_max)
1639  x_max = cell->vertex(vertex_index)[0];
1640 
1641  if (cell->vertex(vertex_index)[1] < y_min)
1642  y_min = cell->vertex(vertex_index)[1];
1643  if (cell->vertex(vertex_index)[1] > y_max)
1644  y_max = cell->vertex(vertex_index)[1];
1645  }
1646 
1647  if (static_cast<unsigned int>(cell->level()) < min_level)
1648  min_level = cell->level();
1649  if (static_cast<unsigned int>(cell->level()) > max_level)
1650  max_level = cell->level();
1651 
1652  materials.insert(cell->material_id());
1653  levels.insert(cell->level());
1654  if (cell->active())
1655  subdomains.insert(cell->subdomain_id() + 2);
1656  level_subdomains.insert(cell->level_subdomain_id() + 2);
1657  }
1658 
1659  x_dimension = x_max - x_min;
1660  y_dimension = y_max - y_min;
1661 
1662  // count the materials being used
1663  const unsigned int n_materials = materials.size();
1664 
1665  // count the levels being used
1666  const unsigned int n_levels = levels.size();
1667 
1668  // count the subdomains being used
1669  const unsigned int n_subdomains = subdomains.size();
1670 
1671  // count the level subdomains being used
1672  const unsigned int n_level_subdomains = level_subdomains.size();
1673 
1674  switch (svg_flags.coloring)
1675  {
1677  n = n_materials;
1678  break;
1680  n = n_levels;
1681  break;
1683  n = n_subdomains;
1684  break;
1686  n = n_level_subdomains;
1687  break;
1688  default:
1689  break;
1690  }
1691 
1692  // set the camera position to top view, targeting at the origin
1693  // vectors and variables for the perspective view
1694  Point<3> camera_position;
1695  camera_position[0] = 0;
1696  camera_position[1] = 0;
1697  camera_position[2] = 2. * std::max(x_dimension, y_dimension);
1698 
1699  Tensor<1, 3> camera_direction;
1700  camera_direction[0] = 0;
1701  camera_direction[1] = 0;
1702  camera_direction[2] = -1;
1703 
1704  Tensor<1, 3> camera_horizontal;
1705  camera_horizontal[0] = 1;
1706  camera_horizontal[1] = 0;
1707  camera_horizontal[2] = 0;
1708 
1709  camera_focus = .5 * std::max(x_dimension, y_dimension);
1710 
1711  Point<3> camera_position_temp;
1712  Point<3> camera_direction_temp;
1713  Point<3> camera_horizontal_temp;
1714 
1715  const double angle_factor = 3.14159265 / 180.;
1716 
1717  // (I) rotate the camera to the chosen polar angle
1718  camera_position_temp[1] =
1719  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[1] -
1720  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[2];
1721  camera_position_temp[2] =
1722  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[1] +
1723  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[2];
1724 
1725  camera_direction_temp[1] =
1726  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[1] -
1727  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1728  camera_direction_temp[2] =
1729  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[1] +
1730  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1731 
1732  camera_horizontal_temp[1] =
1733  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] -
1734  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1735  camera_horizontal_temp[2] =
1736  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] +
1737  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1738 
1739  camera_position[1] = camera_position_temp[1];
1740  camera_position[2] = camera_position_temp[2];
1741 
1742  camera_direction[1] = camera_direction_temp[1];
1743  camera_direction[2] = camera_direction_temp[2];
1744 
1745  camera_horizontal[1] = camera_horizontal_temp[1];
1746  camera_horizontal[2] = camera_horizontal_temp[2];
1747 
1748  // (II) rotate the camera to the chosen azimuth angle
1749  camera_position_temp[0] =
1750  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[0] -
1751  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1752  camera_position_temp[1] =
1753  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[0] +
1754  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1755 
1756  camera_direction_temp[0] =
1757  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] -
1758  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1759  camera_direction_temp[1] =
1760  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] +
1761  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1762 
1763  camera_horizontal_temp[0] =
1764  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] -
1765  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1766  camera_horizontal_temp[1] =
1767  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] +
1768  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1769 
1770  camera_position[0] = camera_position_temp[0];
1771  camera_position[1] = camera_position_temp[1];
1772 
1773  camera_direction[0] = camera_direction_temp[0];
1774  camera_direction[1] = camera_direction_temp[1];
1775 
1776  camera_horizontal[0] = camera_horizontal_temp[0];
1777  camera_horizontal[1] = camera_horizontal_temp[1];
1778 
1779  // translate the camera to the given triangulation
1780  camera_position[0] = x_min + .5 * x_dimension;
1781  camera_position[1] = y_min + .5 * y_dimension;
1782 
1783  camera_position[0] += 2. * std::max(x_dimension, y_dimension) *
1784  std::sin(angle_factor * svg_flags.polar_angle) *
1785  std::sin(angle_factor * svg_flags.azimuth_angle);
1786  camera_position[1] -= 2. * std::max(x_dimension, y_dimension) *
1787  std::sin(angle_factor * svg_flags.polar_angle) *
1788  std::cos(angle_factor * svg_flags.azimuth_angle);
1789 
1790 
1791  // determine the bounding box of the given triangulation on the projection
1792  // plane of the camera viewing system
1793  point[0] = tria.begin()->vertex(0)[0];
1794  point[1] = tria.begin()->vertex(0)[1];
1795  point[2] = 0;
1796 
1797  float min_level_min_vertex_distance = 0;
1798 
1800  {
1801  point[2] = svg_flags.level_height_factor *
1802  (static_cast<float>(tria.begin()->level()) /
1803  static_cast<float>(n_levels)) *
1804  std::max(x_dimension, y_dimension);
1805  }
1806 
1807  projection_decomposition = svg_project_point(
1808  point, camera_position, camera_direction, camera_horizontal, camera_focus);
1809 
1810  x_max_perspective = projection_decomposition[0];
1811  x_min_perspective = projection_decomposition[0];
1812 
1813  y_max_perspective = projection_decomposition[1];
1814  y_min_perspective = projection_decomposition[1];
1815 
1816  for (Triangulation<2, 2>::cell_iterator cell = tria.begin();
1817  cell != tria.end();
1818  ++cell)
1819  {
1820  point[0] = cell->vertex(0)[0];
1821  point[1] = cell->vertex(0)[1];
1822  point[2] = 0;
1823 
1825  {
1826  point[2] =
1828  (static_cast<float>(cell->level()) / static_cast<float>(n_levels)) *
1829  std::max(x_dimension, y_dimension);
1830  }
1831 
1832  projection_decomposition = svg_project_point(point,
1833  camera_position,
1834  camera_direction,
1835  camera_horizontal,
1836  camera_focus);
1837 
1838  if (x_max_perspective < projection_decomposition[0])
1839  x_max_perspective = projection_decomposition[0];
1840  if (x_min_perspective > projection_decomposition[0])
1841  x_min_perspective = projection_decomposition[0];
1842 
1843  if (y_max_perspective < projection_decomposition[1])
1844  y_max_perspective = projection_decomposition[1];
1845  if (y_min_perspective > projection_decomposition[1])
1846  y_min_perspective = projection_decomposition[1];
1847 
1848  point[0] = cell->vertex(1)[0];
1849  point[1] = cell->vertex(1)[1];
1850 
1851  projection_decomposition = svg_project_point(point,
1852  camera_position,
1853  camera_direction,
1854  camera_horizontal,
1855  camera_focus);
1856 
1857  if (x_max_perspective < projection_decomposition[0])
1858  x_max_perspective = projection_decomposition[0];
1859  if (x_min_perspective > projection_decomposition[0])
1860  x_min_perspective = projection_decomposition[0];
1861 
1862  if (y_max_perspective < projection_decomposition[1])
1863  y_max_perspective = projection_decomposition[1];
1864  if (y_min_perspective > projection_decomposition[1])
1865  y_min_perspective = projection_decomposition[1];
1866 
1867  point[0] = cell->vertex(2)[0];
1868  point[1] = cell->vertex(2)[1];
1869 
1870  projection_decomposition = svg_project_point(point,
1871  camera_position,
1872  camera_direction,
1873  camera_horizontal,
1874  camera_focus);
1875 
1876  if (x_max_perspective < projection_decomposition[0])
1877  x_max_perspective = projection_decomposition[0];
1878  if (x_min_perspective > projection_decomposition[0])
1879  x_min_perspective = projection_decomposition[0];
1880 
1881  if (y_max_perspective < projection_decomposition[1])
1882  y_max_perspective = projection_decomposition[1];
1883  if (y_min_perspective > projection_decomposition[1])
1884  y_min_perspective = projection_decomposition[1];
1885 
1886  point[0] = cell->vertex(3)[0];
1887  point[1] = cell->vertex(3)[1];
1888 
1889  projection_decomposition = svg_project_point(point,
1890  camera_position,
1891  camera_direction,
1892  camera_horizontal,
1893  camera_focus);
1894 
1895  if (x_max_perspective < projection_decomposition[0])
1896  x_max_perspective = projection_decomposition[0];
1897  if (x_min_perspective > projection_decomposition[0])
1898  x_min_perspective = projection_decomposition[0];
1899 
1900  if (y_max_perspective < projection_decomposition[1])
1901  y_max_perspective = projection_decomposition[1];
1902  if (y_min_perspective > projection_decomposition[1])
1903  y_min_perspective = projection_decomposition[1];
1904 
1905  if (static_cast<unsigned int>(cell->level()) == min_level)
1906  min_level_min_vertex_distance = cell->minimum_vertex_distance();
1907  }
1908 
1909  x_dimension_perspective = x_max_perspective - x_min_perspective;
1910  y_dimension_perspective = y_max_perspective - y_min_perspective;
1911 
1912  // create the svg file with an internal style sheet
1913  if (width == 0)
1914  width = static_cast<unsigned int>(
1915  .5 + height * (x_dimension_perspective / y_dimension_perspective));
1916  else if (height == 0)
1917  height = static_cast<unsigned int>(
1918  .5 + width * (y_dimension_perspective / x_dimension_perspective));
1919  unsigned int additional_width = 0;
1920  // font size for date, time, legend, and colorbar
1921  unsigned int font_size =
1922  static_cast<unsigned int>(.5 + (height / 100.) * 1.75);
1923  cell_label_font_size = static_cast<unsigned int>(
1924  .5 + (height * .15 * svg_flags.cell_font_scaling *
1925  min_level_min_vertex_distance / std::min(x_dimension, y_dimension)));
1926 
1927  if (svg_flags.draw_legend &&
1931  {
1932  additional_width = static_cast<unsigned int>(
1933  .5 + height * .4); // additional width for legend
1934  }
1935  else if (svg_flags.draw_colorbar && svg_flags.coloring)
1936  {
1937  additional_width = static_cast<unsigned int>(
1938  .5 + height * .175); // additional width for colorbar
1939  }
1940 
1941  // out << "<!-- deal.ii GridOut " << now->tm_mday << '/' << now->tm_mon + 1 <<
1942  // '/' << now->tm_year + 1900
1943  // << ' ' << now->tm_hour << ':';
1944  //
1945  // if (now->tm_min < 10) out << '0';
1946  //
1947  // out << now->tm_min << " -->" << '\n';
1948 
1949  // basic svg header
1950  out << "<svg width=\"" << width + additional_width << "\" height=\"" << height
1951  << "\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" << '\n'
1952  << '\n';
1953 
1954 
1956  {
1957  out
1958  << " <linearGradient id=\"background_gradient\" gradientUnits=\"userSpaceOnUse\" x1=\"0\" y1=\"0\" x2=\"0\" y2=\""
1959  << height << "\">" << '\n'
1960  << " <stop offset=\"0\" style=\"stop-color:white\"/>" << '\n'
1961  << " <stop offset=\"1\" style=\"stop-color:lightsteelblue\"/>" << '\n'
1962  << " </linearGradient>" << '\n';
1963  }
1964 
1965  out << '\n';
1966 
1967  // header for the internal style sheet
1968  out << "<!-- internal style sheet -->" << '\n'
1969  << "<style type=\"text/css\"><![CDATA[" << '\n';
1970 
1971  // set the background of the output graphic
1973  out << " rect.background{fill:url(#background_gradient)}" << '\n';
1975  out << " rect.background{fill:white}" << '\n';
1976  else
1977  out << " rect.background{fill:none}" << '\n';
1978 
1979  // basic svg graphic element styles
1980  out << " rect{fill:none; stroke:rgb(25,25,25); stroke-width:"
1981  << svg_flags.line_thickness << '}' << '\n'
1982  << " text{font-family:Helvetica; text-anchor:middle; fill:rgb(25,25,25)}"
1983  << '\n'
1984  << " line{stroke:rgb(25,25,25); stroke-width:"
1985  << svg_flags.boundary_line_thickness << '}' << '\n'
1986  << " path{fill:none; stroke:rgb(25,25,25); stroke-width:"
1987  << svg_flags.line_thickness << '}' << '\n'
1988  << '\n';
1989 
1990  // polygon styles with respect to the chosen cell coloring
1991  if (svg_flags.coloring)
1992  {
1993  unsigned int labeling_index = 0;
1994  auto materials_it = materials.begin();
1995  auto levels_it = levels.begin();
1996  auto subdomains_it = subdomains.begin();
1997  auto level_subdomains_it = level_subdomains.begin();
1998 
1999  for (unsigned int index = 0; index < n; index++)
2000  {
2001  double h;
2002 
2003  if (n != 1)
2004  h = .6 - (index / (n - 1.)) * .6;
2005  else
2006  h = .6;
2007 
2008  unsigned int r = 0;
2009  unsigned int g = 0;
2010  unsigned int b = 0;
2011 
2012  unsigned int i = static_cast<unsigned int>(h * 6);
2013 
2014  double f = h * 6 - i;
2015  double q = 1 - f;
2016  double t = f;
2017 
2018  switch (i % 6)
2019  {
2020  case 0:
2021  r = 255, g = static_cast<unsigned int>(.5 + 255 * t);
2022  break;
2023  case 1:
2024  r = static_cast<unsigned int>(.5 + 255 * q), g = 255;
2025  break;
2026  case 2:
2027  g = 255, b = static_cast<unsigned int>(.5 + 255 * t);
2028  break;
2029  case 3:
2030  g = static_cast<unsigned int>(.5 + 255 * q), b = 255;
2031  break;
2032  case 4:
2033  r = static_cast<unsigned int>(.5 + 255 * t), b = 255;
2034  break;
2035  case 5:
2036  r = 255, b = static_cast<unsigned int>(.5 + 255 * q);
2037  break;
2038  default:
2039  break;
2040  }
2041 
2042  switch (svg_flags.coloring)
2043  {
2045  labeling_index = *materials_it++;
2046  break;
2048  labeling_index = *levels_it++;
2049  break;
2051  labeling_index = *subdomains_it++;
2052  break;
2054  labeling_index = *level_subdomains_it++;
2055  break;
2056  default:
2057  break;
2058  }
2059 
2060  out << " path.p" << labeling_index << "{fill:rgb(" << r << ',' << g
2061  << ',' << b << "); "
2062  << "stroke:rgb(25,25,25); stroke-width:"
2063  << svg_flags.line_thickness << '}' << '\n';
2064 
2065  out << " path.ps" << labeling_index << "{fill:rgb("
2066  << static_cast<unsigned int>(.5 + .75 * r) << ','
2067  << static_cast<unsigned int>(.5 + .75 * g) << ','
2068  << static_cast<unsigned int>(.5 + .75 * b) << "); "
2069  << "stroke:rgb(20,20,20); stroke-width:"
2070  << svg_flags.line_thickness << '}' << '\n';
2071 
2072  out << " rect.r" << labeling_index << "{fill:rgb(" << r << ',' << g
2073  << ',' << b << "); "
2074  << "stroke:rgb(25,25,25); stroke-width:"
2075  << svg_flags.line_thickness << '}' << '\n';
2076 
2077  labeling_index++;
2078  }
2079  }
2080 
2081  out << "]]></style>" << '\n' << '\n';
2082 
2083  // background rectangle
2084  out << " <rect class=\"background\" width=\"" << width << "\" height=\""
2085  << height << "\"/>" << '\n';
2086 
2088  {
2089  unsigned int x_offset = 0;
2090 
2091  if (svg_flags.margin)
2092  x_offset = static_cast<unsigned int>(.5 + (height / 100.) *
2093  (margin_in_percent / 2.));
2094  else
2095  x_offset = static_cast<unsigned int>(.5 + height * .025);
2096 
2097  out
2098  << " <text x=\"" << x_offset << "\" y=\""
2099  << static_cast<unsigned int>(.5 + height * .0525) << '\"'
2100  << " style=\"font-weight:100; fill:lightsteelblue; text-anchor:start; font-family:Courier; font-size:"
2101  << static_cast<unsigned int>(.5 + height * .045) << "px\">"
2102  << "deal.II"
2103  << "</text>" << '\n';
2104 
2105  // out << " <text x=\"" << x_offset + static_cast<unsigned int>(.5 +
2106  // height * .045 * 4.75) << "\" y=\"" << static_cast<unsigned int>(.5 +
2107  // height * .0525) << '\"'
2108  // << " style=\"fill:lightsteelblue; text-anchor:start; font-size:" <<
2109  // font_size << "\">"
2110  // << now->tm_mday << '/' << now->tm_mon + 1 << '/' << now->tm_year +
2111  // 1900
2112  // << " - " << now->tm_hour << ':';
2113  //
2114  // if(now->tm_min < 10) out << '0';
2115  //
2116  // out << now->tm_min
2117  // << "</text>"<< '\n' << '\n';
2118  }
2119 
2120  // draw the cells, starting out from the minimal level (in order to guaranty a
2121  // correct perspective view)
2122  out << " <!-- cells -->" << '\n';
2123 
2124  for (unsigned int level_index = min_level; level_index <= max_level;
2125  level_index++)
2126  {
2127  Triangulation<2, 2>::cell_iterator cell = tria.begin(level_index),
2128  endc = tria.end(level_index);
2129 
2130  for (; cell != endc; ++cell)
2131  {
2132  if (!svg_flags.convert_level_number_to_height && !cell->active())
2133  continue;
2134 
2135  // draw the current cell
2136  out << " <path";
2137 
2138  if (svg_flags.coloring)
2139  {
2140  out << " class=\"p";
2141 
2142  if (!cell->active() && svg_flags.convert_level_number_to_height)
2143  out << 's';
2144 
2145  switch (svg_flags.coloring)
2146  {
2148  out << cell->material_id();
2149  break;
2151  out << static_cast<unsigned int>(cell->level());
2152  break;
2154  if (cell->active())
2155  out << cell->subdomain_id() + 2;
2156  else
2157  out << 'X';
2158  break;
2160  out << cell->level_subdomain_id() + 2;
2161  break;
2162  default:
2163  break;
2164  }
2165 
2166  out << '\"';
2167  }
2168 
2169  out << " d=\"M ";
2170 
2171  point[0] = cell->vertex(0)[0];
2172  point[1] = cell->vertex(0)[1];
2173  point[2] = 0;
2174 
2176  {
2177  point[2] = svg_flags.level_height_factor *
2178  (static_cast<float>(cell->level()) /
2179  static_cast<float>(n_levels)) *
2180  std::max(x_dimension, y_dimension);
2181  }
2182 
2183  projection_decomposition = svg_project_point(point,
2184  camera_position,
2185  camera_direction,
2186  camera_horizontal,
2187  camera_focus);
2188 
2189  out << static_cast<unsigned int>(
2190  .5 +
2191  ((projection_decomposition[0] - x_min_perspective) /
2192  x_dimension_perspective) *
2193  (width - (width / 100.) * 2. * margin_in_percent) +
2194  ((width / 100.) * margin_in_percent))
2195  << ' '
2196  << static_cast<unsigned int>(
2197  .5 + height - (height / 100.) * margin_in_percent -
2198  ((projection_decomposition[1] - y_min_perspective) /
2199  y_dimension_perspective) *
2200  (height - (height / 100.) * 2. * margin_in_percent));
2201 
2202  out << " L ";
2203 
2204  point[0] = cell->vertex(1)[0];
2205  point[1] = cell->vertex(1)[1];
2206 
2207  projection_decomposition = svg_project_point(point,
2208  camera_position,
2209  camera_direction,
2210  camera_horizontal,
2211  camera_focus);
2212 
2213  out << static_cast<unsigned int>(
2214  .5 +
2215  ((projection_decomposition[0] - x_min_perspective) /
2216  x_dimension_perspective) *
2217  (width - (width / 100.) * 2. * margin_in_percent) +
2218  ((width / 100.) * margin_in_percent))
2219  << ' '
2220  << static_cast<unsigned int>(
2221  .5 + height - (height / 100.) * margin_in_percent -
2222  ((projection_decomposition[1] - y_min_perspective) /
2223  y_dimension_perspective) *
2224  (height - (height / 100.) * 2. * margin_in_percent));
2225 
2226  out << " L ";
2227 
2228  point[0] = cell->vertex(3)[0];
2229  point[1] = cell->vertex(3)[1];
2230 
2231  projection_decomposition = svg_project_point(point,
2232  camera_position,
2233  camera_direction,
2234  camera_horizontal,
2235  camera_focus);
2236 
2237  out << static_cast<unsigned int>(
2238  .5 +
2239  ((projection_decomposition[0] - x_min_perspective) /
2240  x_dimension_perspective) *
2241  (width - (width / 100.) * 2. * margin_in_percent) +
2242  ((width / 100.) * margin_in_percent))
2243  << ' '
2244  << static_cast<unsigned int>(
2245  .5 + height - (height / 100.) * margin_in_percent -
2246  ((projection_decomposition[1] - y_min_perspective) /
2247  y_dimension_perspective) *
2248  (height - (height / 100.) * 2. * margin_in_percent));
2249 
2250  out << " L ";
2251 
2252  point[0] = cell->vertex(2)[0];
2253  point[1] = cell->vertex(2)[1];
2254 
2255  projection_decomposition = svg_project_point(point,
2256  camera_position,
2257  camera_direction,
2258  camera_horizontal,
2259  camera_focus);
2260 
2261  out << static_cast<unsigned int>(
2262  .5 +
2263  ((projection_decomposition[0] - x_min_perspective) /
2264  x_dimension_perspective) *
2265  (width - (width / 100.) * 2. * margin_in_percent) +
2266  ((width / 100.) * margin_in_percent))
2267  << ' '
2268  << static_cast<unsigned int>(
2269  .5 + height - (height / 100.) * margin_in_percent -
2270  ((projection_decomposition[1] - y_min_perspective) /
2271  y_dimension_perspective) *
2272  (height - (height / 100.) * 2. * margin_in_percent));
2273 
2274  out << " L ";
2275 
2276  point[0] = cell->vertex(0)[0];
2277  point[1] = cell->vertex(0)[1];
2278 
2279  projection_decomposition = svg_project_point(point,
2280  camera_position,
2281  camera_direction,
2282  camera_horizontal,
2283  camera_focus);
2284 
2285  out << static_cast<unsigned int>(
2286  .5 +
2287  ((projection_decomposition[0] - x_min_perspective) /
2288  x_dimension_perspective) *
2289  (width - (width / 100.) * 2. * margin_in_percent) +
2290  ((width / 100.) * margin_in_percent))
2291  << ' '
2292  << static_cast<unsigned int>(
2293  .5 + height - (height / 100.) * margin_in_percent -
2294  ((projection_decomposition[1] - y_min_perspective) /
2295  y_dimension_perspective) *
2296  (height - (height / 100.) * 2. * margin_in_percent));
2297 
2298  out << "\"/>" << '\n';
2299 
2300  // label the current cell
2304  {
2305  point[0] = cell->center()[0];
2306  point[1] = cell->center()[1];
2307  point[2] = 0;
2308 
2310  {
2311  point[2] = svg_flags.level_height_factor *
2312  (static_cast<float>(cell->level()) /
2313  static_cast<float>(n_levels)) *
2314  std::max(x_dimension, y_dimension);
2315  }
2316 
2317  float distance_to_camera =
2318  std::sqrt(std::pow(point[0] - camera_position[0], 2.) +
2319  std::pow(point[1] - camera_position[1], 2.) +
2320  std::pow(point[2] - camera_position[2], 2.));
2321  float distance_factor =
2322  distance_to_camera / (2. * std::max(x_dimension, y_dimension));
2323 
2324  projection_decomposition = svg_project_point(point,
2325  camera_position,
2326  camera_direction,
2327  camera_horizontal,
2328  camera_focus);
2329 
2330  const auto font_size_this_cell = static_cast<unsigned int>(
2331  .5 +
2332  cell_label_font_size *
2333  std::pow(.5, cell->level() - 4. + 3.5 * distance_factor));
2334 
2335  out << " <text"
2336  << " x=\""
2337  << static_cast<unsigned int>(
2338  .5 +
2339  ((projection_decomposition[0] - x_min_perspective) /
2340  x_dimension_perspective) *
2341  (width - (width / 100.) * 2. * margin_in_percent) +
2342  ((width / 100.) * margin_in_percent))
2343  << "\" y=\""
2344  << static_cast<unsigned int>(
2345  .5 + height - (height / 100.) * margin_in_percent -
2346  ((projection_decomposition[1] - y_min_perspective) /
2347  y_dimension_perspective) *
2348  (height - (height / 100.) * 2. * margin_in_percent) +
2349  0.5 * font_size_this_cell)
2350  << "\" style=\"font-size:" << font_size_this_cell << "px\">";
2351 
2353  {
2354  out << cell->level();
2355  }
2356 
2358  {
2360  out << '.';
2361  out << cell->index();
2362  }
2363 
2365  {
2368  out << ',';
2369  out
2370  << static_cast<std::make_signed<types::material_id>::type>(
2371  cell->material_id());
2372  }
2373 
2375  {
2378  out << ',';
2379  if (cell->active())
2380  out << static_cast<
2381  std::make_signed<types::subdomain_id>::type>(
2382  cell->subdomain_id());
2383  else
2384  out << 'X';
2385  }
2386 
2388  {
2393  out << ',';
2394  out
2395  << static_cast<std::make_signed<types::subdomain_id>::type>(
2396  cell->level_subdomain_id());
2397  }
2398 
2399  out << "</text>" << '\n';
2400  }
2401 
2402  // if the current cell lies at the boundary of the triangulation, draw
2403  // the additional boundary line
2405  {
2406  for (unsigned int faceIndex = 0; faceIndex < 4; faceIndex++)
2407  {
2408  if (cell->at_boundary(faceIndex))
2409  {
2410  point[0] = cell->face(faceIndex)->vertex(0)[0];
2411  point[1] = cell->face(faceIndex)->vertex(0)[1];
2412  point[2] = 0;
2413 
2415  {
2416  point[2] = svg_flags.level_height_factor *
2417  (static_cast<float>(cell->level()) /
2418  static_cast<float>(n_levels)) *
2419  std::max(x_dimension, y_dimension);
2420  }
2421 
2422  projection_decomposition =
2423  svg_project_point(point,
2424  camera_position,
2425  camera_direction,
2426  camera_horizontal,
2427  camera_focus);
2428 
2429  out << " <line x1=\""
2430  << static_cast<unsigned int>(
2431  .5 +
2432  ((projection_decomposition[0] -
2433  x_min_perspective) /
2434  x_dimension_perspective) *
2435  (width -
2436  (width / 100.) * 2. * margin_in_percent) +
2437  ((width / 100.) * margin_in_percent))
2438  << "\" y1=\""
2439  << static_cast<unsigned int>(
2440  .5 + height -
2441  (height / 100.) * margin_in_percent -
2442  ((projection_decomposition[1] -
2443  y_min_perspective) /
2444  y_dimension_perspective) *
2445  (height -
2446  (height / 100.) * 2. * margin_in_percent));
2447 
2448  point[0] = cell->face(faceIndex)->vertex(1)[0];
2449  point[1] = cell->face(faceIndex)->vertex(1)[1];
2450  point[2] = 0;
2451 
2453  {
2454  point[2] = svg_flags.level_height_factor *
2455  (static_cast<float>(cell->level()) /
2456  static_cast<float>(n_levels)) *
2457  std::max(x_dimension, y_dimension);
2458  }
2459 
2460  projection_decomposition =
2461  svg_project_point(point,
2462  camera_position,
2463  camera_direction,
2464  camera_horizontal,
2465  camera_focus);
2466 
2467  out << "\" x2=\""
2468  << static_cast<unsigned int>(
2469  .5 +
2470  ((projection_decomposition[0] -
2471  x_min_perspective) /
2472  x_dimension_perspective) *
2473  (width -
2474  (width / 100.) * 2. * margin_in_percent) +
2475  ((width / 100.) * margin_in_percent))
2476  << "\" y2=\""
2477  << static_cast<unsigned int>(
2478  .5 + height -
2479  (height / 100.) * margin_in_percent -
2480  ((projection_decomposition[1] -
2481  y_min_perspective) /
2482  y_dimension_perspective) *
2483  (height -
2484  (height / 100.) * 2. * margin_in_percent))
2485  << "\"/>" << '\n';
2486  }
2487  }
2488  }
2489  }
2490  }
2491 
2492 
2493  // draw the legend
2494  if (svg_flags.draw_legend)
2495  out << '\n' << " <!-- legend -->" << '\n';
2496 
2497  additional_width = 0;
2498  if (!svg_flags.margin)
2499  additional_width = static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2500 
2501  // explanation of the cell labeling
2502  if (svg_flags.draw_legend &&
2506  {
2507  unsigned int line_offset = 0;
2508  out << " <rect x=\"" << width + additional_width << "\" y=\""
2509  << static_cast<unsigned int>(.5 + (height / 100.) * margin_in_percent)
2510  << "\" width=\""
2511  << static_cast<unsigned int>(.5 + (height / 100.) *
2512  (40. - margin_in_percent))
2513  << "\" height=\"" << static_cast<unsigned int>(.5 + height * .165)
2514  << "\"/>" << '\n';
2515 
2516  out << " <text x=\""
2517  << width + additional_width +
2518  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2519  << "\" y=\""
2520  << static_cast<unsigned int>(.5 +
2521  (height / 100.) * margin_in_percent +
2522  (++line_offset) * 1.5 * font_size)
2523  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2524  << font_size << "px\">"
2525  << "cell label"
2526  << "</text>" << '\n';
2527 
2529  {
2530  out << " <text x=\""
2531  << width + additional_width +
2532  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2533  << "\" y=\""
2534  << static_cast<unsigned int>(.5 +
2535  (height / 100.) * margin_in_percent +
2536  (++line_offset) * 1.5 * font_size)
2537  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2538  << font_size << "px\">"
2539  << "cell_level";
2540 
2544  out << '.';
2545 
2546  out << "</text>" << '\n';
2547  }
2548 
2550  {
2551  out << " <text x=\""
2552  << width + additional_width +
2553  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2554  << "\" y=\""
2555  << static_cast<unsigned int>(.5 +
2556  (height / 100.) * margin_in_percent +
2557  (++line_offset) * 1.5 * font_size)
2558  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2559  << font_size << "px\">"
2560  << "cell_index";
2561 
2564  out << ',';
2565 
2566  out << "</text>" << '\n';
2567  }
2568 
2570  {
2571  out << " <text x=\""
2572  << width + additional_width +
2573  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2574  << "\" y=\""
2575  << static_cast<unsigned int>(.5 +
2576  (height / 100.) * margin_in_percent +
2577  (++line_offset) * 1.5 * font_size)
2578  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2579  << font_size << "px\">"
2580  << "material_id";
2581 
2584  out << ',';
2585 
2586  out << "</text>" << '\n';
2587  }
2588 
2590  {
2591  out << " <text x= \""
2592  << width + additional_width +
2593  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2594  << "\" y=\""
2595  << static_cast<unsigned int>(.5 +
2596  (height / 100.) * margin_in_percent +
2597  (++line_offset) * 1.5 * font_size)
2598  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2599  << font_size << "px\">"
2600  << "subdomain_id";
2601 
2603  out << ',';
2604 
2605  out << "</text>" << '\n';
2606  }
2607 
2609  {
2610  out << " <text x= \""
2611  << width + additional_width +
2612  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2613  << "\" y=\""
2614  << static_cast<unsigned int>(.5 +
2615  (height / 100.) * margin_in_percent +
2616  (++line_offset) * 1.5 * font_size)
2617  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2618  << font_size << "px\">"
2619  << "level_subdomain_id"
2620  << "</text>" << '\n';
2621  }
2622  }
2623 
2624  // show azimuth angle and polar angle as text below the explanation of the
2625  // cell labeling
2626  if (svg_flags.draw_legend)
2627  {
2628  out << " <text x=\"" << width + additional_width << "\" y=\""
2629  << static_cast<unsigned int>(
2630  .5 + (height / 100.) * margin_in_percent + 10.75 * font_size)
2631  << "\" style=\"text-anchor:start; font-size:" << font_size << "px\">"
2632  << "azimuth: " << svg_flags.azimuth_angle
2633  << "°, polar: " << svg_flags.polar_angle << "°</text>" << '\n';
2634  }
2635 
2636 
2637  // draw the colorbar
2638  if (svg_flags.draw_colorbar && svg_flags.coloring)
2639  {
2640  out << '\n' << " <!-- colorbar -->" << '\n';
2641 
2642  out << " <text x=\"" << width + additional_width << "\" y=\""
2643  << static_cast<unsigned int>(
2644  .5 + (height / 100.) * (margin_in_percent + 29.) -
2645  (font_size / 1.25))
2646  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2647  << font_size << "px\">";
2648 
2649  switch (svg_flags.coloring)
2650  {
2651  case 1:
2652  out << "material_id";
2653  break;
2654  case 2:
2655  out << "level_number";
2656  break;
2657  case 3:
2658  out << "subdomain_id";
2659  break;
2660  case 4:
2661  out << "level_subdomain_id";
2662  break;
2663  default:
2664  break;
2665  }
2666 
2667  out << "</text>" << '\n';
2668 
2669  unsigned int element_height = static_cast<unsigned int>(
2670  ((height / 100.) * (71. - 2. * margin_in_percent)) / n);
2671  unsigned int element_width =
2672  static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2673 
2674  int labeling_index = 0;
2675  auto materials_it = materials.begin();
2676  auto levels_it = levels.begin();
2677  auto subdomains_it = subdomains.begin();
2678  auto level_subdomains_it = level_subdomains.begin();
2679 
2680  for (unsigned int index = 0; index < n; index++)
2681  {
2682  switch (svg_flags.coloring)
2683  {
2685  labeling_index = *materials_it++;
2686  break;
2688  labeling_index = *levels_it++;
2689  break;
2691  labeling_index = *subdomains_it++;
2692  break;
2694  labeling_index = *level_subdomains_it++;
2695  break;
2696  default:
2697  break;
2698  }
2699 
2700  out << " <rect class=\"r" << labeling_index << "\" x=\""
2701  << width + additional_width << "\" y=\""
2702  << static_cast<unsigned int>(.5 + (height / 100.) *
2703  (margin_in_percent + 29)) +
2704  (n - index - 1) * element_height
2705  << "\" width=\"" << element_width << "\" height=\""
2706  << element_height << "\"/>" << '\n';
2707 
2708  out << " <text x=\""
2709  << width + additional_width + 1.5 * element_width << "\" y=\""
2710  << static_cast<unsigned int>(.5 + (height / 100.) *
2711  (margin_in_percent + 29)) +
2712  (n - index - 1 + .5) * element_height +
2713  static_cast<unsigned int>(.5 + font_size * .35)
2714  << "\""
2715  << " style=\"text-anchor:start; font-size:"
2716  << static_cast<unsigned int>(.5 + font_size) << "px";
2717 
2718  if (index == 0 || index == n - 1)
2719  out << "; font-weight:bold";
2720 
2721  out << "\">" << labeling_index;
2722 
2723  if (index == n - 1)
2724  out << " max";
2725  if (index == 0)
2726  out << " min";
2727 
2728  out << "</text>" << '\n';
2729 
2730  labeling_index++;
2731  }
2732  }
2733 
2734 
2735  // finalize the svg file
2736  out << '\n' << "</svg>";
2737  out.flush();
2738 }
2739 
2740 
2741 template <>
2742 void
2743 GridOut::write_mathgl(const Triangulation<1> &, std::ostream &) const
2744 {
2745  // 1d specialization not done yet
2746  Assert(false, ExcNotImplemented());
2747 }
2748 
2749 
2750 template <int dim, int spacedim>
2751 void
2753  std::ostream & out) const
2754 {
2755  AssertThrow(out, ExcIO());
2756 
2757  // (i) write header
2758  {
2759  // block this to have local variables destroyed after use
2760  const std::time_t time1 = std::time(nullptr);
2761  const std::tm * time = std::localtime(&time1);
2762 
2763  out
2764  << "\n#"
2765  << "\n# This file was generated by the deal.II library."
2766  << "\n# Date = " << time->tm_year + 1900 << "/" << std::setfill('0')
2767  << std::setw(2) << time->tm_mon + 1 << "/" << std::setfill('0')
2768  << std::setw(2) << time->tm_mday << "\n# Time = " << std::setfill('0')
2769  << std::setw(2) << time->tm_hour << ":" << std::setfill('0')
2770  << std::setw(2) << time->tm_min << ":" << std::setfill('0')
2771  << std::setw(2) << time->tm_sec << "\n#"
2772  << "\n# For a description of the MathGL script format see the MathGL manual. "
2773  << "\n#"
2774  << "\n# Note: This file is understood by MathGL v2.1 and higher only, and can "
2775  << "\n# be quickly viewed in a graphical environment using \'mglview\'. "
2776  << "\n#"
2777  << "\n";
2778  }
2779 
2780  // define a helper to keep loops approximately dim-independent
2781  // since MathGL labels axes as x, y, z
2782  const std::string axes = "xyz";
2783 
2784  // (ii) write preamble and graphing tweaks
2785  out << "\n#"
2786  << "\n# Preamble."
2787  << "\n#"
2788  << "\n";
2789 
2791  out << "\nbox";
2792 
2793  // deal with dimension dependent preamble; eg. default sizes and
2794  // views for MathGL (cf. gnuplot).
2795  switch (dim)
2796  {
2797  case 2:
2798  out << "\nsetsize 800 800";
2799  out << "\nrotate 0 0";
2800  break;
2801  case 3:
2802  out << "\nsetsize 800 800";
2803  out << "\nrotate 60 40";
2804  break;
2805  default:
2806  Assert(false, ExcNotImplemented());
2807  }
2808  out << "\n";
2809 
2810  // (iii) write vertex ordering
2811  out << "\n#"
2812  << "\n# Vertex ordering."
2813  << "\n# list <vertex order> <vertex indices>"
2814  << "\n#"
2815  << "\n";
2816 
2817  // todo: This denotes the natural ordering of vertices, but it needs
2818  // to check this is really always true for a given grid (it's not
2819  // true in @ref step_1 "step-1" grid-2 for instance).
2820  switch (dim)
2821  {
2822  case 2:
2823  out << "\nlist f 0 1 2 3"
2824  << "\n";
2825  break;
2826  case 3:
2827  out
2828  << "\nlist f 0 2 4 6 | 1 3 5 7 | 0 4 1 5 | 2 6 3 7 | 0 1 2 3 | 4 5 6 7"
2829  << "\n";
2830  break;
2831  default:
2832  Assert(false, ExcNotImplemented());
2833  }
2834 
2835  // (iv) write a list of vertices of cells
2836  out << "\n#"
2837  << "\n# List of vertices."
2838  << "\n# list <id> <vertices>"
2839  << "\n#"
2840  << "\n";
2841 
2842  // run over all active cells and write out a list of
2843  // xyz-coordinates that correspond to vertices
2844  typename ::Triangulation<dim, spacedim>::active_cell_iterator
2845  cell = tria.begin_active(),
2846  endc = tria.end();
2847 
2848  // No global indices in deal.II, so we make one up here.
2849  for (; cell != endc; ++cell)
2850  {
2851  for (unsigned int i = 0; i < dim; ++i)
2852  {
2853  // if (cell->direction_flag ()==true)
2854  // out << "\ntrue";
2855  // else
2856  // out << "\nfalse";
2857 
2858  out << "\nlist " << axes[i] << cell->active_cell_index() << " ";
2859  for (unsigned int j = 0; j < GeometryInfo<dim>::vertices_per_cell;
2860  ++j)
2861  out << cell->vertex(j)[i] << " ";
2862  }
2863  out << '\n';
2864  }
2865 
2866  // (v) write out cells to plot as quadplot objects
2867  out << "\n#"
2868  << "\n# List of cells to quadplot."
2869  << "\n# quadplot <vertex order> <id> <style>"
2870  << "\n#"
2871  << "\n";
2872  for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
2873  {
2874  out << "\nquadplot f ";
2875  for (unsigned int j = 0; j < dim; ++j)
2876  out << axes[j] << i << " ";
2877  out << "\'k#\'";
2878  }
2879  out << "\n";
2880 
2881  // (vi) write footer
2882  out << "\n#"
2883  << "\n#"
2884  << "\n#"
2885  << "\n";
2886 
2887  // make sure everything now gets to the output stream
2888  out.flush();
2889  AssertThrow(out, ExcIO());
2890 }
2891 
2892 
2893 
2894 namespace
2895 {
2902  template <int dim, int spacedim, typename ITERATOR, typename END>
2903  void
2904  generate_triangulation_patches(
2905  std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
2906  ITERATOR cell,
2907  END end)
2908  {
2909  // convert each of the active cells into a patch
2910  for (; cell != end; ++cell)
2911  {
2913  patch.n_subdivisions = 1;
2915 
2916  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2917  {
2918  patch.vertices[v] = cell->vertex(v);
2919  patch.data(0, v) = cell->level();
2920  patch.data(1, v) =
2921  static_cast<std::make_signed<types::manifold_id>::type>(
2922  cell->manifold_id());
2923  patch.data(2, v) =
2924  static_cast<std::make_signed<types::material_id>::type>(
2925  cell->material_id());
2926  if (!cell->has_children())
2927  patch.data(3, v) =
2928  static_cast<std::make_signed<types::subdomain_id>::type>(
2929  cell->subdomain_id());
2930  else
2931  patch.data(3, v) = -1;
2932  patch.data(4, v) =
2933  static_cast<std::make_signed<types::subdomain_id>::type>(
2934  cell->level_subdomain_id());
2935  }
2936  patches.push_back(patch);
2937  }
2938  }
2939 
2940 
2941 
2942  std::vector<std::string>
2943  triangulation_patch_data_names()
2944  {
2945  std::vector<std::string> v(5);
2946  v[0] = "level";
2947  v[1] = "manifold";
2948  v[2] = "material";
2949  v[3] = "subdomain";
2950  v[4] = "level_subdomain";
2951  return v;
2952  }
2953 
2957  std::vector<typename Triangulation<3, 3>::active_line_iterator>
2958  get_boundary_edge_iterators(const Triangulation<3, 3> &tria)
2959  {
2960  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
2961 
2962  std::vector<bool> flags;
2963  tria.save_user_flags_line(flags);
2964  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
2965 
2966  for (auto face : tria.active_face_iterators())
2967  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_face; ++l)
2968  {
2969  const auto line = face->line(l);
2970  if (line->user_flag_set() || line->has_children())
2971  continue;
2972  else
2973  line->set_user_flag();
2974  if (line->at_boundary())
2975  res.emplace_back(line);
2976  }
2977  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
2978  return res;
2979  }
2980 
2981 
2982 
2986  template <int dim, int spacedim>
2987  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
2988  get_boundary_edge_iterators(const Triangulation<dim, spacedim> &)
2989  {
2990  return {};
2991  }
2992 
2993 
2994 
2999  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3000  get_relevant_edge_iterators(const Triangulation<3, 3> &tria)
3001  {
3002  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3003 
3004  std::vector<bool> flags;
3005  tria.save_user_flags_line(flags);
3006  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3007 
3008  for (auto face : tria.active_face_iterators())
3009  for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_face; ++l)
3010  {
3011  const auto line = face->line(l);
3012  if (line->user_flag_set() || line->has_children())
3013  continue;
3014  else
3015  line->set_user_flag();
3016  if (line->manifold_id() != numbers::flat_manifold_id ||
3017  (line->boundary_id() != 0 &&
3018  line->boundary_id() != numbers::invalid_boundary_id))
3019  res.emplace_back(line);
3020  }
3021  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3022  return res;
3023  }
3024 
3025 
3029  template <int dim, int spacedim>
3030  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3031  get_relevant_edge_iterators(const Triangulation<dim, spacedim> &)
3032  {
3033  return {};
3034  }
3035 
3036 
3037 
3041  template <int dim, int spacedim>
3042  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3043  get_boundary_face_iterators(const Triangulation<dim, spacedim> &tria)
3044  {
3045  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3046  res;
3047  if (dim == 1)
3048  return res;
3049  for (auto face : tria.active_face_iterators())
3050  {
3051  if (face->boundary_id() != numbers::invalid_boundary_id)
3052  res.push_back(face);
3053  }
3054  return res;
3055  }
3056 
3057 
3058 
3063  template <int dim, int spacedim>
3064  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3065  get_relevant_face_iterators(const Triangulation<dim, spacedim> &tria)
3066  {
3067  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3068  res;
3069  if (dim == 1)
3070  return res;
3071  for (auto face : tria.active_face_iterators())
3072  {
3073  if (face->manifold_id() != numbers::flat_manifold_id ||
3074  (face->boundary_id() != 0 &&
3075  face->boundary_id() != numbers::invalid_boundary_id))
3076  res.push_back(face);
3077  }
3078  return res;
3079  }
3080 } // namespace
3081 
3082 
3083 
3084 template <int dim, int spacedim>
3085 void
3087  std::ostream & out) const
3088 {
3089  AssertThrow(out, ExcIO());
3090 
3091  // get the positions of the vertices
3092  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
3093 
3094  const auto n_vertices = vertices.size();
3095 
3096  out << "# vtk DataFile Version 3.0\n"
3097  << "Triangulation generated with deal.II\n"
3098  << "ASCII\n"
3099  << "DATASET UNSTRUCTURED_GRID\n"
3100  << "POINTS " << n_vertices << " double\n";
3101 
3102  // actually write the vertices.
3103  for (const auto &v : vertices)
3104  {
3105  out << v;
3106  for (unsigned int d = spacedim + 1; d <= 3; ++d)
3107  out << " 0"; // fill with zeroes
3108  out << '\n';
3109  }
3110 
3111  const auto faces = vtk_flags.output_only_relevant ?
3112  get_relevant_face_iterators(tria) :
3113  get_boundary_face_iterators(tria);
3114  const auto edges = vtk_flags.output_only_relevant ?
3115  get_relevant_edge_iterators(tria) :
3116  get_boundary_edge_iterators(tria);
3117 
3118  AssertThrow(
3119  vtk_flags.output_cells || (dim >= 2 && vtk_flags.output_faces) ||
3120  (dim >= 3 && vtk_flags.output_edges),
3121  ExcMessage(
3122  "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3123 
3124  // Write cells preamble
3125  const int n_cells = (vtk_flags.output_cells ? tria.n_active_cells() : 0) +
3126  (vtk_flags.output_faces ? faces.size() : 0) +
3127  (vtk_flags.output_edges ? edges.size() : 0);
3128 
3129  // VTK now expects a number telling the total storage requirement to read all
3130  // cell connectivity information. The connectivity information is read cell by
3131  // cell, first specifying how many vertices are required to describe the cell,
3132  // and then specifying the index of every vertex. This means that for every
3133  // deal.II object type, we always need n_vertices + 1 integer per cell.
3134  // Compute the total number here.
3135  const int cells_size =
3138  0) +
3140  faces.size() * (GeometryInfo<dim>::vertices_per_face + 1) :
3141  0) +
3142  (vtk_flags.output_edges ? edges.size() * (3) :
3143  0); // only in 3d, otherwise it is always zero.
3144 
3145  AssertThrow(cells_size > 0, ExcMessage("No cells given to be output!"));
3146 
3147  out << "\nCELLS " << n_cells << ' ' << cells_size << '\n';
3148  /*
3149  * VTK cells:
3150  *
3151  * 1 VTK_VERTEX
3152  * 3 VTK_LINE
3153  * 9 VTK_QUAD
3154  * 12 VTK_HEXAHEDRON
3155  * ...
3156  */
3157  const int cell_type = (dim == 1 ? 3 : dim == 2 ? 9 : 12);
3158  const int face_type = (dim == 1 ? 1 : dim == 2 ? 3 : 9);
3159  const int co_face_type = (dim == 1 ? -1 : dim == 2 ? -1 : 3);
3160 
3161  // write cells.
3162  if (vtk_flags.output_cells)
3163  for (const auto &cell : tria.active_cell_iterators())
3164  {
3165  out << GeometryInfo<dim>::vertices_per_cell;
3166  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell; ++i)
3167  {
3168  out << ' ' << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[i]);
3169  }
3170  out << '\n';
3171  }
3172  if (vtk_flags.output_faces)
3173  for (const auto &face : faces)
3174  {
3175  out << GeometryInfo<dim>::vertices_per_face;
3176  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
3177  {
3178  out << ' '
3179  << face->vertex_index(GeometryInfo < (dim > 1) ?
3180  dim - 1 :
3181  dim > ::ucd_to_deal[i]);
3182  }
3183  out << '\n';
3184  }
3185  if (vtk_flags.output_edges)
3186  for (const auto &edge : edges)
3187  {
3188  out << 2;
3189  for (unsigned int i = 0; i < 2; ++i)
3190  out << ' ' << edge->vertex_index(i);
3191  out << '\n';
3192  }
3193 
3194  // write cell types
3195  out << "\nCELL_TYPES " << n_cells << '\n';
3196  if (vtk_flags.output_cells)
3197  {
3198  for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
3199  {
3200  out << cell_type << ' ';
3201  }
3202  out << '\n';
3203  }
3204  if (vtk_flags.output_faces)
3205  {
3206  for (unsigned int i = 0; i < faces.size(); ++i)
3207  {
3208  out << face_type << ' ';
3209  }
3210  out << '\n';
3211  }
3212  if (vtk_flags.output_edges)
3213  {
3214  for (unsigned int i = 0; i < edges.size(); ++i)
3215  {
3216  out << co_face_type << ' ';
3217  }
3218  }
3219  out << "\n\nCELL_DATA " << n_cells << '\n'
3220  << "SCALARS MaterialID int 1\n"
3221  << "LOOKUP_TABLE default\n";
3222 
3223  // Now material id and boundary id
3224  if (vtk_flags.output_cells)
3225  {
3226  for (const auto &cell : tria.active_cell_iterators())
3227  {
3228  out << static_cast<std::make_signed<types::material_id>::type>(
3229  cell->material_id())
3230  << ' ';
3231  }
3232  out << '\n';
3233  }
3234  if (vtk_flags.output_faces)
3235  {
3236  for (const auto &face : faces)
3237  {
3238  out << static_cast<std::make_signed<types::boundary_id>::type>(
3239  face->boundary_id())
3240  << ' ';
3241  }
3242  out << '\n';
3243  }
3244  if (vtk_flags.output_edges)
3245  {
3246  for (const auto &edge : edges)
3247  {
3248  out << static_cast<std::make_signed<types::boundary_id>::type>(
3249  edge->boundary_id())
3250  << ' ';
3251  }
3252  }
3253 
3254  out << "\n\nSCALARS ManifoldID int 1\n"
3255  << "LOOKUP_TABLE default\n";
3256 
3257  // Now material id and boundary id
3258  if (vtk_flags.output_cells)
3259  {
3260  for (const auto &cell : tria.active_cell_iterators())
3261  {
3262  out << static_cast<std::make_signed<types::boundary_id>::type>(
3263  cell->manifold_id())
3264  << ' ';
3265  }
3266  out << '\n';
3267  }
3268  if (vtk_flags.output_faces)
3269  {
3270  for (const auto &face : faces)
3271  {
3272  out << static_cast<std::make_signed<types::boundary_id>::type>(
3273  face->manifold_id())
3274  << ' ';
3275  }
3276  out << '\n';
3277  }
3278  if (vtk_flags.output_edges)
3279  {
3280  for (const auto &edge : edges)
3281  {
3282  out << static_cast<std::make_signed<types::boundary_id>::type>(
3283  edge->manifold_id())
3284  << ' ';
3285  }
3286  out << '\n';
3287  }
3288 
3289  out.flush();
3290 
3291  AssertThrow(out, ExcIO());
3292 }
3293 
3294 
3295 
3296 template <int dim, int spacedim>
3297 void
3299  std::ostream & out) const
3300 {
3301  AssertThrow(out, ExcIO());
3302 
3303  // convert the cells of the triangulation into a set of patches
3304  // and then have them output. since there is no data attached to
3305  // the geometry, we also do not have to provide any names, identifying
3306  // information, etc.
3307  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3308  patches.reserve(tria.n_active_cells());
3309  generate_triangulation_patches(patches, tria.begin_active(), tria.end());
3311  patches,
3312  triangulation_patch_data_names(),
3313  std::vector<
3314  std::tuple<unsigned int,
3315  unsigned int,
3316  std::string,
3318  vtu_flags,
3319  out);
3320 
3321  AssertThrow(out, ExcIO());
3322 }
3323 
3324 
3325 
3326 template <int dim, int spacedim>
3327 void
3329  const Triangulation<dim, spacedim> &tria,
3330  const std::string & filename_without_extension,
3331  const bool view_levels,
3332  const bool include_artificial) const
3333 {
3334  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3335  const unsigned int n_datasets = 4;
3336  std::vector<std::string> data_names;
3337  data_names.emplace_back("level");
3338  data_names.emplace_back("subdomain");
3339  data_names.emplace_back("level_subdomain");
3340  data_names.emplace_back("proc_writing");
3341 
3342  const unsigned int n_q_points = GeometryInfo<dim>::vertices_per_cell;
3343 
3344  typename Triangulation<dim, spacedim>::cell_iterator cell, endc;
3345  for (cell = tria.begin(), endc = tria.end(); cell != endc; ++cell)
3346  {
3347  if (!view_levels)
3348  {
3349  if (cell->has_children())
3350  continue;
3351  if (!include_artificial &&
3352  cell->subdomain_id() == numbers::artificial_subdomain_id)
3353  continue;
3354  }
3355  else if (!include_artificial)
3356  {
3357  if (cell->has_children() &&
3358  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
3359  continue;
3360  else if (!cell->has_children() &&
3361  cell->level_subdomain_id() ==
3363  cell->subdomain_id() == numbers::artificial_subdomain_id)
3364  continue;
3365  }
3366 
3368  patch.data.reinit(n_datasets, n_q_points);
3369  patch.points_are_available = false;
3370 
3371  for (unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3372  {
3373  patch.vertices[vertex] = cell->vertex(vertex);
3374  patch.data(0, vertex) = cell->level();
3375  if (!cell->has_children())
3376  patch.data(1, vertex) = static_cast<double>(
3377  static_cast<std::make_signed<types::subdomain_id>::type>(
3378  cell->subdomain_id()));
3379  else
3380  patch.data(1, vertex) = -1.0;
3381  patch.data(2, vertex) = static_cast<double>(
3382  static_cast<std::make_signed<types::subdomain_id>::type>(
3383  cell->level_subdomain_id()));
3384  patch.data(3, vertex) = tria.locally_owned_subdomain();
3385  }
3386 
3387  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
3389  patches.push_back(patch);
3390  }
3391 
3392  // only create .pvtu file if running in parallel
3393  // if not, just create a .vtu file with no reference
3394  // to the processor number
3395  std::string new_file = filename_without_extension + ".vtu";
3397  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
3398  {
3399  new_file = filename_without_extension + ".proc" +
3400  Utilities::int_to_string(tr->locally_owned_subdomain(), 4) +
3401  ".vtu";
3402 
3403  // create .pvtu record
3404  if (tr->locally_owned_subdomain() == 0)
3405  {
3406  std::vector<std::string> filenames;
3407 
3408  // .pvtu needs to reference the files without a relative path because
3409  // it will be written in the same directory. For this, remove any
3410  // paths from filename.
3411  std::size_t pos = filename_without_extension.find_last_of('/');
3412  if (pos == std::string::npos)
3413  pos = 0;
3414  else
3415  pos += 1;
3416  const unsigned int n_procs =
3417  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3418  for (unsigned int i = 0; i < n_procs; ++i)
3419  filenames.push_back(filename_without_extension.substr(pos) +
3420  ".proc" + Utilities::int_to_string(i, 4) +
3421  ".vtu");
3422 
3423  const std::string pvtu_master_filename =
3424  (filename_without_extension + ".pvtu");
3425  std::ofstream pvtu_master(pvtu_master_filename.c_str());
3426 
3428  data_out.attach_triangulation(*tr);
3429 
3430  // We need a dummy vector with the names of the data values in the
3431  // .vtu files in order that the .pvtu contains reference these values
3432  Vector<float> dummy_vector(tr->n_active_cells());
3433  data_out.add_data_vector(dummy_vector, "level");
3434  data_out.add_data_vector(dummy_vector, "subdomain");
3435  data_out.add_data_vector(dummy_vector, "level_subdomain");
3436  data_out.add_data_vector(dummy_vector, "proc_writing");
3437 
3438  data_out.build_patches();
3439 
3440  data_out.write_pvtu_record(pvtu_master, filenames);
3441  }
3442  }
3443 
3444  std::ofstream out(new_file.c_str());
3445  std::vector<
3446  std::tuple<unsigned int,
3447  unsigned int,
3448  std::string,
3450  vector_data_ranges;
3451  DataOutBase::VtkFlags flags;
3452  DataOutBase::write_vtu(patches, data_names, vector_data_ranges, flags, out);
3453 }
3454 
3455 
3456 
3457 unsigned int
3459 {
3460  return 0;
3461 }
3462 
3463 unsigned int
3465 {
3466  return 0;
3467 }
3468 
3469 
3470 unsigned int
3472 {
3473  return 0;
3474 }
3475 
3476 unsigned int
3478 {
3479  return 0;
3480 }
3481 
3482 unsigned int
3484 {
3485  return 0;
3486 }
3487 
3488 unsigned int
3490 {
3491  return 0;
3492 }
3493 
3494 unsigned int
3496 {
3497  return 0;
3498 }
3499 
3500 unsigned int
3502 {
3503  return 0;
3504 }
3505 
3506 
3507 
3508 template <int dim, int spacedim>
3509 unsigned int
3511 {
3513  unsigned int n_faces = 0;
3514 
3515  for (face = tria.begin_active_face(), endf = tria.end_face(); face != endf;
3516  ++face)
3517  if ((face->at_boundary()) && (face->boundary_id() != 0))
3518  n_faces++;
3519 
3520  return n_faces;
3521 }
3522 
3523 
3524 
3525 template <int dim, int spacedim>
3526 unsigned int
3528 {
3529  // save the user flags for lines so
3530  // we can use these flags to track
3531  // which ones we've already counted
3532  std::vector<bool> line_flags;
3533  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3534  line_flags);
3535  const_cast<::Triangulation<dim, spacedim> &>(tria)
3536  .clear_user_flags_line();
3537 
3538  unsigned int n_lines = 0;
3539 
3541 
3542  for (cell = tria.begin_active(), endc = tria.end(); cell != endc; ++cell)
3543  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3544  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3545  (cell->line(l)->user_flag_set() == false))
3546  {
3547  ++n_lines;
3548  cell->line(l)->set_user_flag();
3549  }
3550 
3551  // at the end, restore the user
3552  // flags for the lines
3553  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3554  line_flags);
3555 
3556  return n_lines;
3557 }
3558 
3559 
3560 
3561 unsigned int
3563  const unsigned int next_element_index,
3564  std::ostream &) const
3565 {
3566  return next_element_index;
3567 }
3568 
3569 
3570 unsigned int
3572  const unsigned int next_element_index,
3573  std::ostream &) const
3574 {
3575  return next_element_index;
3576 }
3577 
3578 unsigned int
3580  const unsigned int next_element_index,
3581  std::ostream &) const
3582 {
3583  return next_element_index;
3584 }
3585 
3586 
3587 unsigned int
3589  const unsigned int next_element_index,
3590  std::ostream &) const
3591 {
3592  return next_element_index;
3593 }
3594 
3595 unsigned int
3597  const unsigned int next_element_index,
3598  std::ostream &) const
3599 {
3600  return next_element_index;
3601 }
3602 
3603 
3604 unsigned int
3606  const unsigned int next_element_index,
3607  std::ostream &) const
3608 {
3609  return next_element_index;
3610 }
3611 
3612 
3613 unsigned int
3615  const unsigned int next_element_index,
3616  std::ostream &) const
3617 {
3618  return next_element_index;
3619 }
3620 
3621 unsigned int
3623  const unsigned int next_element_index,
3624  std::ostream &) const
3625 {
3626  return next_element_index;
3627 }
3628 
3629 
3630 
3631 template <int dim, int spacedim>
3632 unsigned int
3634  const unsigned int next_element_index,
3635  std::ostream & out) const
3636 {
3637  unsigned int current_element_index = next_element_index;
3639 
3640  for (face = tria.begin_active_face(), endf = tria.end_face(); face != endf;
3641  ++face)
3642  if (face->at_boundary() && (face->boundary_id() != 0))
3643  {
3644  out << current_element_index << ' ';
3645  switch (dim)
3646  {
3647  case 2:
3648  out << 1 << ' ';
3649  break;
3650  case 3:
3651  out << 3 << ' ';
3652  break;
3653  default:
3654  Assert(false, ExcNotImplemented());
3655  }
3656  out << static_cast<unsigned int>(face->boundary_id()) << ' '
3657  << static_cast<unsigned int>(face->boundary_id()) << ' '
3659  // note: vertex numbers are 1-base
3660  for (unsigned int vertex = 0;
3661  vertex < GeometryInfo<dim>::vertices_per_face;
3662  ++vertex)
3663  out << ' '
3664  << face->vertex_index(
3666  1;
3667  out << '\n';
3668 
3669  ++current_element_index;
3670  }
3671  return current_element_index;
3672 }
3673 
3674 
3675 template <int dim, int spacedim>
3676 unsigned int
3678  const unsigned int next_element_index,
3679  std::ostream & out) const
3680 {
3681  unsigned int current_element_index = next_element_index;
3682  // save the user flags for lines so
3683  // we can use these flags to track
3684  // which ones we've already taken
3685  // care of
3686  std::vector<bool> line_flags;
3687  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3688  line_flags);
3689  const_cast<::Triangulation<dim, spacedim> &>(tria)
3690  .clear_user_flags_line();
3691 
3693 
3694  for (cell = tria.begin_active(), endc = tria.end(); cell != endc; ++cell)
3695  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3696  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3697  (cell->line(l)->user_flag_set() == false))
3698  {
3699  out << next_element_index << " 1 ";
3700  out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
3701  << static_cast<unsigned int>(cell->line(l)->boundary_id())
3702  << " 2 ";
3703  // note: vertex numbers are 1-base
3704  for (unsigned int vertex = 0; vertex < 2; ++vertex)
3705  out << ' '
3706  << cell->line(l)->vertex_index(
3708  1;
3709  out << '\n';
3710 
3711  // move on to the next line
3712  // but mark the current one
3713  // as taken care of
3714  ++current_element_index;
3715  cell->line(l)->set_user_flag();
3716  }
3717 
3718  // at the end, restore the user
3719  // flags for the lines
3720  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3721  line_flags);
3722 
3723  return current_element_index;
3724 }
3725 
3726 
3727 
3728 unsigned int
3730  const unsigned int next_element_index,
3731  std::ostream &) const
3732 {
3733  return next_element_index;
3734 }
3735 
3736 unsigned int
3738  const unsigned int next_element_index,
3739  std::ostream &) const
3740 {
3741  return next_element_index;
3742 }
3743 
3744 unsigned int
3746  const unsigned int next_element_index,
3747  std::ostream &) const
3748 {
3749  return next_element_index;
3750 }
3751 
3752 unsigned int
3754  const unsigned int next_element_index,
3755  std::ostream &) const
3756 {
3757  return next_element_index;
3758 }
3759 
3760 unsigned int
3762  const unsigned int next_element_index,
3763  std::ostream &) const
3764 {
3765  return next_element_index;
3766 }
3767 
3768 
3769 unsigned int
3771  const unsigned int next_element_index,
3772  std::ostream &) const
3773 {
3774  return next_element_index;
3775 }
3776 
3777 
3778 unsigned int
3780  const unsigned int next_element_index,
3781  std::ostream &) const
3782 {
3783  return next_element_index;
3784 }
3785 
3786 unsigned int
3788  const unsigned int next_element_index,
3789  std::ostream &) const
3790 {
3791  return next_element_index;
3792 }
3793 
3794 
3795 
3796 template <int dim, int spacedim>
3797 unsigned int
3799  const unsigned int next_element_index,
3800  std::ostream & out) const
3801 {
3802  unsigned int current_element_index = next_element_index;
3804 
3805  for (face = tria.begin_active_face(), endf = tria.end_face(); face != endf;
3806  ++face)
3807  if (face->at_boundary() && (face->boundary_id() != 0))
3808  {
3809  out << current_element_index << " "
3810  << static_cast<unsigned int>(face->boundary_id()) << " ";
3811  switch (dim)
3812  {
3813  case 2:
3814  out << "line ";
3815  break;
3816  case 3:
3817  out << "quad ";
3818  break;
3819  default:
3820  Assert(false, ExcNotImplemented());
3821  }
3822  // note: vertex numbers are 1-base
3823  for (unsigned int vertex = 0;
3824  vertex < GeometryInfo<dim>::vertices_per_face;
3825  ++vertex)
3826  out << face->vertex_index(
3828  1
3829  << ' ';
3830  out << '\n';
3831 
3832  ++current_element_index;
3833  }
3834  return current_element_index;
3835 }
3836 
3837 
3838 
3839 template <int dim, int spacedim>
3840 unsigned int
3842  const unsigned int next_element_index,
3843  std::ostream & out) const
3844 {
3845  unsigned int current_element_index = next_element_index;
3846  // save the user flags for lines so
3847  // we can use these flags to track
3848  // which ones we've already taken
3849  // care of
3850  std::vector<bool> line_flags;
3851  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3852  line_flags);
3853  const_cast<::Triangulation<dim, spacedim> &>(tria)
3854  .clear_user_flags_line();
3855 
3857 
3858  for (cell = tria.begin_active(), endc = tria.end(); cell != endc; ++cell)
3859  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3860  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3861  (cell->line(l)->user_flag_set() == false))
3862  {
3863  out << current_element_index << " "
3864  << static_cast<unsigned int>(cell->line(l)->boundary_id())
3865  << " line ";
3866  // note: vertex numbers in ucd format are 1-base
3867  for (unsigned int vertex = 0; vertex < 2; ++vertex)
3868  out << cell->line(l)->vertex_index(
3870  1
3871  << ' ';
3872  out << '\n';
3873 
3874  // move on to the next line
3875  // but mark the current one
3876  // as taken care of
3877  ++current_element_index;
3878  cell->line(l)->set_user_flag();
3879  }
3880 
3881  // at the end, restore the user
3882  // flags for the lines
3883  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3884  line_flags);
3885  return current_element_index;
3886 }
3887 
3888 
3889 
3890 namespace internal
3891 {
3892  namespace
3893  {
3902  template <int spacedim>
3903  void
3904  remove_colinear_points(std::vector<Point<spacedim>> &points)
3905  {
3906  while (points.size() > 2)
3907  {
3908  Tensor<1, spacedim> first_difference = points[1] - points[0];
3909  first_difference /= first_difference.norm();
3910  Tensor<1, spacedim> second_difference = points[2] - points[1];
3911  second_difference /= second_difference.norm();
3912  // If the three points are colinear then remove the middle one.
3913  if ((first_difference - second_difference).norm() < 1e-10)
3914  points.erase(points.begin() + 1);
3915  else
3916  break;
3917  }
3918  }
3919 
3920 
3921 
3922  template <int spacedim>
3923  void
3924  write_gnuplot(const ::Triangulation<1, spacedim> &tria,
3925  std::ostream & out,
3926  const Mapping<1, spacedim> *,
3928  {
3929  AssertThrow(out, ExcIO());
3930 
3931  const int dim = 1;
3932 
3933  typename ::Triangulation<dim, spacedim>::active_cell_iterator cell =
3934  tria.begin_active();
3935  const typename ::Triangulation<dim, spacedim>::active_cell_iterator
3936  endc = tria.end();
3937  for (; cell != endc; ++cell)
3938  {
3939  if (gnuplot_flags.write_cell_numbers)
3940  out << "# cell " << cell << '\n';
3941 
3942  out << cell->vertex(0) << ' ' << cell->level() << ' '
3943  << static_cast<unsigned int>(cell->material_id()) << '\n'
3944  << cell->vertex(1) << ' ' << cell->level() << ' '
3945  << static_cast<unsigned int>(cell->material_id()) << '\n'
3946  << "\n\n";
3947  }
3948 
3949  // make sure everything now gets to
3950  // disk
3951  out.flush();
3952 
3953  AssertThrow(out, ExcIO());
3954  }
3955 
3956 
3957 
3958  template <int spacedim>
3959  void
3960  write_gnuplot(const ::Triangulation<2, spacedim> &tria,
3961  std::ostream & out,
3962  const Mapping<2, spacedim> * mapping,
3963  const GridOutFlags::Gnuplot & gnuplot_flags)
3964  {
3965  AssertThrow(out, ExcIO());
3966 
3967  const int dim = 2;
3968 
3969  const unsigned int n_additional_points =
3970  gnuplot_flags.n_boundary_face_points;
3971  const unsigned int n_points = 2 + n_additional_points;
3972 
3973  typename ::Triangulation<dim, spacedim>::active_cell_iterator cell =
3974  tria.begin_active();
3975  const typename ::Triangulation<dim, spacedim>::active_cell_iterator
3976  endc = tria.end();
3977 
3978  // If we need to plot curved lines then generate a quadrature formula to
3979  // place points via the mapping
3980  Quadrature<dim> * q_projector = nullptr;
3981  std::vector<Point<dim - 1>> boundary_points;
3982  if (mapping != nullptr)
3983  {
3984  boundary_points.resize(n_points);
3985  boundary_points[0][0] = 0;
3986  boundary_points[n_points - 1][0] = 1;
3987  for (unsigned int i = 1; i < n_points - 1; ++i)
3988  boundary_points[i](0) = 1. * i / (n_points - 1);
3989 
3990  std::vector<double> dummy_weights(n_points, 1. / n_points);
3991  Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
3992 
3993  q_projector = new Quadrature<dim>(
3995  }
3996 
3997  for (; cell != endc; ++cell)
3998  {
3999  if (gnuplot_flags.write_cell_numbers)
4000  out << "# cell " << cell << '\n';
4001 
4002  if (mapping == nullptr ||
4003  (dim == spacedim ?
4004  (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4005  // ignore checking for boundary or interior cells in the codim
4006  // 1 case: 'or false' is a no-op
4007  false))
4008  {
4009  // write out the four sides of this cell by putting the four
4010  // points (+ the initial point again) in a row and lifting the
4011  // drawing pencil at the end
4012  for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_cell;
4013  ++i)
4014  out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i]) << ' '
4015  << cell->level() << ' '
4016  << static_cast<unsigned int>(cell->material_id()) << '\n';
4017  out << cell->vertex(0) << ' ' << cell->level() << ' '
4018  << static_cast<unsigned int>(cell->material_id()) << '\n'
4019  << '\n' // double new line for gnuplot 3d plots
4020  << '\n';
4021  }
4022  else
4023  // cell is at boundary and we are to treat curved boundaries. so
4024  // loop over all faces and draw them as small pieces of lines
4025  {
4026  for (unsigned int face_no = 0;
4027  face_no < GeometryInfo<dim>::faces_per_cell;
4028  ++face_no)
4029  {
4030  const typename ::Triangulation<dim,
4031  spacedim>::face_iterator
4032  face = cell->face(face_no);
4033  if (dim != spacedim || face->at_boundary() ||
4034  gnuplot_flags.curved_inner_cells)
4035  {
4036  // Save the points on each face to a vector and then try
4037  // to remove colinear points that won't show up in the
4038  // generated plot.
4039  std::vector<Point<spacedim>> line_points;
4040  // compute offset of quadrature points within set of
4041  // projected points
4042  const unsigned int offset = face_no * n_points;
4043  for (unsigned int i = 0; i < n_points; ++i)
4044  line_points.push_back(
4045  mapping->transform_unit_to_real_cell(
4046  cell, q_projector->point(offset + i)));
4047  internal::remove_colinear_points(line_points);
4048 
4049  for (const Point<spacedim> &point : line_points)
4050  out << point << ' ' << cell->level() << ' '
4051  << static_cast<unsigned int>(cell->material_id())
4052  << '\n';
4053 
4054  out << '\n' << '\n';
4055  }
4056  else
4057  {
4058  // if, however, the face is not at the boundary and we
4059  // don't want to curve anything, then draw it as usual
4060  out << face->vertex(0) << ' ' << cell->level() << ' '
4061  << static_cast<unsigned int>(cell->material_id())
4062  << '\n'
4063  << face->vertex(1) << ' ' << cell->level() << ' '
4064  << static_cast<unsigned int>(cell->material_id())
4065  << '\n'
4066  << '\n'
4067  << '\n';
4068  }
4069  }
4070  }
4071  }
4072 
4073  if (q_projector != nullptr)
4074  delete q_projector;
4075 
4076  // make sure everything now gets to disk
4077  out.flush();
4078 
4079  AssertThrow(out, ExcIO());
4080  }
4081 
4082 
4083 
4084  template <int spacedim>
4085  void
4086  write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4087  std::ostream & out,
4088  const Mapping<3, spacedim> * mapping,
4089  const GridOutFlags::Gnuplot & gnuplot_flags)
4090  {
4091  AssertThrow(out, ExcIO());
4092 
4093  const int dim = 3;
4094 
4095  const unsigned int n_additional_points =
4096  gnuplot_flags.n_boundary_face_points;
4097  const unsigned int n_points = 2 + n_additional_points;
4098 
4099  typename ::Triangulation<dim, spacedim>::active_cell_iterator cell =
4100  tria.begin_active();
4101  const typename ::Triangulation<dim, spacedim>::active_cell_iterator
4102  endc = tria.end();
4103 
4104  // If we need to plot curved lines then generate a quadrature formula to
4105  // place points via the mapping
4106  Quadrature<dim> * q_projector = nullptr;
4107  std::vector<Point<1>> boundary_points;
4108  if (mapping != nullptr)
4109  {
4110  boundary_points.resize(n_points);
4111  boundary_points[0][0] = 0;
4112  boundary_points[n_points - 1][0] = 1;
4113  for (unsigned int i = 1; i < n_points - 1; ++i)
4114  boundary_points[i](0) = 1. * i / (n_points - 1);
4115 
4116  std::vector<double> dummy_weights(n_points, 1. / n_points);
4117  Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4118 
4119  // tensor product of points, only one copy
4120  QIterated<dim - 1> quadrature(quadrature1d, 1);
4121  q_projector = new Quadrature<dim>(
4123  }
4124 
4125  for (; cell != endc; ++cell)
4126  {
4127  if (gnuplot_flags.write_cell_numbers)
4128  out << "# cell " << cell << '\n';
4129 
4130  if (mapping == nullptr || n_points == 2 ||
4131  (!cell->has_boundary_lines() &&
4132  !gnuplot_flags.curved_inner_cells))
4133  {
4134  // front face
4135  out << cell->vertex(0) << ' ' << cell->level() << ' '
4136  << static_cast<unsigned int>(cell->material_id()) << '\n'
4137  << cell->vertex(1) << ' ' << cell->level() << ' '
4138  << static_cast<unsigned int>(cell->material_id()) << '\n'
4139  << cell->vertex(5) << ' ' << cell->level() << ' '
4140  << static_cast<unsigned int>(cell->material_id()) << '\n'
4141  << cell->vertex(4) << ' ' << cell->level() << ' '
4142  << static_cast<unsigned int>(cell->material_id()) << '\n'
4143  << cell->vertex(0) << ' ' << cell->level() << ' '
4144  << static_cast<unsigned int>(cell->material_id()) << '\n'
4145  << '\n';
4146  // back face
4147  out << cell->vertex(2) << ' ' << cell->level() << ' '
4148  << static_cast<unsigned int>(cell->material_id()) << '\n'
4149  << cell->vertex(3) << ' ' << cell->level() << ' '
4150  << static_cast<unsigned int>(cell->material_id()) << '\n'
4151  << cell->vertex(7) << ' ' << cell->level() << ' '
4152  << static_cast<unsigned int>(cell->material_id()) << '\n'
4153  << cell->vertex(6) << ' ' << cell->level() << ' '
4154  << static_cast<unsigned int>(cell->material_id()) << '\n'
4155  << cell->vertex(2) << ' ' << cell->level() << ' '
4156  << static_cast<unsigned int>(cell->material_id()) << '\n'
4157  << '\n';
4158 
4159  // now for the four connecting lines
4160  out << cell->vertex(0) << ' ' << cell->level() << ' '
4161  << static_cast<unsigned int>(cell->material_id()) << '\n'
4162  << cell->vertex(2) << ' ' << cell->level() << ' '
4163  << static_cast<unsigned int>(cell->material_id()) << '\n'
4164  << '\n';
4165  out << cell->vertex(1) << ' ' << cell->level() << ' '
4166  << static_cast<unsigned int>(cell->material_id()) << '\n'
4167  << cell->vertex(3) << ' ' << cell->level() << ' '
4168  << static_cast<unsigned int>(cell->material_id()) << '\n'
4169  << '\n';
4170  out << cell->vertex(5) << ' ' << cell->level() << ' '
4171  << static_cast<unsigned int>(cell->material_id()) << '\n'
4172  << cell->vertex(7) << ' ' << cell->level() << ' '
4173  << static_cast<unsigned int>(cell->material_id()) << '\n'
4174  << '\n';
4175  out << cell->vertex(4) << ' ' << cell->level() << ' '
4176  << static_cast<unsigned int>(cell->material_id()) << '\n'
4177  << cell->vertex(6) << ' ' << cell->level() << ' '
4178  << static_cast<unsigned int>(cell->material_id()) << '\n'
4179  << '\n';
4180  }
4181  else
4182  {
4183  for (unsigned int face_no = 0;
4184  face_no < GeometryInfo<dim>::faces_per_cell;
4185  ++face_no)
4186  {
4187  const typename ::Triangulation<dim,
4188  spacedim>::face_iterator
4189  face = cell->face(face_no);
4190 
4191  if (face->at_boundary() &&
4192  gnuplot_flags.write_additional_boundary_lines)
4193  {
4194  const unsigned int offset = face_no * n_points * n_points;
4195  for (unsigned int i = 0; i < n_points - 1; ++i)
4196  for (unsigned int j = 0; j < n_points - 1; ++j)
4197  {
4198  const Point<spacedim> p0 =
4199  mapping->transform_unit_to_real_cell(
4200  cell,
4201  q_projector->point(offset + i * n_points + j));
4202  out
4203  << p0 << ' ' << cell->level() << ' '
4204  << static_cast<unsigned int>(cell->material_id())
4205  << '\n';
4206  out
4207  << (mapping->transform_unit_to_real_cell(
4208  cell,
4209  q_projector->point(offset +
4210  (i + 1) * n_points + j)))
4211  << ' ' << cell->level() << ' '
4212  << static_cast<unsigned int>(cell->material_id())
4213  << '\n';
4214  out
4215  << (mapping->transform_unit_to_real_cell(
4216  cell,
4217  q_projector->point(
4218  offset + (i + 1) * n_points + j + 1)))
4219  << ' ' << cell->level() << ' '
4220  << static_cast<unsigned int>(cell->material_id())
4221  << '\n';
4222  out
4223  << (mapping->transform_unit_to_real_cell(
4224  cell,
4225  q_projector->point(offset + i * n_points +
4226  j + 1)))
4227  << ' ' << cell->level() << ' '
4228  << static_cast<unsigned int>(cell->material_id())
4229  << '\n';
4230  // and the first point again
4231  out
4232  << p0 << ' ' << cell->level() << ' '
4233  << static_cast<unsigned int>(cell->material_id())
4234  << '\n';
4235  out << '\n' << '\n';
4236  }
4237  }
4238  else
4239  {
4240  for (unsigned int l = 0;
4241  l < GeometryInfo<dim>::lines_per_face;
4242  ++l)
4243  {
4244  const typename ::Triangulation<dim, spacedim>::
4245  line_iterator line = face->line(l);
4246 
4247  const Point<spacedim> &v0 = line->vertex(0),
4248  &v1 = line->vertex(1);
4249  if (line->at_boundary() ||
4250  gnuplot_flags.curved_inner_cells)
4251  {
4252  // Save the points on each face to a vector and
4253  // then try to remove colinear points that won't
4254  // show up in the generated plot.
4255  std::vector<Point<spacedim>> line_points;
4256  // transform_real_to_unit_cell could be replaced
4257  // by using QProjector<dim>::project_to_line
4258  // which is not yet implemented
4259  const Point<spacedim>
4260  u0 = mapping->transform_real_to_unit_cell(cell,
4261  v0),
4262  u1 = mapping->transform_real_to_unit_cell(cell,
4263  v1);
4264  for (unsigned int i = 0; i < n_points; ++i)
4265  line_points.push_back(
4266  mapping->transform_unit_to_real_cell(
4267  cell,
4268  (1 - boundary_points[i][0]) * u0 +
4269  boundary_points[i][0] * u1));
4270  internal::remove_colinear_points(line_points);
4271  for (const Point<spacedim> &point : line_points)
4272  out << point << ' ' << cell->level() << ' '
4273  << static_cast<unsigned int>(
4274  cell->material_id())
4275  << '\n';
4276  }
4277  else
4278  out
4279  << v0 << ' ' << cell->level() << ' '
4280  << static_cast<unsigned int>(cell->material_id())
4281  << '\n'
4282  << v1 << ' ' << cell->level() << ' '
4283  << static_cast<unsigned int>(cell->material_id())
4284  << '\n';
4285 
4286  out << '\n' << '\n';
4287  }
4288  }
4289  }
4290  }
4291  }
4292 
4293  if (q_projector != nullptr)
4294  delete q_projector;
4295 
4296 
4297  // make sure everything now gets to disk
4298  out.flush();
4299 
4300  AssertThrow(out, ExcIO());
4301  }
4302  } // namespace
4303 } // namespace internal
4304 
4305 
4306 
4307 template <int dim, int spacedim>
4308 void
4310  std::ostream & out,
4311  const Mapping<dim, spacedim> * mapping) const
4312 {
4313  internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4314 }
4315 
4316 
4317 
4318 namespace internal
4319 {
4320  namespace
4321  {
4322  struct LineEntry
4323  {
4324  Point<2> first;
4325  Point<2> second;
4326  bool colorize;
4327  unsigned int level;
4328  LineEntry(const Point<2> & f,
4329  const Point<2> & s,
4330  const bool c,
4331  const unsigned int l)
4332  : first(f)
4333  , second(s)
4334  , colorize(c)
4335  , level(l)
4336  {}
4337  };
4338 
4339 
4340  void
4341  write_eps(const ::Triangulation<1> &,
4342  std::ostream &,
4343  const Mapping<1> *,
4344  const GridOutFlags::Eps<2> &,
4345  const GridOutFlags::Eps<3> &)
4346  {
4347  Assert(false, ExcNotImplemented());
4348  }
4349 
4350  void
4351  write_eps(const ::Triangulation<1, 2> &,
4352  std::ostream &,
4353  const Mapping<1, 2> *,
4354  const GridOutFlags::Eps<2> &,
4355  const GridOutFlags::Eps<3> &)
4356  {
4357  Assert(false, ExcNotImplemented());
4358  }
4359 
4360  void
4361  write_eps(const ::Triangulation<1, 3> &,
4362  std::ostream &,
4363  const Mapping<1, 3> *,
4364  const GridOutFlags::Eps<2> &,
4365  const GridOutFlags::Eps<3> &)
4366  {
4367  Assert(false, ExcNotImplemented());
4368  }
4369 
4370  void
4371  write_eps(const ::Triangulation<2, 3> &,
4372  std::ostream &,
4373  const Mapping<2, 3> *,
4374  const GridOutFlags::Eps<2> &,
4375  const GridOutFlags::Eps<3> &)
4376  {
4377  Assert(false, ExcNotImplemented());
4378  }
4379 
4380 
4381 
4382  template <int dim, int spacedim>
4383  void
4384  write_eps(const ::Triangulation<dim, spacedim> &tria,
4385  std::ostream & out,
4386  const Mapping<dim, spacedim> * mapping,
4389  {
4390  using LineList = std::list<LineEntry>;
4391 
4392  // We should never get here in 1D since this function is overloaded for
4393  // all dim == 1 cases.
4394  Assert(dim == 2 || dim == 3, ExcInternalError());
4395 
4396  // Copy, with an object slice, something containing the flags common to
4397  // all dimensions in order to avoid the recurring distinctions between
4398  // the different eps_flags present.
4399  const GridOutFlags::EpsFlagsBase eps_flags_base =
4400  dim == 2 ?
4401  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4402  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4403 
4404  AssertThrow(out, ExcIO());
4405  const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4406 
4407  // make up a list of lines by which
4408  // we will construct the triangulation
4409  //
4410  // this part unfortunately is a bit
4411  // dimension dependent, so we have to
4412  // treat every dimension different.
4413  // however, by directly producing
4414  // the lines to be printed, i.e. their
4415  // 2d images, we can later do the
4416  // actual output dimension independent
4417  // again
4418  LineList line_list;
4419 
4420  switch (dim)
4421  {
4422  case 1:
4423  {
4424  Assert(false, ExcInternalError());
4425  break;
4426  }
4427 
4428  case 2:
4429  {
4430  for (typename ::Triangulation<dim, spacedim>::
4431  active_cell_iterator cell = tria.begin_active();
4432  cell != tria.end();
4433  ++cell)
4434  for (unsigned int line_no = 0;
4435  line_no < GeometryInfo<dim>::lines_per_cell;
4436  ++line_no)
4437  {
4438  typename ::Triangulation<dim, spacedim>::line_iterator
4439  line = cell->line(line_no);
4440 
4441  // first treat all
4442  // interior lines and
4443  // make up a list of
4444  // them. if curved
4445  // lines shall not be
4446  // supported (i.e. no
4447  // mapping is
4448  // provided), then also
4449  // treat all other
4450  // lines
4451  if (!line->has_children() &&
4452  (mapping == nullptr || !line->at_boundary()))
4453  // one would expect
4454  // make_pair(line->vertex(0),
4455  // line->vertex(1))
4456  // here, but that is
4457  // not dimension
4458  // independent, since
4459  // vertex(i) is
4460  // Point<dim>, but we
4461  // want a Point<2>.
4462  // in fact, whenever
4463  // we're here, the
4464  // vertex is a
4465  // Point<dim>, but
4466  // the compiler does
4467  // not know
4468  // this. hopefully,
4469  // the compiler will
4470  // optimize away this
4471  // little kludge
4472  line_list.emplace_back(
4473  Point<2>(line->vertex(0)(0), line->vertex(0)(1)),
4474  Point<2>(line->vertex(1)(0), line->vertex(1)(1)),
4475  line->user_flag_set(),
4476  cell->level());
4477  }
4478 
4479  // next if we are to treat
4480  // curved boundaries
4481  // specially, then add lines
4482  // to the list consisting of
4483  // pieces of the boundary
4484  // lines
4485  if (mapping != nullptr)
4486  {
4487  // to do so, first
4488  // generate a sequence of
4489  // points on a face and
4490  // project them onto the
4491  // faces of a unit cell
4492  std::vector<Point<dim - 1>> boundary_points(n_points);
4493 
4494  for (unsigned int i = 0; i < n_points; ++i)
4495  boundary_points[i](0) = 1. * (i + 1) / (n_points + 1);
4496 
4497  Quadrature<dim - 1> quadrature(boundary_points);
4498  Quadrature<dim> q_projector(
4500 
4501  // next loop over all
4502  // boundary faces and
4503  // generate the info from
4504  // them
4505  for (typename ::Triangulation<dim, spacedim>::
4506  active_cell_iterator cell = tria.begin_active();
4507  cell != tria.end();
4508  ++cell)
4509  for (unsigned int face_no = 0;
4510  face_no < GeometryInfo<dim>::faces_per_cell;
4511  ++face_no)
4512  {
4513  const typename ::Triangulation<dim, spacedim>::
4514  face_iterator face = cell->face(face_no);
4515 
4516  if (face->at_boundary())
4517  {
4518  Point<dim> p0_dim(face->vertex(0));
4519  Point<2> p0(p0_dim(0), p0_dim(1));
4520 
4521  // loop over
4522  // all pieces
4523  // of the line
4524  // and generate
4525  // line-lets
4526  const unsigned int offset = face_no * n_points;
4527  for (unsigned int i = 0; i < n_points; ++i)
4528  {
4529  const Point<dim> p1_dim(
4530  mapping->transform_unit_to_real_cell(
4531  cell, q_projector.point(offset + i)));
4532  const Point<2> p1(p1_dim(0), p1_dim(1));
4533 
4534  line_list.emplace_back(p0,
4535  p1,
4536  face->user_flag_set(),
4537  cell->level());
4538  p0 = p1;
4539  }
4540 
4541  // generate last piece
4542  const Point<dim> p1_dim(face->vertex(1));
4543  const Point<2> p1(p1_dim(0), p1_dim(1));
4544  line_list.emplace_back(p0,
4545  p1,
4546  face->user_flag_set(),
4547  cell->level());
4548  }
4549  }
4550  }
4551 
4552  break;
4553  }
4554 
4555  case 3:
4556  {
4557  // curved boundary output
4558  // presently not supported
4559  Assert(mapping == nullptr, ExcNotImplemented());
4560 
4561  typename ::Triangulation<dim,
4562  spacedim>::active_cell_iterator
4563  cell = tria.begin_active(),
4564  endc = tria.end();
4565 
4566  // loop over all lines and compute their
4567  // projection on the plane perpendicular
4568  // to the direction of sight
4569 
4570  // direction of view equals the unit
4571  // vector of the position of the
4572  // spectator to the origin.
4573  //
4574  // we chose here the viewpoint as in
4575  // gnuplot as default.
4576  //
4577  // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4578  // GridOut
4579  // note: the following might be wrong
4580  // if one of the base vectors below
4581  // is in direction of the viewer, but
4582  // I am too tired at present to fix
4583  // this
4584  const double pi = numbers::PI;
4585  const double z_angle = eps_flags_3.azimut_angle;
4586  const double turn_angle = eps_flags_3.turn_angle;
4587  const Point<dim> view_direction(
4588  -std::sin(z_angle * 2. * pi / 360.) *
4589  std::sin(turn_angle * 2. * pi / 360.),
4590  +std::sin(z_angle * 2. * pi / 360.) *
4591  std::cos(turn_angle * 2. * pi / 360.),
4592  -std::cos(z_angle * 2. * pi / 360.));
4593 
4594  // decide about the two unit vectors
4595  // in this plane. we chose the first one
4596  // to be the projection of the z-axis
4597  // to this plane
4598  const Tensor<1, dim> vector1 =
4599  Point<dim>(0, 0, 1) -
4600  ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4601  const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4602 
4603  // now the third vector is fixed. we
4604  // chose the projection of a more or
4605  // less arbitrary vector to the plane
4606  // perpendicular to the first one
4607  const Tensor<1, dim> vector2 =
4608  (Point<dim>(1, 0, 0) -
4609  ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4610  ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4611  const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4612 
4613 
4614  for (; cell != endc; ++cell)
4615  for (unsigned int line_no = 0;
4616  line_no < GeometryInfo<dim>::lines_per_cell;
4617  ++line_no)
4618  {
4619  typename ::Triangulation<dim, spacedim>::line_iterator
4620  line = cell->line(line_no);
4621  line_list.emplace_back(
4622  Point<2>(line->vertex(0) * unit_vector2,
4623  line->vertex(0) * unit_vector1),
4624  Point<2>(line->vertex(1) * unit_vector2,
4625  line->vertex(1) * unit_vector1),
4626  line->user_flag_set(),
4627  cell->level());
4628  }
4629 
4630  break;
4631  }
4632 
4633  default:
4634  Assert(false, ExcNotImplemented());
4635  }
4636 
4637 
4638 
4639  // find out minimum and maximum x and
4640  // y coordinates to compute offsets
4641  // and scaling factors
4642  double x_min = tria.begin_active()->vertex(0)(0);
4643  double x_max = x_min;
4644  double y_min = tria.begin_active()->vertex(0)(1);
4645  double y_max = y_min;
4646  unsigned int max_level = line_list.begin()->level;
4647 
4648  for (LineList::const_iterator line = line_list.begin();
4649  line != line_list.end();
4650  ++line)
4651  {
4652  x_min = std::min(x_min, line->first(0));
4653  x_min = std::min(x_min, line->second(0));
4654 
4655  x_max = std::max(x_max, line->first(0));
4656  x_max = std::max(x_max, line->second(0));
4657 
4658  y_min = std::min(y_min, line->first(1));
4659  y_min = std::min(y_min, line->second(1));
4660 
4661  y_max = std::max(y_max, line->first(1));
4662  y_max = std::max(y_max, line->second(1));
4663 
4664  max_level = std::max(max_level, line->level);
4665  }
4666 
4667  // scale in x-direction such that
4668  // in the output 0 <= x <= 300.
4669  // don't scale in y-direction to
4670  // preserve the shape of the
4671  // triangulation
4672  const double scale =
4673  (eps_flags_base.size /
4674  (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4675  x_max - x_min :
4676  y_min - y_max));
4677 
4678 
4679  // now write preamble
4680  {
4681  // block this to have local
4682  // variables destroyed after
4683  // use
4684  std::time_t time1 = std::time(nullptr);
4685  std::tm * time = std::localtime(&time1);
4686  out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4687  << "%%Title: deal.II Output" << '\n'
4688  << "%%Creator: the deal.II library" << '\n'
4689  << "%%Creation Date: " << time->tm_year + 1900 << "/"
4690  << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4691  << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4692  << std::setw(2) << time->tm_sec << '\n'
4693  << "%%BoundingBox: "
4694  // lower left corner
4695  << "0 0 "
4696  // upper right corner
4697  << static_cast<unsigned int>(
4698  std::floor(((x_max - x_min) * scale) + 1))
4699  << ' '
4700  << static_cast<unsigned int>(
4701  std::floor(((y_max - y_min) * scale) + 1))
4702  << '\n';
4703 
4704  // define some abbreviations to keep
4705  // the output small:
4706  // m=move turtle to
4707  // x=execute line stroke
4708  // b=black pen
4709  // r=red pen
4710  out << "/m {moveto} bind def" << '\n'
4711  << "/x {lineto stroke} bind def" << '\n'
4712  << "/b {0 0 0 setrgbcolor} def" << '\n'
4713  << "/r {1 0 0 setrgbcolor} def" << '\n';
4714 
4715  // calculate colors for level
4716  // coloring; level 0 is black,
4717  // other levels are blue
4718  // ... red
4719  if (eps_flags_base.color_lines_level)
4720  out << "/l { neg " << (max_level) << " add "
4721  << (0.66666 / std::max(1U, (max_level - 1)))
4722  << " mul 1 0.8 sethsbcolor} def" << '\n';
4723 
4724  // in 2d, we can also plot cell
4725  // and vertex numbers, but this
4726  // requires a somewhat more
4727  // lengthy preamble. please
4728  // don't ask me what most of
4729  // this means, it is reverse
4730  // engineered from what GNUPLOT
4731  // uses in its output
4732  if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
4733  eps_flags_2.write_vertex_numbers))
4734  {
4735  out
4736  << ("/R {rmoveto} bind def\n"
4737  "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
4738  "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
4739  "currentdict end definefont\n"
4740  "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
4741  "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
4742  "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
4743  "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
4744  "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
4745  "5 get stringwidth pop add}\n"
4746  "{pop} ifelse} forall} bind def\n"
4747  "/MCshow { currentpoint stroke m\n"
4748  "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
4749  << '\n';
4750  }
4751 
4752  out << "%%EndProlog" << '\n' << '\n';
4753 
4754  // set fine lines
4755  out << eps_flags_base.line_width << " setlinewidth" << '\n';
4756  }
4757 
4758  // now write the lines
4759  const Point<2> offset(x_min, y_min);
4760 
4761  for (LineList::const_iterator line = line_list.begin();
4762  line != line_list.end();
4763  ++line)
4764  if (eps_flags_base.color_lines_level && (line->level > 0))
4765  out << line->level << " l " << (line->first - offset) * scale << " m "
4766  << (line->second - offset) * scale << " x" << '\n';
4767  else
4768  out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
4769  "r " :
4770  "b ")
4771  << (line->first - offset) * scale << " m "
4772  << (line->second - offset) * scale << " x" << '\n';
4773 
4774  // finally write the cell numbers
4775  // in 2d, if that is desired
4776  if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
4777  {
4778  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
4779 
4780  typename ::Triangulation<dim, spacedim>::active_cell_iterator
4781  cell = tria.begin_active(),
4782  endc = tria.end();
4783  for (; cell != endc; ++cell)
4784  {
4785  out << (cell->center()(0) - offset(0)) * scale << ' '
4786  << (cell->center()(1) - offset(1)) * scale << " m" << '\n'
4787  << "[ [(Helvetica) 12.0 0.0 true true (";
4788  if (eps_flags_2.write_cell_number_level)
4789  out << cell;
4790  else
4791  out << cell->index();
4792 
4793  out << ")] "
4794  << "] -6 MCshow" << '\n';
4795  }
4796  }
4797 
4798  // and the vertex numbers
4799  if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
4800  {
4801  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
4802 
4803  // have a list of those
4804  // vertices which we have
4805  // already tracked, to avoid
4806  // doing this multiply
4807  std::set<unsigned int> treated_vertices;
4808  typename ::Triangulation<dim, spacedim>::active_cell_iterator
4809  cell = tria.begin_active(),
4810  endc = tria.end();
4811  for (; cell != endc; ++cell)
4812  for (unsigned int vertex = 0;
4813  vertex < GeometryInfo<dim>::vertices_per_cell;
4814  ++vertex)
4815  if (treated_vertices.find(cell->vertex_index(vertex)) ==
4816  treated_vertices.end())
4817  {
4818  treated_vertices.insert(cell->vertex_index(vertex));
4819 
4820  out << (cell->vertex(vertex)(0) - offset(0)) * scale << ' '
4821  << (cell->vertex(vertex)(1) - offset(1)) * scale << " m"
4822  << '\n'
4823  << "[ [(Helvetica) 10.0 0.0 true true ("
4824  << cell->vertex_index(vertex) << ")] "
4825  << "] -6 MCshow" << '\n';
4826  }
4827  }
4828 
4829  out << "showpage" << '\n';
4830 
4831  // make sure everything now gets to
4832  // disk
4833  out.flush();
4834 
4835  AssertThrow(out, ExcIO());
4836  }
4837  } // namespace
4838 } // namespace internal
4839 
4840 
4841 template <int dim, int spacedim>
4842 void
4844  std::ostream & out,
4845  const Mapping<dim, spacedim> * mapping) const
4846 {
4847  internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
4848 }
4849 
4850 
4851 template <int dim, int spacedim>
4852 void
4854  std::ostream & out,
4855  const OutputFormat output_format,
4856  const Mapping<dim, spacedim> * mapping) const
4857 {
4858  switch (output_format)
4859  {
4860  case none:
4861  return;
4862 
4863  case dx:
4864  write_dx(tria, out);
4865  return;
4866 
4867  case ucd:
4868  write_ucd(tria, out);
4869  return;
4870 
4871  case gnuplot:
4872  write_gnuplot(tria, out, mapping);
4873  return;
4874 
4875  case eps:
4876  write_eps(tria, out, mapping);
4877  return;
4878 
4879  case xfig:
4880  write_xfig(tria, out, mapping);
4881  return;
4882 
4883  case msh:
4884  write_msh(tria, out);
4885  return;
4886 
4887  case svg:
4888  write_svg(tria, out);
4889  return;
4890 
4891  case mathgl:
4892  write_mathgl(tria, out);
4893  return;
4894 
4895  case vtk:
4896  write_vtk(tria, out);
4897  return;
4898 
4899  case vtu:
4900  write_vtu(tria, out);
4901  return;
4902  }
4903 
4904  Assert(false, ExcInternalError());
4905 }
4906 
4907 
4908 template <int dim, int spacedim>
4909 void
4911  std::ostream & out,
4912  const Mapping<dim, spacedim> * mapping) const
4913 {
4914  write(tria, out, default_format, mapping);
4915 }
4916 
4917 
4918 // explicit instantiations
4919 #include "grid_out.inst"
4920 
4921 
4922 DEAL_II_NAMESPACE_CLOSE
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:193
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:716
unsigned int n_boundary_faces(const Triangulation< dim, spacedim > &tria) const
Definition: grid_out.cc:3510
unsigned int n_active_cells() const
Definition: tria.cc:12605
long int get_integer(const std::string &entry_string) const
unsigned int n_used_vertices() const
Definition: tria.cc:13189
const types::manifold_id flat_manifold_id
Definition: types.h:267
static const unsigned int invalid_unsigned_int
Definition: types.h:187
OutputFormat default_format
Definition: grid_out.h:1503
static void declare_parameters(ParameterHandler &prm)
DX(const bool write_cells=true, const bool write_faces=false, const bool write_diameter=false, const bool write_measure=false, const bool write_all_faces=true)
Definition: grid_out.cc:47
active_face_iterator begin_active_face() const
Definition: tria.cc:12172
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="")
static OutputFormat parse_output_format(const std::string &format_name)
Definition: grid_out.cc:617
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:375
bool write_cells
Definition: grid_out.h:61
write() calls write_dx()
Definition: grid_out.h:1005
void attach_triangulation(const Triangulation< DoFHandlerType::dimension, DoFHandlerType::space_dimension > &)
bool write_additional_boundary_lines
Definition: grid_out.h:284
void write_xfig(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:1316
static ::ExceptionBase & ExcIO()
unsigned int n_boundary_face_points
Definition: grid_out.h:607
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:467
GridOutFlags::Eps< 2 > eps_flags_2
Definition: grid_out.h:1538
unsigned int height
Definition: grid_out.h:679
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:12115
OutputFormat
Definition: grid_out.h:1000
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3298
std::string get(const std::string &entry_string) const
GridOut()
Definition: grid_out.cc:481
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:249
void write_mesh_per_processor_as_vtu(const Triangulation< dim, spacedim > &tria, const std::string &filename_without_extension, const bool view_levels=false, const bool include_artificial=false) const
Definition: grid_out.cc:3328
const Point< dim > & point(const unsigned int i) const
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:399
Convert the level number into the cell color.
Definition: grid_out.h:751
write() calls write_eps()
Definition: grid_out.h:1009
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: tria.cc:11943
unsigned int line_thickness
Definition: grid_out.h:690
bool convert_level_number_to_height
Definition: grid_out.h:762
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
bool write_diameter
Definition: grid_out.h:71
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1426
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:2752
virtual Point< dim > transform_real_to_unit_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p) const =0
unsigned int boundary_line_thickness
Definition: grid_out.h:694
Point< spacedim > vertices[GeometryInfo< dim >::vertices_per_cell]
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:11923
Svg(const unsigned int line_thickness=2, const unsigned int boundary_line_thickness=4, const bool margin=true, const Background background=white, const int azimuth_angle=0, const int polar_angle=0, const Coloring coloring=level_number, const bool convert_level_number_to_height=false, const bool label_level_number=false, const bool label_cell_index=false, const bool label_material_id=false, const bool label_subdomain_id=false, const bool draw_colorbar=false, const bool draw_legend=false)
Definition: grid_out.cc:427
unsigned int write_msh_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3633
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:88
__global__ void scale(Number *val, const Number *V_val, const size_type N)
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:332
unsigned int n_extra_curved_line_points
Definition: grid_out.h:253
bool write_all_faces
Definition: grid_out.h:82
write() calls write_ucd()
Definition: grid_out.h:1011
bool output_only_relevant
Definition: grid_out.h:901
void write_gnuplot(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:4309
cell_iterator end() const
Definition: tria.cc:12009
unsigned int n_boundary_lines(const Triangulation< dim, spacedim > &tria) const
Definition: grid_out.cc:3527
void enter_subsection(const std::string &subsection)
GridOutFlags::Vtk vtk_flags
Definition: grid_out.h:1564
Point< 2 > scaling
Definition: grid_out.h:612
bool label_cell_index
Definition: grid_out.h:782
GridOutFlags::Gnuplot gnuplot_flags
Definition: grid_out.h:1526
void write_svg(const Triangulation< 2, 2 > &tria, std::ostream &out) const
Definition: grid_out.cc:1567
Convert the global subdomain id into the cell color.
Definition: grid_out.h:591
unsigned int & n_boundary_face_points
Definition: grid_out.h:265
Convert the material id into the cell color.
Definition: grid_out.h:587
static ::ExceptionBase & ExcMessage(std::string arg1)
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:185
std::size_t memory_consumption() const
Definition: grid_out.cc:762
double get_double(const std::string &entry_name) const
GridOutFlags::MathGL mathgl_flags
Definition: grid_out.h:1559
void save_user_flags_line(std::ostream &out) const
Definition: tria.cc:11234
write() calls write_mathgl()
Definition: grid_out.h:1019
GridOutFlags::DX dx_flags
Definition: grid_out.h:1508
GridOutFlags::Msh msh_flags
Definition: grid_out.h:1514
GridOutFlags::Eps< 1 > eps_flags_1
Definition: grid_out.h:1532
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
write() calls write_gnuplot()
Definition: grid_out.h:1007
#define Assert(cond, exc)
Definition: exceptions.h:1411
void parse_parameters(const ParameterHandler &prm)
bool label_material_id
Definition: grid_out.h:787
void write_vtu(const std::vector< Patch< dim, spacedim >> &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string >> &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
EpsFlagsBase(const SizeType size_type=width, const unsigned int size=300, const double line_width=0.5, const bool color_lines_on_user_flag=false, const unsigned int n_boundary_face_points=2, const bool color_lines_level=false)
Definition: grid_out.cc:200
Abstract base class for mapping classes.
Definition: mapping.h:302
const types::boundary_id invalid_boundary_id
Definition: types.h:228
GridOutFlags::XFig xfig_flags
Definition: grid_out.h:1549
bool label_subdomain_id
Definition: grid_out.h:792
Ucd(const bool write_preamble=false, const bool write_faces=false, const bool write_lines=false)
Definition: grid_out.cc:119
const std::vector< Point< spacedim > > & get_vertices() const
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:667
bool label_level_subdomain_id
Definition: grid_out.h:797
static std::string get_output_format_names()
Definition: grid_out.cc:660
Convert the subdomain id into the cell color.
Definition: grid_out.h:753
Msh(const bool write_faces=false, const bool write_lines=false)
Definition: grid_out.cc:98
unsigned int n_subdivisions
bool get_bool(const std::string &entry_name) const
void write_vtk(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3086
IteratorRange< active_face_iterator > active_face_iterators() const
Definition: tria.cc:12214
GridOutFlags::Ucd ucd_flags
Definition: grid_out.h:1520
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:130
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:383
Background background
Definition: grid_out.h:726
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:216
Do nothing in write()
Definition: grid_out.h:1003
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:75
float cell_font_scaling
Definition: grid_out.h:773
Convert the level subdomain id into the cell color.
Definition: grid_out.h:593
write() calls write_xfig()
Definition: grid_out.h:1013
Convert the level subdomain id into the cell color.
Definition: grid_out.h:755
static Quadrature< dim > project_to_all_faces(const SubQuadrature &quadrature)
bool write_faces
Definition: grid_out.h:66
const types::subdomain_id artificial_subdomain_id
Definition: types.h:296
Convert the level into the cell color.
Definition: grid_out.h:589
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:112
unsigned int n_boundary_face_points
Definition: grid_out.h:387
Point< 2 > offset
Definition: grid_out.h:618
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:104
Definition: tensor.h:422
unsigned int write_ucd_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3798
write() calls write_msh()
Definition: grid_out.h:1015
static constexpr double PI
Definition: numbers.h:238
const std::vector< bool > & get_used_vertices() const
Definition: tria.cc:13198
void write_ucd(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:1191
unsigned int write_ucd_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3841
write() calls write_svg()
Definition: grid_out.h:1017
unsigned int write_msh_lines(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3677
GridOutFlags::Svg svg_flags
Definition: grid_out.h:1554
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:415
std::array< unsigned int, GeometryInfo< dim >::faces_per_cell > neighbors
void write_eps(const Triangulation< dim, spacedim > &tria, std::ostream &out, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:4843
bool write_measure
Definition: grid_out.h:76
static ::ExceptionBase & ExcNotImplemented()
Convert the material id into the cell color (default)
Definition: grid_out.h:749
IteratorRange< cell_iterator > cell_iterators() const
Definition: tria.cc:12106
float level_height_factor
Definition: grid_out.h:768
GridOutFlags::Eps< 3 > eps_flags_3
Definition: grid_out.h:1544
face_iterator end_face() const
Definition: tria.cc:12193
static ::ExceptionBase & ExcInvalidState()
const types::boundary_id internal_face_boundary_id
Definition: types.h:244
write() calls write_vtu()
Definition: grid_out.h:1023
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:60
void set_flags(const GridOutFlags::DX &flags)
Definition: grid_out.cc:487
bool label_level_number
Definition: grid_out.h:777
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:139
void write_dx(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:797
void write(const Triangulation< dim, spacedim > &tria, std::ostream &out, const OutputFormat output_format, const Mapping< dim, spacedim > *mapping=nullptr) const
Definition: grid_out.cc:4853
virtual Point< spacedim > transform_unit_to_real_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &p) const =0
unsigned int boundary_id
Definition: types.h:125
Table< 2, float > data
Gnuplot & operator=(const Gnuplot &flags)
Definition: grid_out.cc:172
GridOutFlags::Vtu vtu_flags
Definition: grid_out.h:1569
unsigned int width
Definition: grid_out.h:685
write() calls write_vtk()
Definition: grid_out.h:1021
virtual types::subdomain_id locally_owned_subdomain() const
Definition: tria.cc:13267
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:473
void parse_parameters(ParameterHandler &param)
Definition: grid_out.cc:283
Gnuplot(const bool write_cell_number=false, const unsigned int n_extra_curved_line_points=2, const bool curved_inner_cells=false, const bool write_additional_boundary_lines=true)
Definition: grid_out.cc:148
static ::ExceptionBase & ExcInternalError()
void write_msh(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:1044
std::string default_suffix() const
Definition: grid_out.cc:609