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