Reference documentation for deal.II version Git 9297d75edf 2020-11-26 18:52:14 +0100
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 unsigned int 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 (unsigned int 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 < cell->n_vertices();
1588  ++vertex_index)
1589  {
1590  if (cell->vertex(vertex_index)[0] < x_min)
1591  x_min = cell->vertex(vertex_index)[0];
1592  if (cell->vertex(vertex_index)[0] > x_max)
1593  x_max = cell->vertex(vertex_index)[0];
1594 
1595  if (cell->vertex(vertex_index)[1] < y_min)
1596  y_min = cell->vertex(vertex_index)[1];
1597  if (cell->vertex(vertex_index)[1] > y_max)
1598  y_max = cell->vertex(vertex_index)[1];
1599  }
1600 
1601  if (static_cast<unsigned int>(cell->level()) < min_level)
1602  min_level = cell->level();
1603  if (static_cast<unsigned int>(cell->level()) > max_level)
1604  max_level = cell->level();
1605 
1606  materials.insert(cell->material_id());
1607  levels.insert(cell->level());
1608  if (cell->is_active())
1609  subdomains.insert(cell->subdomain_id() + 2);
1610  level_subdomains.insert(cell->level_subdomain_id() + 2);
1611  }
1612 
1613  x_dimension = x_max - x_min;
1614  y_dimension = y_max - y_min;
1615 
1616  // count the materials being used
1617  const unsigned int n_materials = materials.size();
1618 
1619  // count the levels being used
1620  const unsigned int n_levels = levels.size();
1621 
1622  // count the subdomains being used
1623  const unsigned int n_subdomains = subdomains.size();
1624 
1625  // count the level subdomains being used
1626  const unsigned int n_level_subdomains = level_subdomains.size();
1627 
1628  switch (svg_flags.coloring)
1629  {
1631  n = n_materials;
1632  break;
1634  n = n_levels;
1635  break;
1637  n = n_subdomains;
1638  break;
1640  n = n_level_subdomains;
1641  break;
1642  default:
1643  break;
1644  }
1645 
1646  // set the camera position to top view, targeting at the origin
1647  // vectors and variables for the perspective view
1648  Point<3> camera_position;
1649  camera_position[0] = 0;
1650  camera_position[1] = 0;
1651  camera_position[2] = 2. * std::max(x_dimension, y_dimension);
1652 
1653  Tensor<1, 3> camera_direction;
1654  camera_direction[0] = 0;
1655  camera_direction[1] = 0;
1656  camera_direction[2] = -1;
1657 
1658  Tensor<1, 3> camera_horizontal;
1659  camera_horizontal[0] = 1;
1660  camera_horizontal[1] = 0;
1661  camera_horizontal[2] = 0;
1662 
1663  camera_focus = .5 * std::max(x_dimension, y_dimension);
1664 
1665  Point<3> camera_position_temp;
1666  Point<3> camera_direction_temp;
1667  Point<3> camera_horizontal_temp;
1668 
1669  const double angle_factor = 3.14159265 / 180.;
1670 
1671  // (I) rotate the camera to the chosen polar angle
1672  camera_position_temp[1] =
1673  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[1] -
1674  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[2];
1675  camera_position_temp[2] =
1676  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[1] +
1677  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[2];
1678 
1679  camera_direction_temp[1] =
1680  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[1] -
1681  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1682  camera_direction_temp[2] =
1683  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[1] +
1684  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1685 
1686  camera_horizontal_temp[1] =
1687  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] -
1688  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1689  camera_horizontal_temp[2] =
1690  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] +
1691  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1692 
1693  camera_position[1] = camera_position_temp[1];
1694  camera_position[2] = camera_position_temp[2];
1695 
1696  camera_direction[1] = camera_direction_temp[1];
1697  camera_direction[2] = camera_direction_temp[2];
1698 
1699  camera_horizontal[1] = camera_horizontal_temp[1];
1700  camera_horizontal[2] = camera_horizontal_temp[2];
1701 
1702  // (II) rotate the camera to the chosen azimuth angle
1703  camera_position_temp[0] =
1704  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[0] -
1705  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1706  camera_position_temp[1] =
1707  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[0] +
1708  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1709 
1710  camera_direction_temp[0] =
1711  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] -
1712  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1713  camera_direction_temp[1] =
1714  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] +
1715  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1716 
1717  camera_horizontal_temp[0] =
1718  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] -
1719  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1720  camera_horizontal_temp[1] =
1721  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] +
1722  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1723 
1724  camera_position[0] = camera_position_temp[0];
1725  camera_position[1] = camera_position_temp[1];
1726 
1727  camera_direction[0] = camera_direction_temp[0];
1728  camera_direction[1] = camera_direction_temp[1];
1729 
1730  camera_horizontal[0] = camera_horizontal_temp[0];
1731  camera_horizontal[1] = camera_horizontal_temp[1];
1732 
1733  // translate the camera to the given triangulation
1734  camera_position[0] = x_min + .5 * x_dimension;
1735  camera_position[1] = y_min + .5 * y_dimension;
1736 
1737  camera_position[0] += 2. * std::max(x_dimension, y_dimension) *
1738  std::sin(angle_factor * svg_flags.polar_angle) *
1739  std::sin(angle_factor * svg_flags.azimuth_angle);
1740  camera_position[1] -= 2. * std::max(x_dimension, y_dimension) *
1741  std::sin(angle_factor * svg_flags.polar_angle) *
1742  std::cos(angle_factor * svg_flags.azimuth_angle);
1743 
1744 
1745  // determine the bounding box of the given triangulation on the projection
1746  // plane of the camera viewing system
1747  point[0] = tria.begin()->vertex(0)[0];
1748  point[1] = tria.begin()->vertex(0)[1];
1749  point[2] = 0;
1750 
1751  float min_level_min_vertex_distance = 0;
1752 
1754  {
1755  point[2] = svg_flags.level_height_factor *
1756  (static_cast<float>(tria.begin()->level()) /
1757  static_cast<float>(n_levels)) *
1758  std::max(x_dimension, y_dimension);
1759  }
1760 
1761  projection_decomposition = svg_project_point(
1762  point, camera_position, camera_direction, camera_horizontal, camera_focus);
1763 
1764  x_max_perspective = projection_decomposition[0];
1765  x_min_perspective = projection_decomposition[0];
1766 
1767  y_max_perspective = projection_decomposition[1];
1768  y_min_perspective = projection_decomposition[1];
1769 
1770  for (const auto &cell : tria.cell_iterators())
1771  {
1772  point[0] = cell->vertex(0)[0];
1773  point[1] = cell->vertex(0)[1];
1774  point[2] = 0;
1775 
1777  {
1778  point[2] =
1780  (static_cast<float>(cell->level()) / static_cast<float>(n_levels)) *
1781  std::max(x_dimension, y_dimension);
1782  }
1783 
1784  projection_decomposition = svg_project_point(point,
1785  camera_position,
1786  camera_direction,
1787  camera_horizontal,
1788  camera_focus);
1789 
1790  if (x_max_perspective < projection_decomposition[0])
1791  x_max_perspective = projection_decomposition[0];
1792  if (x_min_perspective > projection_decomposition[0])
1793  x_min_perspective = projection_decomposition[0];
1794 
1795  if (y_max_perspective < projection_decomposition[1])
1796  y_max_perspective = projection_decomposition[1];
1797  if (y_min_perspective > projection_decomposition[1])
1798  y_min_perspective = projection_decomposition[1];
1799 
1800  point[0] = cell->vertex(1)[0];
1801  point[1] = cell->vertex(1)[1];
1802 
1803  projection_decomposition = svg_project_point(point,
1804  camera_position,
1805  camera_direction,
1806  camera_horizontal,
1807  camera_focus);
1808 
1809  if (x_max_perspective < projection_decomposition[0])
1810  x_max_perspective = projection_decomposition[0];
1811  if (x_min_perspective > projection_decomposition[0])
1812  x_min_perspective = projection_decomposition[0];
1813 
1814  if (y_max_perspective < projection_decomposition[1])
1815  y_max_perspective = projection_decomposition[1];
1816  if (y_min_perspective > projection_decomposition[1])
1817  y_min_perspective = projection_decomposition[1];
1818 
1819  point[0] = cell->vertex(2)[0];
1820  point[1] = cell->vertex(2)[1];
1821 
1822  projection_decomposition = svg_project_point(point,
1823  camera_position,
1824  camera_direction,
1825  camera_horizontal,
1826  camera_focus);
1827 
1828  if (x_max_perspective < projection_decomposition[0])
1829  x_max_perspective = projection_decomposition[0];
1830  if (x_min_perspective > projection_decomposition[0])
1831  x_min_perspective = projection_decomposition[0];
1832 
1833  if (y_max_perspective < projection_decomposition[1])
1834  y_max_perspective = projection_decomposition[1];
1835  if (y_min_perspective > projection_decomposition[1])
1836  y_min_perspective = projection_decomposition[1];
1837 
1838  if (cell->n_vertices() == 4) // in case of quadrilateral
1839  {
1840  point[0] = cell->vertex(3)[0];
1841  point[1] = cell->vertex(3)[1];
1842 
1843  projection_decomposition = svg_project_point(point,
1844  camera_position,
1845  camera_direction,
1846  camera_horizontal,
1847  camera_focus);
1848 
1849  if (x_max_perspective < projection_decomposition[0])
1850  x_max_perspective = projection_decomposition[0];
1851  if (x_min_perspective > projection_decomposition[0])
1852  x_min_perspective = projection_decomposition[0];
1853 
1854  if (y_max_perspective < projection_decomposition[1])
1855  y_max_perspective = projection_decomposition[1];
1856  if (y_min_perspective > projection_decomposition[1])
1857  y_min_perspective = projection_decomposition[1];
1858  }
1859 
1860  if (static_cast<unsigned int>(cell->level()) == min_level)
1861  min_level_min_vertex_distance = cell->minimum_vertex_distance();
1862  }
1863 
1864  x_dimension_perspective = x_max_perspective - x_min_perspective;
1865  y_dimension_perspective = y_max_perspective - y_min_perspective;
1866 
1867  // create the svg file with an internal style sheet
1868  if (width == 0)
1869  width = static_cast<unsigned int>(
1870  .5 + height * (x_dimension_perspective / y_dimension_perspective));
1871  else if (height == 0)
1872  height = static_cast<unsigned int>(
1873  .5 + width * (y_dimension_perspective / x_dimension_perspective));
1874  unsigned int additional_width = 0;
1875  // font size for date, time, legend, and colorbar
1876  unsigned int font_size =
1877  static_cast<unsigned int>(.5 + (height / 100.) * 1.75);
1878  cell_label_font_size = static_cast<unsigned int>(
1879  .5 + (height * .15 * svg_flags.cell_font_scaling *
1880  min_level_min_vertex_distance / std::min(x_dimension, y_dimension)));
1881 
1882  if (svg_flags.draw_legend &&
1886  {
1887  additional_width = static_cast<unsigned int>(
1888  .5 + height * .4); // additional width for legend
1889  }
1891  {
1892  additional_width = static_cast<unsigned int>(
1893  .5 + height * .175); // additional width for colorbar
1894  }
1895 
1896  // out << "<!-- deal.ii GridOut " << now->tm_mday << '/' << now->tm_mon + 1 <<
1897  // '/' << now->tm_year + 1900
1898  // << ' ' << now->tm_hour << ':';
1899  //
1900  // if (now->tm_min < 10) out << '0';
1901  //
1902  // out << now->tm_min << " -->" << '\n';
1903 
1904  // basic svg header
1905  out << "<svg width=\"" << width + additional_width << "\" height=\"" << height
1906  << "\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" << '\n'
1907  << '\n';
1908 
1909 
1911  {
1912  out
1913  << " <linearGradient id=\"background_gradient\" gradientUnits=\"userSpaceOnUse\" x1=\"0\" y1=\"0\" x2=\"0\" y2=\""
1914  << height << "\">" << '\n'
1915  << " <stop offset=\"0\" style=\"stop-color:white\"/>" << '\n'
1916  << " <stop offset=\"1\" style=\"stop-color:lightsteelblue\"/>" << '\n'
1917  << " </linearGradient>" << '\n';
1918  }
1919 
1920  out << '\n';
1921 
1922  // header for the internal style sheet
1923  out << "<!-- internal style sheet -->" << '\n'
1924  << "<style type=\"text/css\"><![CDATA[" << '\n';
1925 
1926  // set the background of the output graphic
1928  out << " rect.background{fill:url(#background_gradient)}" << '\n';
1930  out << " rect.background{fill:white}" << '\n';
1931  else
1932  out << " rect.background{fill:none}" << '\n';
1933 
1934  // basic svg graphic element styles
1935  out << " rect{fill:none; stroke:rgb(25,25,25); stroke-width:"
1936  << svg_flags.line_thickness << '}' << '\n'
1937  << " text{font-family:Helvetica; text-anchor:middle; fill:rgb(25,25,25)}"
1938  << '\n'
1939  << " line{stroke:rgb(25,25,25); stroke-width:"
1940  << svg_flags.boundary_line_thickness << '}' << '\n'
1941  << " path{fill:none; stroke:rgb(25,25,25); stroke-width:"
1942  << svg_flags.line_thickness << '}' << '\n'
1943  << " circle{fill:white; stroke:black; stroke-width:2}" << '\n'
1944  << '\n';
1945 
1946  // polygon styles with respect to the chosen cell coloring
1947  if (svg_flags.coloring)
1948  {
1949  unsigned int labeling_index = 0;
1950  auto materials_it = materials.begin();
1951  auto levels_it = levels.begin();
1952  auto subdomains_it = subdomains.begin();
1953  auto level_subdomains_it = level_subdomains.begin();
1954 
1955  for (unsigned int index = 0; index < n; index++)
1956  {
1957  double h;
1958 
1959  if (n != 1)
1960  h = .6 - (index / (n - 1.)) * .6;
1961  else
1962  h = .6;
1963 
1964  unsigned int r = 0;
1965  unsigned int g = 0;
1966  unsigned int b = 0;
1967 
1968  unsigned int i = static_cast<unsigned int>(h * 6);
1969 
1970  double f = h * 6 - i;
1971  double q = 1 - f;
1972  double t = f;
1973 
1974  switch (i % 6)
1975  {
1976  case 0:
1977  r = 255, g = static_cast<unsigned int>(.5 + 255 * t);
1978  break;
1979  case 1:
1980  r = static_cast<unsigned int>(.5 + 255 * q), g = 255;
1981  break;
1982  case 2:
1983  g = 255, b = static_cast<unsigned int>(.5 + 255 * t);
1984  break;
1985  case 3:
1986  g = static_cast<unsigned int>(.5 + 255 * q), b = 255;
1987  break;
1988  case 4:
1989  r = static_cast<unsigned int>(.5 + 255 * t), b = 255;
1990  break;
1991  case 5:
1992  r = 255, b = static_cast<unsigned int>(.5 + 255 * q);
1993  break;
1994  default:
1995  break;
1996  }
1997 
1998  switch (svg_flags.coloring)
1999  {
2001  labeling_index = *materials_it++;
2002  break;
2004  labeling_index = *levels_it++;
2005  break;
2007  labeling_index = *subdomains_it++;
2008  break;
2010  labeling_index = *level_subdomains_it++;
2011  break;
2012  default:
2013  break;
2014  }
2015 
2016  out << " path.p" << labeling_index << "{fill:rgb(" << r << ',' << g
2017  << ',' << b << "); "
2018  << "stroke:rgb(25,25,25); stroke-width:"
2019  << svg_flags.line_thickness << '}' << '\n';
2020 
2021  out << " path.ps" << labeling_index << "{fill:rgb("
2022  << static_cast<unsigned int>(.5 + .75 * r) << ','
2023  << static_cast<unsigned int>(.5 + .75 * g) << ','
2024  << static_cast<unsigned int>(.5 + .75 * b) << "); "
2025  << "stroke:rgb(20,20,20); stroke-width:"
2026  << svg_flags.line_thickness << '}' << '\n';
2027 
2028  out << " rect.r" << labeling_index << "{fill:rgb(" << r << ',' << g
2029  << ',' << b << "); "
2030  << "stroke:rgb(25,25,25); stroke-width:"
2031  << svg_flags.line_thickness << '}' << '\n';
2032 
2033  labeling_index++;
2034  }
2035  }
2036 
2037  out << "]]></style>" << '\n' << '\n';
2038 
2039  // background rectangle
2040  out << " <rect class=\"background\" width=\"" << width << "\" height=\""
2041  << height << "\"/>" << '\n';
2042 
2044  {
2045  unsigned int x_offset = 0;
2046 
2047  if (svg_flags.margin)
2048  x_offset = static_cast<unsigned int>(.5 + (height / 100.) *
2049  (margin_in_percent / 2.));
2050  else
2051  x_offset = static_cast<unsigned int>(.5 + height * .025);
2052 
2053  out
2054  << " <text x=\"" << x_offset << "\" y=\""
2055  << static_cast<unsigned int>(.5 + height * .0525) << '\"'
2056  << " style=\"font-weight:100; fill:lightsteelblue; text-anchor:start; font-family:Courier; font-size:"
2057  << static_cast<unsigned int>(.5 + height * .045) << "px\">"
2058  << "deal.II"
2059  << "</text>" << '\n';
2060 
2061  // out << " <text x=\"" << x_offset + static_cast<unsigned int>(.5 +
2062  // height * .045 * 4.75) << "\" y=\"" << static_cast<unsigned int>(.5 +
2063  // height * .0525) << '\"'
2064  // << " style=\"fill:lightsteelblue; text-anchor:start; font-size:" <<
2065  // font_size << "\">"
2066  // << now->tm_mday << '/' << now->tm_mon + 1 << '/' << now->tm_year +
2067  // 1900
2068  // << " - " << now->tm_hour << ':';
2069  //
2070  // if(now->tm_min < 10) out << '0';
2071  //
2072  // out << now->tm_min
2073  // << "</text>"<< '\n' << '\n';
2074  }
2075 
2076  // draw the cells, starting out from the minimal level (in order to guaranty a
2077  // correct perspective view)
2078  out << " <!-- cells -->" << '\n';
2079 
2080  for (unsigned int level_index = min_level; level_index <= max_level;
2081  level_index++)
2082  {
2083  for (const auto &cell : tria.cell_iterators_on_level(level_index))
2084  {
2085  if (!svg_flags.convert_level_number_to_height && !cell->is_active())
2086  continue;
2087 
2088  // draw the current cell
2089  out << " <path";
2090 
2091  if (svg_flags.coloring)
2092  {
2093  out << " class=\"p";
2094 
2095  if (!cell->is_active() &&
2097  out << 's';
2098 
2099  switch (svg_flags.coloring)
2100  {
2102  out << cell->material_id();
2103  break;
2105  out << static_cast<unsigned int>(cell->level());
2106  break;
2108  if (cell->is_active())
2109  out << cell->subdomain_id() + 2;
2110  else
2111  out << 'X';
2112  break;
2114  out << cell->level_subdomain_id() + 2;
2115  break;
2116  default:
2117  break;
2118  }
2119 
2120  out << '\"';
2121  }
2122 
2123  out << " d=\"M ";
2124 
2125  point[0] = cell->vertex(0)[0];
2126  point[1] = cell->vertex(0)[1];
2127  point[2] = 0;
2128 
2130  {
2131  point[2] = svg_flags.level_height_factor *
2132  (static_cast<float>(cell->level()) /
2133  static_cast<float>(n_levels)) *
2134  std::max(x_dimension, y_dimension);
2135  }
2136 
2137  projection_decomposition = svg_project_point(point,
2138  camera_position,
2139  camera_direction,
2140  camera_horizontal,
2141  camera_focus);
2142 
2143  out << static_cast<unsigned int>(
2144  .5 +
2145  ((projection_decomposition[0] - x_min_perspective) /
2146  x_dimension_perspective) *
2147  (width - (width / 100.) * 2. * margin_in_percent) +
2148  ((width / 100.) * margin_in_percent))
2149  << ' '
2150  << static_cast<unsigned int>(
2151  .5 + height - (height / 100.) * margin_in_percent -
2152  ((projection_decomposition[1] - y_min_perspective) /
2153  y_dimension_perspective) *
2154  (height - (height / 100.) * 2. * margin_in_percent));
2155 
2156  out << " L ";
2157 
2158  point[0] = cell->vertex(1)[0];
2159  point[1] = cell->vertex(1)[1];
2160 
2161  projection_decomposition = svg_project_point(point,
2162  camera_position,
2163  camera_direction,
2164  camera_horizontal,
2165  camera_focus);
2166 
2167  out << static_cast<unsigned int>(
2168  .5 +
2169  ((projection_decomposition[0] - x_min_perspective) /
2170  x_dimension_perspective) *
2171  (width - (width / 100.) * 2. * margin_in_percent) +
2172  ((width / 100.) * margin_in_percent))
2173  << ' '
2174  << static_cast<unsigned int>(
2175  .5 + height - (height / 100.) * margin_in_percent -
2176  ((projection_decomposition[1] - y_min_perspective) /
2177  y_dimension_perspective) *
2178  (height - (height / 100.) * 2. * margin_in_percent));
2179 
2180  out << " L ";
2181 
2182  if (cell->n_vertices() == 4) // in case of quadrilateral
2183  {
2184  point[0] = cell->vertex(3)[0];
2185  point[1] = cell->vertex(3)[1];
2186 
2187  projection_decomposition = svg_project_point(point,
2188  camera_position,
2189  camera_direction,
2190  camera_horizontal,
2191  camera_focus);
2192 
2193  out << static_cast<unsigned int>(
2194  .5 +
2195  ((projection_decomposition[0] - x_min_perspective) /
2196  x_dimension_perspective) *
2197  (width - (width / 100.) * 2. * margin_in_percent) +
2198  ((width / 100.) * margin_in_percent))
2199  << ' '
2200  << static_cast<unsigned int>(
2201  .5 + height - (height / 100.) * margin_in_percent -
2202  ((projection_decomposition[1] - y_min_perspective) /
2203  y_dimension_perspective) *
2204  (height - (height / 100.) * 2. * margin_in_percent));
2205 
2206  out << " L ";
2207  }
2208 
2209  point[0] = cell->vertex(2)[0];
2210  point[1] = cell->vertex(2)[1];
2211 
2212  projection_decomposition = svg_project_point(point,
2213  camera_position,
2214  camera_direction,
2215  camera_horizontal,
2216  camera_focus);
2217 
2218  out << static_cast<unsigned int>(
2219  .5 +
2220  ((projection_decomposition[0] - x_min_perspective) /
2221  x_dimension_perspective) *
2222  (width - (width / 100.) * 2. * margin_in_percent) +
2223  ((width / 100.) * margin_in_percent))
2224  << ' '
2225  << static_cast<unsigned int>(
2226  .5 + height - (height / 100.) * margin_in_percent -
2227  ((projection_decomposition[1] - y_min_perspective) /
2228  y_dimension_perspective) *
2229  (height - (height / 100.) * 2. * margin_in_percent));
2230 
2231  out << " L ";
2232 
2233  point[0] = cell->vertex(0)[0];
2234  point[1] = cell->vertex(0)[1];
2235 
2236  projection_decomposition = svg_project_point(point,
2237  camera_position,
2238  camera_direction,
2239  camera_horizontal,
2240  camera_focus);
2241 
2242  out << static_cast<unsigned int>(
2243  .5 +
2244  ((projection_decomposition[0] - x_min_perspective) /
2245  x_dimension_perspective) *
2246  (width - (width / 100.) * 2. * margin_in_percent) +
2247  ((width / 100.) * margin_in_percent))
2248  << ' '
2249  << static_cast<unsigned int>(
2250  .5 + height - (height / 100.) * margin_in_percent -
2251  ((projection_decomposition[1] - y_min_perspective) /
2252  y_dimension_perspective) *
2253  (height - (height / 100.) * 2. * margin_in_percent));
2254 
2255  out << "\"/>" << '\n';
2256 
2257  // label the current cell
2261  {
2262  point[0] = cell->center()[0];
2263  point[1] = cell->center()[1];
2264  point[2] = 0;
2265 
2267  {
2268  point[2] = svg_flags.level_height_factor *
2269  (static_cast<float>(cell->level()) /
2270  static_cast<float>(n_levels)) *
2271  std::max(x_dimension, y_dimension);
2272  }
2273 
2274  const double distance_to_camera =
2275  std::sqrt(std::pow(point[0] - camera_position[0], 2.) +
2276  std::pow(point[1] - camera_position[1], 2.) +
2277  std::pow(point[2] - camera_position[2], 2.));
2278  const double distance_factor =
2279  distance_to_camera / (2. * std::max(x_dimension, y_dimension));
2280 
2281  projection_decomposition = svg_project_point(point,
2282  camera_position,
2283  camera_direction,
2284  camera_horizontal,
2285  camera_focus);
2286 
2287  const unsigned int font_size_this_cell =
2288  static_cast<unsigned int>(
2289  .5 +
2290  cell_label_font_size *
2291  std::pow(.5, cell->level() - 4. + 3.5 * distance_factor));
2292 
2293  out << " <text"
2294  << " x=\""
2295  << static_cast<unsigned int>(
2296  .5 +
2297  ((projection_decomposition[0] - x_min_perspective) /
2298  x_dimension_perspective) *
2299  (width - (width / 100.) * 2. * margin_in_percent) +
2300  ((width / 100.) * margin_in_percent))
2301  << "\" y=\""
2302  << static_cast<unsigned int>(
2303  .5 + height - (height / 100.) * margin_in_percent -
2304  ((projection_decomposition[1] - y_min_perspective) /
2305  y_dimension_perspective) *
2306  (height - (height / 100.) * 2. * margin_in_percent) +
2307  0.5 * font_size_this_cell)
2308  << "\" style=\"font-size:" << font_size_this_cell << "px\">";
2309 
2311  {
2312  out << cell->level();
2313  }
2314 
2316  {
2318  out << '.';
2319  out << cell->index();
2320  }
2321 
2323  {
2326  out << ',';
2327  out
2328  << static_cast<std::make_signed<types::material_id>::type>(
2329  cell->material_id());
2330  }
2331 
2333  {
2336  out << ',';
2337  if (cell->is_active())
2338  out << static_cast<
2339  std::make_signed<types::subdomain_id>::type>(
2340  cell->subdomain_id());
2341  else
2342  out << 'X';
2343  }
2344 
2346  {
2351  out << ',';
2352  out
2353  << static_cast<std::make_signed<types::subdomain_id>::type>(
2354  cell->level_subdomain_id());
2355  }
2356 
2357  out << "</text>" << '\n';
2358  }
2359 
2360  // if the current cell lies at the boundary of the triangulation, draw
2361  // the additional boundary line
2363  {
2364  for (auto faceIndex : cell->face_indices())
2365  {
2366  if (cell->at_boundary(faceIndex))
2367  {
2368  point[0] = cell->face(faceIndex)->vertex(0)[0];
2369  point[1] = cell->face(faceIndex)->vertex(0)[1];
2370  point[2] = 0;
2371 
2373  {
2374  point[2] = svg_flags.level_height_factor *
2375  (static_cast<float>(cell->level()) /
2376  static_cast<float>(n_levels)) *
2377  std::max(x_dimension, y_dimension);
2378  }
2379 
2380  projection_decomposition =
2381  svg_project_point(point,
2382  camera_position,
2383  camera_direction,
2384  camera_horizontal,
2385  camera_focus);
2386 
2387  out << " <line x1=\""
2388  << static_cast<unsigned int>(
2389  .5 +
2390  ((projection_decomposition[0] -
2391  x_min_perspective) /
2392  x_dimension_perspective) *
2393  (width -
2394  (width / 100.) * 2. * margin_in_percent) +
2395  ((width / 100.) * margin_in_percent))
2396  << "\" y1=\""
2397  << static_cast<unsigned int>(
2398  .5 + height -
2399  (height / 100.) * margin_in_percent -
2400  ((projection_decomposition[1] -
2401  y_min_perspective) /
2402  y_dimension_perspective) *
2403  (height -
2404  (height / 100.) * 2. * margin_in_percent));
2405 
2406  point[0] = cell->face(faceIndex)->vertex(1)[0];
2407  point[1] = cell->face(faceIndex)->vertex(1)[1];
2408  point[2] = 0;
2409 
2411  {
2412  point[2] = svg_flags.level_height_factor *
2413  (static_cast<float>(cell->level()) /
2414  static_cast<float>(n_levels)) *
2415  std::max(x_dimension, y_dimension);
2416  }
2417 
2418  projection_decomposition =
2419  svg_project_point(point,
2420  camera_position,
2421  camera_direction,
2422  camera_horizontal,
2423  camera_focus);
2424 
2425  out << "\" x2=\""
2426  << static_cast<unsigned int>(
2427  .5 +
2428  ((projection_decomposition[0] -
2429  x_min_perspective) /
2430  x_dimension_perspective) *
2431  (width -
2432  (width / 100.) * 2. * margin_in_percent) +
2433  ((width / 100.) * margin_in_percent))
2434  << "\" y2=\""
2435  << static_cast<unsigned int>(
2436  .5 + height -
2437  (height / 100.) * margin_in_percent -
2438  ((projection_decomposition[1] -
2439  y_min_perspective) /
2440  y_dimension_perspective) *
2441  (height -
2442  (height / 100.) * 2. * margin_in_percent))
2443  << "\"/>" << '\n';
2444 
2445 
2447  {
2448  const double distance_to_camera = std::sqrt(
2449  std::pow(point[0] - camera_position[0], 2.) +
2450  std::pow(point[1] - camera_position[1], 2.) +
2451  std::pow(point[2] - camera_position[2], 2.));
2452  const double distance_factor =
2453  distance_to_camera /
2454  (2. * std::max(x_dimension, y_dimension));
2455 
2456  const unsigned int font_size_this_edge =
2457  static_cast<unsigned int>(
2458  .5 + .5 * cell_label_font_size *
2459  std::pow(.5,
2460  cell->level() - 4. +
2461  3.5 * distance_factor));
2462 
2463  point[0] = cell->face(faceIndex)->center()[0];
2464  point[1] = cell->face(faceIndex)->center()[1];
2465  point[2] = 0;
2466 
2468  {
2469  point[2] = svg_flags.level_height_factor *
2470  (static_cast<float>(cell->level()) /
2471  static_cast<float>(n_levels)) *
2472  std::max(x_dimension, y_dimension);
2473  }
2474 
2475  projection_decomposition =
2476  svg_project_point(point,
2477  camera_position,
2478  camera_direction,
2479  camera_horizontal,
2480  camera_focus);
2481 
2482  const unsigned int xc = static_cast<unsigned int>(
2483  .5 +
2484  ((projection_decomposition[0] - x_min_perspective) /
2485  x_dimension_perspective) *
2486  (width -
2487  (width / 100.) * 2. * margin_in_percent) +
2488  ((width / 100.) * margin_in_percent));
2489  const unsigned int yc = static_cast<unsigned int>(
2490  .5 + height - (height / 100.) * margin_in_percent -
2491  ((projection_decomposition[1] - y_min_perspective) /
2492  y_dimension_perspective) *
2493  (height -
2494  (height / 100.) * 2. * margin_in_percent));
2495 
2496  out << " <circle cx=\"" << xc << "\" cy=\"" << yc
2497  << "\" r=\"" << font_size_this_edge << "\" />"
2498  << '\n';
2499 
2500  out << " <text x=\"" << xc << "\" y=\"" << yc
2501  << "\" style=\"font-size:" << font_size_this_edge
2502  << "px\" dominant-baseline=\"middle\">"
2503  << static_cast<int>(
2504  cell->face(faceIndex)->boundary_id())
2505  << "</text>" << '\n';
2506  }
2507  }
2508  }
2509  }
2510  }
2511  }
2512 
2513 
2514 
2515  // draw the legend
2516  if (svg_flags.draw_legend)
2517  out << '\n' << " <!-- legend -->" << '\n';
2518 
2519  additional_width = 0;
2520  if (!svg_flags.margin)
2521  additional_width = static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2522 
2523  // explanation of the cell labeling
2524  if (svg_flags.draw_legend &&
2528  {
2529  unsigned int line_offset = 0;
2530  out << " <rect x=\"" << width + additional_width << "\" y=\""
2531  << static_cast<unsigned int>(.5 + (height / 100.) * margin_in_percent)
2532  << "\" width=\""
2533  << static_cast<unsigned int>(.5 + (height / 100.) *
2534  (40. - margin_in_percent))
2535  << "\" height=\"" << static_cast<unsigned int>(.5 + height * .215)
2536  << "\"/>" << '\n';
2537 
2538  out << " <text x=\""
2539  << width + additional_width +
2540  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2541  << "\" y=\""
2542  << static_cast<unsigned int>(.5 +
2543  (height / 100.) * margin_in_percent +
2544  (++line_offset) * 1.5 * font_size)
2545  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2546  << font_size << "px\">"
2547  << "cell label"
2548  << "</text>" << '\n';
2549 
2551  {
2552  out << " <text x=\""
2553  << width + additional_width +
2554  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2555  << "\" y=\""
2556  << static_cast<unsigned int>(.5 +
2557  (height / 100.) * margin_in_percent +
2558  (++line_offset) * 1.5 * font_size)
2559  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2560  << font_size << "px\">"
2561  << "cell_level";
2562 
2566  out << '.';
2567 
2568  out << "</text>" << '\n';
2569  }
2570 
2572  {
2573  out << " <text x=\""
2574  << width + additional_width +
2575  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2576  << "\" y=\""
2577  << static_cast<unsigned int>(.5 +
2578  (height / 100.) * margin_in_percent +
2579  (++line_offset) * 1.5 * font_size)
2580  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2581  << font_size << "px\">"
2582  << "cell_index";
2583 
2586  out << ',';
2587 
2588  out << "</text>" << '\n';
2589  }
2590 
2592  {
2593  out << " <text x=\""
2594  << width + additional_width +
2595  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2596  << "\" y=\""
2597  << static_cast<unsigned int>(.5 +
2598  (height / 100.) * margin_in_percent +
2599  (++line_offset) * 1.5 * font_size)
2600  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2601  << font_size << "px\">"
2602  << "material_id";
2603 
2606  out << ',';
2607 
2608  out << "</text>" << '\n';
2609  }
2610 
2612  {
2613  out << " <text x= \""
2614  << width + additional_width +
2615  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2616  << "\" y=\""
2617  << static_cast<unsigned int>(.5 +
2618  (height / 100.) * margin_in_percent +
2619  (++line_offset) * 1.5 * font_size)
2620  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2621  << font_size << "px\">"
2622  << "subdomain_id";
2623 
2625  out << ',';
2626 
2627  out << "</text>" << '\n';
2628  }
2629 
2631  {
2632  out << " <text x= \""
2633  << width + additional_width +
2634  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2635  << "\" y=\""
2636  << static_cast<unsigned int>(.5 +
2637  (height / 100.) * margin_in_percent +
2638  (++line_offset) * 1.5 * font_size)
2639  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2640  << font_size << "px\">"
2641  << "level_subdomain_id"
2642  << "</text>" << '\n';
2643  }
2644 
2646  {
2647  out << " <text x=\""
2648  << width + additional_width +
2649  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2650  << "\" y=\""
2651  << static_cast<unsigned int>(.5 +
2652  (height / 100.) * margin_in_percent +
2653  (++line_offset) * 1.5 * font_size)
2654  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2655  << font_size << "px\">"
2656  << "edge label"
2657  << "</text>" << '\n';
2658 
2659  out << " <text x= \""
2660  << width + additional_width +
2661  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2662  << "\" y=\""
2663  << static_cast<unsigned int>(.5 +
2664  (height / 100.) * margin_in_percent +
2665  (++line_offset) * 1.5 * font_size)
2666  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2667  << font_size << "px\">"
2668  << "boundary_id"
2669  << "</text>" << '\n';
2670  }
2671  }
2672 
2673  // show azimuth angle and polar angle as text below the explanation of the
2674  // cell labeling
2675  if (svg_flags.draw_legend)
2676  {
2677  out << " <text x=\"" << width + additional_width << "\" y=\""
2678  << static_cast<unsigned int>(
2679  .5 + (height / 100.) * margin_in_percent + 13.75 * font_size)
2680  << "\" style=\"text-anchor:start; font-size:" << font_size << "px\">"
2681  << "azimuth: " << svg_flags.azimuth_angle
2682  << "°, polar: " << svg_flags.polar_angle << "°</text>" << '\n';
2683  }
2684 
2685 
2686  // draw the colorbar
2688  {
2689  out << '\n' << " <!-- colorbar -->" << '\n';
2690 
2691  out << " <text x=\"" << width + additional_width << "\" y=\""
2692  << static_cast<unsigned int>(
2693  .5 + (height / 100.) * (margin_in_percent + 29.) -
2694  (font_size / 1.25))
2695  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2696  << font_size << "px\">";
2697 
2698  switch (svg_flags.coloring)
2699  {
2700  case 1:
2701  out << "material_id";
2702  break;
2703  case 2:
2704  out << "level_number";
2705  break;
2706  case 3:
2707  out << "subdomain_id";
2708  break;
2709  case 4:
2710  out << "level_subdomain_id";
2711  break;
2712  default:
2713  break;
2714  }
2715 
2716  out << "</text>" << '\n';
2717 
2718  unsigned int element_height = static_cast<unsigned int>(
2719  ((height / 100.) * (71. - 2. * margin_in_percent)) / n);
2720  unsigned int element_width =
2721  static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2722 
2723  int labeling_index = 0;
2724  auto materials_it = materials.begin();
2725  auto levels_it = levels.begin();
2726  auto subdomains_it = subdomains.begin();
2727  auto level_subdomains_it = level_subdomains.begin();
2728 
2729  for (unsigned int index = 0; index < n; index++)
2730  {
2731  switch (svg_flags.coloring)
2732  {
2734  labeling_index = *materials_it++;
2735  break;
2737  labeling_index = *levels_it++;
2738  break;
2740  labeling_index = *subdomains_it++;
2741  break;
2743  labeling_index = *level_subdomains_it++;
2744  break;
2745  default:
2746  break;
2747  }
2748 
2749  out << " <rect class=\"r" << labeling_index << "\" x=\""
2750  << width + additional_width << "\" y=\""
2751  << static_cast<unsigned int>(.5 + (height / 100.) *
2752  (margin_in_percent + 29)) +
2753  (n - index - 1) * element_height
2754  << "\" width=\"" << element_width << "\" height=\""
2755  << element_height << "\"/>" << '\n';
2756 
2757  out << " <text x=\""
2758  << width + additional_width + 1.5 * element_width << "\" y=\""
2759  << static_cast<unsigned int>(.5 + (height / 100.) *
2760  (margin_in_percent + 29)) +
2761  (n - index - 1 + .5) * element_height +
2762  static_cast<unsigned int>(.5 + font_size * .35)
2763  << "\""
2764  << " style=\"text-anchor:start; font-size:"
2765  << static_cast<unsigned int>(.5 + font_size) << "px";
2766 
2767  if (index == 0 || index == n - 1)
2768  out << "; font-weight:bold";
2769 
2770  out << "\">" << labeling_index;
2771 
2772  if (index == n - 1)
2773  out << " max";
2774  if (index == 0)
2775  out << " min";
2776 
2777  out << "</text>" << '\n';
2778 
2779  labeling_index++;
2780  }
2781  }
2782 
2783 
2784  // finalize the svg file
2785  out << '\n' << "</svg>";
2786  out.flush();
2787 }
2788 
2789 
2790 template <>
2791 void
2792 GridOut::write_mathgl(const Triangulation<1> &, std::ostream &) const
2793 {
2794  // 1d specialization not done yet
2795  Assert(false, ExcNotImplemented());
2796 }
2797 
2798 
2799 template <int dim, int spacedim>
2800 void
2802  std::ostream & out) const
2803 {
2804  AssertThrow(out, ExcIO());
2805 
2806  // (i) write header
2807  {
2808  // block this to have local variables destroyed after use
2809  const std::time_t time1 = std::time(nullptr);
2810  const std::tm * time = std::localtime(&time1);
2811 
2812  out
2813  << "\n#"
2814  << "\n# This file was generated by the deal.II library."
2815  << "\n# Date = " << time->tm_year + 1900 << "/" << std::setfill('0')
2816  << std::setw(2) << time->tm_mon + 1 << "/" << std::setfill('0')
2817  << std::setw(2) << time->tm_mday << "\n# Time = " << std::setfill('0')
2818  << std::setw(2) << time->tm_hour << ":" << std::setfill('0')
2819  << std::setw(2) << time->tm_min << ":" << std::setfill('0')
2820  << std::setw(2) << time->tm_sec << "\n#"
2821  << "\n# For a description of the MathGL script format see the MathGL manual. "
2822  << "\n#"
2823  << "\n# Note: This file is understood by MathGL v2.1 and higher only, and can "
2824  << "\n# be quickly viewed in a graphical environment using \'mglview\'. "
2825  << "\n#"
2826  << "\n";
2827  }
2828 
2829  // define a helper to keep loops approximately dim-independent
2830  // since MathGL labels axes as x, y, z
2831  const std::string axes = "xyz";
2832 
2833  // (ii) write preamble and graphing tweaks
2834  out << "\n#"
2835  << "\n# Preamble."
2836  << "\n#"
2837  << "\n";
2838 
2840  out << "\nbox";
2841 
2842  // deal with dimension dependent preamble; eg. default sizes and
2843  // views for MathGL (cf. gnuplot).
2844  switch (dim)
2845  {
2846  case 2:
2847  out << "\nsetsize 800 800";
2848  out << "\nrotate 0 0";
2849  break;
2850  case 3:
2851  out << "\nsetsize 800 800";
2852  out << "\nrotate 60 40";
2853  break;
2854  default:
2855  Assert(false, ExcNotImplemented());
2856  }
2857  out << "\n";
2858 
2859  // (iii) write vertex ordering
2860  out << "\n#"
2861  << "\n# Vertex ordering."
2862  << "\n# list <vertex order> <vertex indices>"
2863  << "\n#"
2864  << "\n";
2865 
2866  // todo: This denotes the natural ordering of vertices, but it needs
2867  // to check this is really always true for a given grid (it's not
2868  // true in @ref step_1 "step-1" grid-2 for instance).
2869  switch (dim)
2870  {
2871  case 2:
2872  out << "\nlist f 0 1 2 3"
2873  << "\n";
2874  break;
2875  case 3:
2876  out
2877  << "\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"
2878  << "\n";
2879  break;
2880  default:
2881  Assert(false, ExcNotImplemented());
2882  }
2883 
2884  // (iv) write a list of vertices of cells
2885  out << "\n#"
2886  << "\n# List of vertices."
2887  << "\n# list <id> <vertices>"
2888  << "\n#"
2889  << "\n";
2890 
2891  // run over all active cells and write out a list of
2892  // xyz-coordinates that correspond to vertices
2893  // No global indices in deal.II, so we make one up here.
2894  for (const auto &cell : tria.active_cell_iterators())
2895  {
2896  for (unsigned int i = 0; i < dim; ++i)
2897  {
2898  // if (cell->direction_flag ()==true)
2899  // out << "\ntrue";
2900  // else
2901  // out << "\nfalse";
2902 
2903  out << "\nlist " << axes[i] << cell->active_cell_index() << " ";
2904  for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
2905  out << cell->vertex(j)[i] << " ";
2906  }
2907  out << '\n';
2908  }
2909 
2910  // (v) write out cells to plot as quadplot objects
2911  out << "\n#"
2912  << "\n# List of cells to quadplot."
2913  << "\n# quadplot <vertex order> <id> <style>"
2914  << "\n#"
2915  << "\n";
2916  for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
2917  {
2918  out << "\nquadplot f ";
2919  for (unsigned int j = 0; j < dim; ++j)
2920  out << axes[j] << i << " ";
2921  out << "\'k#\'";
2922  }
2923  out << "\n";
2924 
2925  // (vi) write footer
2926  out << "\n#"
2927  << "\n#"
2928  << "\n#"
2929  << "\n";
2930 
2931  // make sure everything now gets to the output stream
2932  out.flush();
2933  AssertThrow(out, ExcIO());
2934 }
2935 
2936 
2937 
2938 namespace
2939 {
2946  template <int dim, int spacedim, typename ITERATOR, typename END>
2947  void
2948  generate_triangulation_patches(
2949  std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
2950  ITERATOR cell,
2951  END end)
2952  {
2953  // convert each of the active cells into a patch
2954  for (; cell != end; ++cell)
2955  {
2957  patch.reference_cell_type = cell->reference_cell_type();
2958  patch.n_subdivisions = 1;
2959  patch.data.reinit(5, cell->n_vertices());
2960 
2961  for (const unsigned int v : cell->vertex_indices())
2962  {
2963  patch.vertices[v] = cell->vertex(v);
2964  patch.data(0, v) = cell->level();
2965  patch.data(1, v) =
2966  static_cast<std::make_signed<types::manifold_id>::type>(
2967  cell->manifold_id());
2968  patch.data(2, v) =
2969  static_cast<std::make_signed<types::material_id>::type>(
2970  cell->material_id());
2971  if (cell->is_active())
2972  patch.data(3, v) =
2973  static_cast<std::make_signed<types::subdomain_id>::type>(
2974  cell->subdomain_id());
2975  else
2976  patch.data(3, v) = -1;
2977  patch.data(4, v) =
2978  static_cast<std::make_signed<types::subdomain_id>::type>(
2979  cell->level_subdomain_id());
2980  }
2981  patches.push_back(patch);
2982  }
2983  }
2984 
2985 
2986 
2987  std::vector<std::string>
2988  triangulation_patch_data_names()
2989  {
2990  std::vector<std::string> v(5);
2991  v[0] = "level";
2992  v[1] = "manifold";
2993  v[2] = "material";
2994  v[3] = "subdomain";
2995  v[4] = "level_subdomain";
2996  return v;
2997  }
2998 
3002  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3003  get_boundary_edge_iterators(const Triangulation<3, 3> &tria)
3004  {
3005  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3006 
3007  std::vector<bool> flags;
3008  tria.save_user_flags_line(flags);
3009  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3010 
3011  for (auto face : tria.active_face_iterators())
3012  for (const auto l : face->line_indices())
3013  {
3014  const auto line = face->line(l);
3015  if (line->user_flag_set() || line->has_children())
3016  continue;
3017  else
3018  line->set_user_flag();
3019  if (line->at_boundary())
3020  res.emplace_back(line);
3021  }
3022  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3023  return res;
3024  }
3025 
3026 
3027 
3031  template <int dim, int spacedim>
3032  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3033  get_boundary_edge_iterators(const Triangulation<dim, spacedim> &)
3034  {
3035  return {};
3036  }
3037 
3038 
3039 
3044  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3045  get_relevant_edge_iterators(const Triangulation<3, 3> &tria)
3046  {
3047  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3048 
3049  std::vector<bool> flags;
3050  tria.save_user_flags_line(flags);
3051  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3052 
3053  for (auto face : tria.active_face_iterators())
3054  for (const auto l : face->line_indices())
3055  {
3056  const auto line = face->line(l);
3057  if (line->user_flag_set() || line->has_children())
3058  continue;
3059  else
3060  line->set_user_flag();
3061  if (line->manifold_id() != numbers::flat_manifold_id ||
3062  (line->boundary_id() != 0 &&
3063  line->boundary_id() != numbers::invalid_boundary_id))
3064  res.emplace_back(line);
3065  }
3066  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3067  return res;
3068  }
3069 
3070 
3074  template <int dim, int spacedim>
3075  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3076  get_relevant_edge_iterators(const Triangulation<dim, spacedim> &)
3077  {
3078  return {};
3079  }
3080 
3081 
3082 
3086  template <int dim, int spacedim>
3087  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3088  get_boundary_face_iterators(const Triangulation<dim, spacedim> &tria)
3089  {
3090  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3091  res;
3092  if (dim == 1)
3093  return res;
3094  for (auto face : tria.active_face_iterators())
3095  {
3096  if (face->boundary_id() != numbers::invalid_boundary_id)
3097  res.push_back(face);
3098  }
3099  return res;
3100  }
3101 
3102 
3103 
3108  template <int dim, int spacedim>
3109  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3110  get_relevant_face_iterators(const Triangulation<dim, spacedim> &tria)
3111  {
3112  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3113  res;
3114  if (dim == 1)
3115  return res;
3116  for (auto face : tria.active_face_iterators())
3117  {
3118  if (face->manifold_id() != numbers::flat_manifold_id ||
3119  (face->boundary_id() != 0 &&
3120  face->boundary_id() != numbers::invalid_boundary_id))
3121  res.push_back(face);
3122  }
3123  return res;
3124  }
3125 } // namespace
3126 
3127 
3128 
3129 template <int dim, int spacedim>
3130 void
3132  std::ostream & out) const
3133 {
3134  AssertThrow(out, ExcIO());
3135 
3136  // get the positions of the vertices
3137  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
3138 
3139  const auto n_vertices = vertices.size();
3140 
3141  out << "# vtk DataFile Version 3.0\n"
3142  << "Triangulation generated with deal.II\n"
3143  << "ASCII\n"
3144  << "DATASET UNSTRUCTURED_GRID\n"
3145  << "POINTS " << n_vertices << " double\n";
3146 
3147  // actually write the vertices.
3148  for (const auto &v : vertices)
3149  {
3150  out << v;
3151  for (unsigned int d = spacedim + 1; d <= 3; ++d)
3152  out << " 0"; // fill with zeroes
3153  out << '\n';
3154  }
3155 
3156  const auto faces = vtk_flags.output_only_relevant ?
3157  get_relevant_face_iterators(tria) :
3158  get_boundary_face_iterators(tria);
3159  const auto edges = vtk_flags.output_only_relevant ?
3160  get_relevant_edge_iterators(tria) :
3161  get_boundary_edge_iterators(tria);
3162 
3163  AssertThrow(
3164  vtk_flags.output_cells || (dim >= 2 && vtk_flags.output_faces) ||
3165  (dim >= 3 && vtk_flags.output_edges),
3166  ExcMessage(
3167  "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3168 
3169  // Write cells preamble
3170  const int n_cells = (vtk_flags.output_cells ? tria.n_active_cells() : 0) +
3171  (vtk_flags.output_faces ? faces.size() : 0) +
3172  (vtk_flags.output_edges ? edges.size() : 0);
3173 
3174  // VTK now expects a number telling the total storage requirement to read all
3175  // cell connectivity information. The connectivity information is read cell by
3176  // cell, first specifying how many vertices are required to describe the cell,
3177  // and then specifying the index of every vertex. This means that for every
3178  // deal.II object type, we always need n_vertices + 1 integer per cell.
3179  // Compute the total number here.
3180  int cells_size = 0;
3181 
3182  if (vtk_flags.output_cells)
3183  for (const auto &cell : tria.active_cell_iterators())
3184  cells_size += cell->n_vertices() + 1;
3185 
3186  if (vtk_flags.output_faces)
3187  for (const auto &face : faces)
3188  cells_size += face->n_vertices() + 1;
3189 
3190  if (vtk_flags.output_edges)
3191  for (const auto &edge : edges)
3192  cells_size += edge->n_vertices() + 1;
3193 
3194  AssertThrow(cells_size > 0, ExcMessage("No cells given to be output!"));
3195 
3196  out << "\nCELLS " << n_cells << ' ' << cells_size << '\n';
3197  /*
3198  * VTK cells:
3199  *
3200  * 1 VTK_VERTEX
3201  * 3 VTK_LINE
3202  * 5 VTK_TRIANGLE
3203  * 9 VTK_QUAD
3204  * 10 VTK_TETRA
3205  * 14 VTK_PYRAMID
3206  * 13 VTK_WEDGE
3207  * 12 VTK_HEXAHEDRON
3208  *
3209  * see also: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
3210  */
3211  static const std::array<int, 8> table = {{1, 3, 5, 9, 10, 14, 13, 12}};
3212 
3213  // write cells.
3214  if (vtk_flags.output_cells)
3215  for (const auto &cell : tria.active_cell_iterators())
3216  {
3217  out << cell->n_vertices();
3218  for (const unsigned int i : cell->vertex_indices())
3219  {
3220  out << ' '
3221  << cell->vertex_index(GeometryInfo<dim>::vertices_per_cell ==
3222  cell->n_vertices() ?
3224  i);
3225  }
3226  out << '\n';
3227  }
3228  if (vtk_flags.output_faces)
3229  for (const auto &face : faces)
3230  {
3231  out << face->n_vertices();
3232  constexpr int face_dim = dim > 1 ? dim - 1 : 1;
3233  for (const unsigned int i : face->vertex_indices())
3234  {
3235  out << ' '
3236  << face->vertex_index(GeometryInfo<dim>::vertices_per_face ==
3237  face->n_vertices() ?
3239  i);
3240  }
3241  out << '\n';
3242  }
3243  if (vtk_flags.output_edges)
3244  for (const auto &edge : edges)
3245  {
3246  out << 2;
3247  for (const unsigned int i : edge->vertex_indices())
3248  out << ' ' << edge->vertex_index(i);
3249  out << '\n';
3250  }
3251 
3252  // write cell types
3253  out << "\nCELL_TYPES " << n_cells << '\n';
3254  if (vtk_flags.output_cells)
3255  {
3256  for (const auto &cell : tria.active_cell_iterators())
3257  out << table[static_cast<int>(cell->reference_cell_type())] << ' ';
3258  out << '\n';
3259  }
3260  if (vtk_flags.output_faces)
3261  {
3262  for (const auto &face : faces)
3263  out << table[static_cast<int>(face->reference_cell_type())] << ' ';
3264  out << '\n';
3265  }
3266  if (vtk_flags.output_edges)
3267  {
3268  for (const auto &edge : edges)
3269  out << table[static_cast<int>(edge->reference_cell_type())] << ' ';
3270  }
3271  out << "\n\nCELL_DATA " << n_cells << '\n'
3272  << "SCALARS MaterialID int 1\n"
3273  << "LOOKUP_TABLE default\n";
3274 
3275  // Now material id and boundary id
3276  if (vtk_flags.output_cells)
3277  {
3278  for (const auto &cell : tria.active_cell_iterators())
3279  {
3280  out << static_cast<std::make_signed<types::material_id>::type>(
3281  cell->material_id())
3282  << ' ';
3283  }
3284  out << '\n';
3285  }
3286  if (vtk_flags.output_faces)
3287  {
3288  for (const auto &face : faces)
3289  {
3290  out << static_cast<std::make_signed<types::boundary_id>::type>(
3291  face->boundary_id())
3292  << ' ';
3293  }
3294  out << '\n';
3295  }
3296  if (vtk_flags.output_edges)
3297  {
3298  for (const auto &edge : edges)
3299  {
3300  out << static_cast<std::make_signed<types::boundary_id>::type>(
3301  edge->boundary_id())
3302  << ' ';
3303  }
3304  }
3305 
3306  out << "\n\nSCALARS ManifoldID int 1\n"
3307  << "LOOKUP_TABLE default\n";
3308 
3309  // Now manifold id
3310  if (vtk_flags.output_cells)
3311  {
3312  for (const auto &cell : tria.active_cell_iterators())
3313  {
3314  out << static_cast<std::make_signed<types::manifold_id>::type>(
3315  cell->manifold_id())
3316  << ' ';
3317  }
3318  out << '\n';
3319  }
3320  if (vtk_flags.output_faces)
3321  {
3322  for (const auto &face : faces)
3323  {
3324  out << static_cast<std::make_signed<types::manifold_id>::type>(
3325  face->manifold_id())
3326  << ' ';
3327  }
3328  out << '\n';
3329  }
3330  if (vtk_flags.output_edges)
3331  {
3332  for (const auto &edge : edges)
3333  {
3334  out << static_cast<std::make_signed<types::manifold_id>::type>(
3335  edge->manifold_id())
3336  << ' ';
3337  }
3338  out << '\n';
3339  }
3340 
3341  out.flush();
3342 
3343  AssertThrow(out, ExcIO());
3344 }
3345 
3346 
3347 
3348 template <int dim, int spacedim>
3349 void
3351  std::ostream & out) const
3352 {
3353  AssertThrow(out, ExcIO());
3354 
3355  // convert the cells of the triangulation into a set of patches
3356  // and then have them output. since there is no data attached to
3357  // the geometry, we also do not have to provide any names, identifying
3358  // information, etc.
3359  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3360  patches.reserve(tria.n_active_cells());
3361  generate_triangulation_patches(patches, tria.begin_active(), tria.end());
3362 
3365  patches,
3366  triangulation_patch_data_names(),
3367  std::vector<
3368  std::tuple<unsigned int,
3369  unsigned int,
3370  std::string,
3372  vtu_flags,
3373  out);
3375  {
3376  out << " </UnstructuredGrid>\n";
3377  out << "<dealiiData encoding=\"base64\">";
3378  std::stringstream outstring;
3379  boost::archive::binary_oarchive ia(outstring);
3380  tria.save(ia, 0);
3381  const auto compressed = Utilities::compress(outstring.str());
3382  out << Utilities::encode_base64({compressed.begin(), compressed.end()});
3383  out << "\n</dealiiData>\n";
3384  out << "</VTKFile>\n";
3385  }
3386  else
3388 
3389  out << std::flush;
3390  AssertThrow(out, ExcIO());
3391 }
3392 
3393 
3394 
3395 template <int dim, int spacedim>
3396 void
3398  const Triangulation<dim, spacedim> &tria,
3399  const std::string & filename_without_extension,
3400  const bool view_levels,
3401  const bool include_artificial) const
3402 {
3403  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3404  const unsigned int n_datasets = 4;
3405  std::vector<std::string> data_names;
3406  data_names.emplace_back("level");
3407  data_names.emplace_back("subdomain");
3408  data_names.emplace_back("level_subdomain");
3409  data_names.emplace_back("proc_writing");
3410 
3411  const unsigned int n_q_points = GeometryInfo<dim>::vertices_per_cell;
3412 
3413  for (const auto &cell : tria.cell_iterators())
3414  {
3415  if (!view_levels)
3416  {
3417  if (cell->has_children())
3418  continue;
3419  if (!include_artificial &&
3420  cell->subdomain_id() == numbers::artificial_subdomain_id)
3421  continue;
3422  }
3423  else if (!include_artificial)
3424  {
3425  if (cell->has_children() &&
3426  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
3427  continue;
3428  else if (cell->is_active() &&
3429  cell->level_subdomain_id() ==
3431  cell->subdomain_id() == numbers::artificial_subdomain_id)
3432  continue;
3433  }
3434 
3436  patch.data.reinit(n_datasets, n_q_points);
3437  patch.points_are_available = false;
3438 
3439  for (unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3440  {
3441  patch.vertices[vertex] = cell->vertex(vertex);
3442  patch.data(0, vertex) = cell->level();
3443  if (cell->is_active())
3444  patch.data(1, vertex) = static_cast<double>(
3445  static_cast<std::make_signed<types::subdomain_id>::type>(
3446  cell->subdomain_id()));
3447  else
3448  patch.data(1, vertex) = -1.0;
3449  patch.data(2, vertex) = static_cast<double>(
3450  static_cast<std::make_signed<types::subdomain_id>::type>(
3451  cell->level_subdomain_id()));
3452  patch.data(3, vertex) = tria.locally_owned_subdomain();
3453  }
3454 
3455  for (auto f : GeometryInfo<dim>::face_indices())
3457  patches.push_back(patch);
3458  }
3459 
3460  // only create .pvtu file if running in parallel
3461  // if not, just create a .vtu file with no reference
3462  // to the processor number
3463  std::string new_file = filename_without_extension + ".vtu";
3465  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
3466  {
3467  new_file = filename_without_extension + ".proc" +
3468  Utilities::int_to_string(tr->locally_owned_subdomain(), 4) +
3469  ".vtu";
3470 
3471  // create .pvtu record
3472  if (tr->locally_owned_subdomain() == 0)
3473  {
3474  std::vector<std::string> filenames;
3475 
3476  // .pvtu needs to reference the files without a relative path because
3477  // it will be written in the same directory. For this, remove any
3478  // paths from filename.
3479  std::size_t pos = filename_without_extension.find_last_of('/');
3480  if (pos == std::string::npos)
3481  pos = 0;
3482  else
3483  pos += 1;
3484  const unsigned int n_procs =
3485  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3486  for (unsigned int i = 0; i < n_procs; ++i)
3487  filenames.push_back(filename_without_extension.substr(pos) +
3488  ".proc" + Utilities::int_to_string(i, 4) +
3489  ".vtu");
3490 
3491  const std::string pvtu_filename =
3492  (filename_without_extension + ".pvtu");
3493  std::ofstream pvtu_output(pvtu_filename.c_str());
3494 
3496  data_out.attach_triangulation(*tr);
3497 
3498  // We need a dummy vector with the names of the data values in the
3499  // .vtu files in order that the .pvtu contains reference these values
3500  Vector<float> dummy_vector(tr->n_active_cells());
3501  data_out.add_data_vector(dummy_vector, "level");
3502  data_out.add_data_vector(dummy_vector, "subdomain");
3503  data_out.add_data_vector(dummy_vector, "level_subdomain");
3504  data_out.add_data_vector(dummy_vector, "proc_writing");
3505 
3506  data_out.build_patches();
3507 
3508  data_out.write_pvtu_record(pvtu_output, filenames);
3509  }
3510  }
3511 
3512  std::ofstream out(new_file.c_str());
3513  std::vector<
3514  std::tuple<unsigned int,
3515  unsigned int,
3516  std::string,
3518  vector_data_ranges;
3519  DataOutBase::VtkFlags flags;
3520  DataOutBase::write_vtu(patches, data_names, vector_data_ranges, flags, out);
3521 }
3522 
3523 
3524 
3525 unsigned int
3527 {
3528  return 0;
3529 }
3530 
3531 unsigned int
3533 {
3534  return 0;
3535 }
3536 
3537 
3538 unsigned int
3540 {
3541  return 0;
3542 }
3543 
3544 unsigned int
3546 {
3547  return 0;
3548 }
3549 
3550 unsigned int
3552 {
3553  return 0;
3554 }
3555 
3556 unsigned int
3558 {
3559  return 0;
3560 }
3561 
3562 unsigned int
3564 {
3565  return 0;
3566 }
3567 
3568 unsigned int
3570 {
3571  return 0;
3572 }
3573 
3574 
3575 
3576 template <int dim, int spacedim>
3577 unsigned int
3579 {
3581  unsigned int n_faces = 0;
3582 
3583  for (const auto &face : tria.active_face_iterators())
3584  if ((face->at_boundary()) && (face->boundary_id() != 0))
3585  n_faces++;
3586 
3587  return n_faces;
3588 }
3589 
3590 
3591 
3592 template <int dim, int spacedim>
3593 unsigned int
3595 {
3596  // save the user flags for lines so
3597  // we can use these flags to track
3598  // which ones we've already counted
3599  std::vector<bool> line_flags;
3600  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3601  line_flags);
3602  const_cast<::Triangulation<dim, spacedim> &>(tria)
3603  .clear_user_flags_line();
3604 
3605  unsigned int n_lines = 0;
3606 
3607  for (const auto &cell : tria.active_cell_iterators())
3608  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3609  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3610  (cell->line(l)->user_flag_set() == false))
3611  {
3612  ++n_lines;
3613  cell->line(l)->set_user_flag();
3614  }
3615 
3616  // at the end, restore the user
3617  // flags for the lines
3618  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3619  line_flags);
3620 
3621  return n_lines;
3622 }
3623 
3624 
3625 
3626 unsigned int
3628  const unsigned int next_element_index,
3629  std::ostream &) const
3630 {
3631  return next_element_index;
3632 }
3633 
3634 
3635 unsigned int
3637  const unsigned int next_element_index,
3638  std::ostream &) const
3639 {
3640  return next_element_index;
3641 }
3642 
3643 unsigned int
3645  const unsigned int next_element_index,
3646  std::ostream &) const
3647 {
3648  return next_element_index;
3649 }
3650 
3651 
3652 unsigned int
3654  const unsigned int next_element_index,
3655  std::ostream &) const
3656 {
3657  return next_element_index;
3658 }
3659 
3660 unsigned int
3662  const unsigned int next_element_index,
3663  std::ostream &) const
3664 {
3665  return next_element_index;
3666 }
3667 
3668 
3669 unsigned int
3671  const unsigned int next_element_index,
3672  std::ostream &) const
3673 {
3674  return next_element_index;
3675 }
3676 
3677 
3678 unsigned int
3680  const unsigned int next_element_index,
3681  std::ostream &) const
3682 {
3683  return next_element_index;
3684 }
3685 
3686 unsigned int
3688  const unsigned int next_element_index,
3689  std::ostream &) const
3690 {
3691  return next_element_index;
3692 }
3693 
3694 
3695 
3696 template <int dim, int spacedim>
3697 unsigned int
3699  const unsigned int next_element_index,
3700  std::ostream & out) const
3701 {
3702  unsigned int current_element_index = next_element_index;
3703 
3704  for (const auto &face : tria.active_face_iterators())
3705  if (face->at_boundary() && (face->boundary_id() != 0))
3706  {
3707  out << current_element_index << ' ';
3708  switch (dim)
3709  {
3710  case 2:
3711  out << 1 << ' ';
3712  break;
3713  case 3:
3714  out << 3 << ' ';
3715  break;
3716  default:
3717  Assert(false, ExcNotImplemented());
3718  }
3719  out << static_cast<unsigned int>(face->boundary_id()) << ' '
3720  << static_cast<unsigned int>(face->boundary_id()) << ' '
3722  // note: vertex numbers are 1-base
3723  for (unsigned int vertex = 0;
3724  vertex < GeometryInfo<dim>::vertices_per_face;
3725  ++vertex)
3726  out << ' '
3727  << face->vertex_index(
3729  1;
3730  out << '\n';
3731 
3732  ++current_element_index;
3733  }
3734  return current_element_index;
3735 }
3736 
3737 
3738 template <int dim, int spacedim>
3739 unsigned int
3741  const unsigned int next_element_index,
3742  std::ostream & out) const
3743 {
3744  unsigned int current_element_index = next_element_index;
3745  // save the user flags for lines so
3746  // we can use these flags to track
3747  // which ones we've already taken
3748  // care of
3749  std::vector<bool> line_flags;
3750  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3751  line_flags);
3752  const_cast<::Triangulation<dim, spacedim> &>(tria)
3753  .clear_user_flags_line();
3754 
3755  for (const auto &cell : tria.active_cell_iterators())
3756  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3757  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3758  (cell->line(l)->user_flag_set() == false))
3759  {
3760  out << next_element_index << " 1 ";
3761  out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
3762  << static_cast<unsigned int>(cell->line(l)->boundary_id())
3763  << " 2 ";
3764  // note: vertex numbers are 1-base
3765  for (unsigned int vertex = 0; vertex < 2; ++vertex)
3766  out << ' '
3767  << cell->line(l)->vertex_index(
3769  1;
3770  out << '\n';
3771 
3772  // move on to the next line
3773  // but mark the current one
3774  // as taken care of
3775  ++current_element_index;
3776  cell->line(l)->set_user_flag();
3777  }
3778 
3779  // at the end, restore the user
3780  // flags for the lines
3781  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3782  line_flags);
3783 
3784  return current_element_index;
3785 }
3786 
3787 
3788 
3789 unsigned int
3791  const unsigned int next_element_index,
3792  std::ostream &) const
3793 {
3794  return next_element_index;
3795 }
3796 
3797 unsigned int
3799  const unsigned int next_element_index,
3800  std::ostream &) const
3801 {
3802  return next_element_index;
3803 }
3804 
3805 unsigned int
3807  const unsigned int next_element_index,
3808  std::ostream &) const
3809 {
3810  return next_element_index;
3811 }
3812 
3813 unsigned int
3815  const unsigned int next_element_index,
3816  std::ostream &) const
3817 {
3818  return next_element_index;
3819 }
3820 
3821 unsigned int
3823  const unsigned int next_element_index,
3824  std::ostream &) const
3825 {
3826  return next_element_index;
3827 }
3828 
3829 
3830 unsigned int
3832  const unsigned int next_element_index,
3833  std::ostream &) const
3834 {
3835  return next_element_index;
3836 }
3837 
3838 
3839 unsigned int
3841  const unsigned int next_element_index,
3842  std::ostream &) const
3843 {
3844  return next_element_index;
3845 }
3846 
3847 unsigned int
3849  const unsigned int next_element_index,
3850  std::ostream &) const
3851 {
3852  return next_element_index;
3853 }
3854 
3855 
3856 
3857 template <int dim, int spacedim>
3858 unsigned int
3860  const unsigned int next_element_index,
3861  std::ostream & out) const
3862 {
3863  unsigned int current_element_index = next_element_index;
3865 
3866  for (const auto &face : tria.active_face_iterators())
3867  if (face->at_boundary() && (face->boundary_id() != 0))
3868  {
3869  out << current_element_index << " "
3870  << static_cast<unsigned int>(face->boundary_id()) << " ";
3871  switch (dim)
3872  {
3873  case 2:
3874  out << "line ";
3875  break;
3876  case 3:
3877  out << "quad ";
3878  break;
3879  default:
3880  Assert(false, ExcNotImplemented());
3881  }
3882  // note: vertex numbers are 1-base
3883  for (unsigned int vertex = 0;
3884  vertex < GeometryInfo<dim>::vertices_per_face;
3885  ++vertex)
3886  out << face->vertex_index(
3888  1
3889  << ' ';
3890  out << '\n';
3891 
3892  ++current_element_index;
3893  }
3894  return current_element_index;
3895 }
3896 
3897 
3898 
3899 template <int dim, int spacedim>
3900 unsigned int
3902  const unsigned int next_element_index,
3903  std::ostream & out) const
3904 {
3905  unsigned int current_element_index = next_element_index;
3906  // save the user flags for lines so
3907  // we can use these flags to track
3908  // which ones we've already taken
3909  // care of
3910  std::vector<bool> line_flags;
3911  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3912  line_flags);
3913  const_cast<::Triangulation<dim, spacedim> &>(tria)
3914  .clear_user_flags_line();
3915 
3916  for (const auto &cell : tria.active_cell_iterators())
3917  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3918  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3919  (cell->line(l)->user_flag_set() == false))
3920  {
3921  out << current_element_index << " "
3922  << static_cast<unsigned int>(cell->line(l)->boundary_id())
3923  << " line ";
3924  // note: vertex numbers in ucd format are 1-base
3925  for (unsigned int vertex = 0; vertex < 2; ++vertex)
3926  out << cell->line(l)->vertex_index(
3928  1
3929  << ' ';
3930  out << '\n';
3931 
3932  // move on to the next line
3933  // but mark the current one
3934  // as taken care of
3935  ++current_element_index;
3936  cell->line(l)->set_user_flag();
3937  }
3938 
3939  // at the end, restore the user
3940  // flags for the lines
3941  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3942  line_flags);
3943  return current_element_index;
3944 }
3945 
3946 
3947 
3948 namespace internal
3949 {
3950  namespace
3951  {
3960  template <int spacedim>
3961  void
3962  remove_colinear_points(std::vector<Point<spacedim>> &points)
3963  {
3964  while (points.size() > 2)
3965  {
3966  Tensor<1, spacedim> first_difference = points[1] - points[0];
3967  first_difference /= first_difference.norm();
3968  Tensor<1, spacedim> second_difference = points[2] - points[1];
3969  second_difference /= second_difference.norm();
3970  // If the three points are colinear then remove the middle one.
3971  if ((first_difference - second_difference).norm() < 1e-10)
3972  points.erase(points.begin() + 1);
3973  else
3974  break;
3975  }
3976  }
3977 
3978 
3979 
3980  template <int spacedim>
3981  void
3982  write_gnuplot(const ::Triangulation<1, spacedim> &tria,
3983  std::ostream & out,
3984  const Mapping<1, spacedim> *,
3986  {
3987  AssertThrow(out, ExcIO());
3988 
3989  for (const auto &cell : tria.active_cell_iterators())
3990  {
3991  if (gnuplot_flags.write_cell_numbers)
3992  out << "# cell " << cell << '\n';
3993 
3994  out << cell->vertex(0) << ' ' << cell->level() << ' '
3995  << cell->material_id() << '\n'
3996  << cell->vertex(1) << ' ' << cell->level() << ' '
3997  << cell->material_id() << '\n'
3998  << "\n\n";
3999  }
4000 
4001  // make sure everything now gets to
4002  // disk
4003  out.flush();
4004 
4005  AssertThrow(out, ExcIO());
4006  }
4007 
4008 
4009 
4010  template <int spacedim>
4011  void
4012  write_gnuplot(const ::Triangulation<2, spacedim> &tria,
4013  std::ostream & out,
4014  const Mapping<2, spacedim> * mapping,
4015  const GridOutFlags::Gnuplot & gnuplot_flags)
4016  {
4017  AssertThrow(out, ExcIO());
4018 
4019  const int dim = 2;
4020 
4021  const unsigned int n_additional_points =
4022  gnuplot_flags.n_extra_curved_line_points;
4023  const unsigned int n_points = 2 + n_additional_points;
4024 
4025  // If we need to plot curved lines then generate a quadrature formula to
4026  // place points via the mapping
4027  Quadrature<dim> q_projector;
4028  std::vector<Point<dim - 1>> boundary_points;
4029  if (mapping != nullptr)
4030  {
4031  boundary_points.resize(n_points);
4032  boundary_points[0][0] = 0;
4033  boundary_points[n_points - 1][0] = 1;
4034  for (unsigned int i = 1; i < n_points - 1; ++i)
4035  boundary_points[i](0) = 1. * i / (n_points - 1);
4036 
4037  std::vector<double> dummy_weights(n_points, 1. / n_points);
4038  Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4039 
4040  q_projector =
4042  quadrature);
4043  }
4044 
4045  for (const auto &cell : tria.active_cell_iterators())
4046  {
4047  if (gnuplot_flags.write_cell_numbers)
4048  out << "# cell " << cell << '\n';
4049 
4050  if (mapping == nullptr ||
4051  (dim == spacedim ?
4052  (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4053  // ignore checking for boundary or interior cells in the codim
4054  // 1 case: 'or false' is a no-op
4055  false))
4056  {
4057  // write out the four sides of this cell by putting the four
4058  // points (+ the initial point again) in a row and lifting the
4059  // drawing pencil at the end
4060  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4061  out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i]) << ' '
4062  << cell->level() << ' ' << cell->material_id() << '\n';
4063  out << cell->vertex(0) << ' ' << cell->level() << ' '
4064  << cell->material_id() << '\n'
4065  << '\n' // double new line for gnuplot 3d plots
4066  << '\n';
4067  }
4068  else
4069  // cell is at boundary and we are to treat curved boundaries. so
4070  // loop over all faces and draw them as small pieces of lines
4071  {
4072  for (const unsigned int face_no :
4074  {
4075  const typename ::Triangulation<dim,
4076  spacedim>::face_iterator
4077  face = cell->face(face_no);
4078  if (dim != spacedim || face->at_boundary() ||
4079  gnuplot_flags.curved_inner_cells)
4080  {
4081  // Save the points on each face to a vector and then try
4082  // to remove colinear points that won't show up in the
4083  // generated plot.
4084  std::vector<Point<spacedim>> line_points;
4085  // compute offset of quadrature points within set of
4086  // projected points
4087  const unsigned int offset = face_no * n_points;
4088  for (unsigned int i = 0; i < n_points; ++i)
4089  line_points.push_back(
4090  mapping->transform_unit_to_real_cell(
4091  cell, q_projector.point(offset + i)));
4092  internal::remove_colinear_points(line_points);
4093 
4094  for (const Point<spacedim> &point : line_points)
4095  out << point << ' ' << cell->level() << ' '
4096  << cell->material_id() << '\n';
4097 
4098  out << '\n' << '\n';
4099  }
4100  else
4101  {
4102  // if, however, the face is not at the boundary and we
4103  // don't want to curve anything, then draw it as usual
4104  out << face->vertex(0) << ' ' << cell->level() << ' '
4105  << cell->material_id() << '\n'
4106  << face->vertex(1) << ' ' << cell->level() << ' '
4107  << cell->material_id() << '\n'
4108  << '\n'
4109  << '\n';
4110  }
4111  }
4112  }
4113  }
4114 
4115  // make sure everything now gets to disk
4116  out.flush();
4117 
4118  AssertThrow(out, ExcIO());
4119  }
4120 
4121 
4122 
4123  template <int spacedim>
4124  void
4125  write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4126  std::ostream & out,
4127  const Mapping<3, spacedim> * mapping,
4128  const GridOutFlags::Gnuplot & gnuplot_flags)
4129  {
4130  AssertThrow(out, ExcIO());
4131 
4132  const int dim = 3;
4133 
4134  const unsigned int n_additional_points =
4135  gnuplot_flags.n_extra_curved_line_points;
4136  const unsigned int n_points = 2 + n_additional_points;
4137 
4138  // If we need to plot curved lines then generate a quadrature formula to
4139  // place points via the mapping
4140  std::unique_ptr<Quadrature<dim>> q_projector;
4141  std::vector<Point<1>> boundary_points;
4142  if (mapping != nullptr)
4143  {
4144  boundary_points.resize(n_points);
4145  boundary_points[0][0] = 0;
4146  boundary_points[n_points - 1][0] = 1;
4147  for (unsigned int i = 1; i < n_points - 1; ++i)
4148  boundary_points[i](0) = 1. * i / (n_points - 1);
4149 
4150  std::vector<double> dummy_weights(n_points, 1. / n_points);
4151  Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4152 
4153  // tensor product of points, only one copy
4154  QIterated<dim - 1> quadrature(quadrature1d, 1);
4155  q_projector = std::make_unique<Quadrature<dim>>(
4157  }
4158 
4159  for (const auto &cell : tria.active_cell_iterators())
4160  {
4161  if (gnuplot_flags.write_cell_numbers)
4162  out << "# cell " << cell << '\n';
4163 
4164  if (mapping == nullptr || n_points == 2 ||
4165  (!cell->has_boundary_lines() &&
4166  !gnuplot_flags.curved_inner_cells))
4167  {
4168  // front face
4169  out << cell->vertex(0) << ' ' << cell->level() << ' '
4170  << cell->material_id() << '\n'
4171  << cell->vertex(1) << ' ' << cell->level() << ' '
4172  << cell->material_id() << '\n'
4173  << cell->vertex(5) << ' ' << cell->level() << ' '
4174  << cell->material_id() << '\n'
4175  << cell->vertex(4) << ' ' << cell->level() << ' '
4176  << cell->material_id() << '\n'
4177  << cell->vertex(0) << ' ' << cell->level() << ' '
4178  << cell->material_id() << '\n'
4179  << '\n';
4180  // back face
4181  out << cell->vertex(2) << ' ' << cell->level() << ' '
4182  << cell->material_id() << '\n'
4183  << cell->vertex(3) << ' ' << cell->level() << ' '
4184  << cell->material_id() << '\n'
4185  << cell->vertex(7) << ' ' << cell->level() << ' '
4186  << cell->material_id() << '\n'
4187  << cell->vertex(6) << ' ' << cell->level() << ' '
4188  << cell->material_id() << '\n'
4189  << cell->vertex(2) << ' ' << cell->level() << ' '
4190  << cell->material_id() << '\n'
4191  << '\n';
4192 
4193  // now for the four connecting lines
4194  out << cell->vertex(0) << ' ' << cell->level() << ' '
4195  << cell->material_id() << '\n'
4196  << cell->vertex(2) << ' ' << cell->level() << ' '
4197  << cell->material_id() << '\n'
4198  << '\n';
4199  out << cell->vertex(1) << ' ' << cell->level() << ' '
4200  << cell->material_id() << '\n'
4201  << cell->vertex(3) << ' ' << cell->level() << ' '
4202  << cell->material_id() << '\n'
4203  << '\n';
4204  out << cell->vertex(5) << ' ' << cell->level() << ' '
4205  << cell->material_id() << '\n'
4206  << cell->vertex(7) << ' ' << cell->level() << ' '
4207  << cell->material_id() << '\n'
4208  << '\n';
4209  out << cell->vertex(4) << ' ' << cell->level() << ' '
4210  << cell->material_id() << '\n'
4211  << cell->vertex(6) << ' ' << cell->level() << ' '
4212  << cell->material_id() << '\n'
4213  << '\n';
4214  }
4215  else
4216  {
4217  for (const unsigned int face_no :
4219  {
4220  const typename ::Triangulation<dim,
4221  spacedim>::face_iterator
4222  face = cell->face(face_no);
4223 
4224  if (face->at_boundary() &&
4225  gnuplot_flags.write_additional_boundary_lines)
4226  {
4227  const unsigned int offset = face_no * n_points * n_points;
4228  for (unsigned int i = 0; i < n_points - 1; ++i)
4229  for (unsigned int j = 0; j < n_points - 1; ++j)
4230  {
4231  const Point<spacedim> p0 =
4232  mapping->transform_unit_to_real_cell(
4233  cell,
4234  q_projector->point(offset + i * n_points + j));
4235  out << p0 << ' ' << cell->level() << ' '
4236  << cell->material_id() << '\n';
4237  out << (mapping->transform_unit_to_real_cell(
4238  cell,
4239  q_projector->point(
4240  offset + (i + 1) * n_points + j)))
4241  << ' ' << cell->level() << ' '
4242  << cell->material_id() << '\n';
4243  out << (mapping->transform_unit_to_real_cell(
4244  cell,
4245  q_projector->point(
4246  offset + (i + 1) * n_points + j + 1)))
4247  << ' ' << cell->level() << ' '
4248  << cell->material_id() << '\n';
4249  out << (mapping->transform_unit_to_real_cell(
4250  cell,
4251  q_projector->point(offset + i * n_points +
4252  j + 1)))
4253  << ' ' << cell->level() << ' '
4254  << cell->material_id() << '\n';
4255  // and the first point again
4256  out << p0 << ' ' << cell->level() << ' '
4257  << cell->material_id() << '\n';
4258  out << '\n' << '\n';
4259  }
4260  }
4261  else
4262  {
4263  for (unsigned int l = 0;
4264  l < GeometryInfo<dim>::lines_per_face;
4265  ++l)
4266  {
4267  const typename ::Triangulation<dim, spacedim>::
4268  line_iterator line = face->line(l);
4269 
4270  const Point<spacedim> &v0 = line->vertex(0),
4271  &v1 = line->vertex(1);
4272  if (line->at_boundary() ||
4273  gnuplot_flags.curved_inner_cells)
4274  {
4275  // Save the points on each face to a vector and
4276  // then try to remove colinear points that won't
4277  // show up in the generated plot.
4278  std::vector<Point<spacedim>> line_points;
4279  // transform_real_to_unit_cell could be replaced
4280  // by using QProjector<dim>::project_to_line
4281  // which is not yet implemented
4282  const Point<spacedim>
4283  u0 = mapping->transform_real_to_unit_cell(cell,
4284  v0),
4285  u1 = mapping->transform_real_to_unit_cell(cell,
4286  v1);
4287  for (unsigned int i = 0; i < n_points; ++i)
4288  line_points.push_back(
4289  mapping->transform_unit_to_real_cell(
4290  cell,
4291  (1 - boundary_points[i][0]) * u0 +
4292  boundary_points[i][0] * u1));
4293  internal::remove_colinear_points(line_points);
4294  for (const Point<spacedim> &point : line_points)
4295  out << point << ' ' << cell->level() << ' '
4296  << static_cast<unsigned int>(
4297  cell->material_id())
4298  << '\n';
4299  }
4300  else
4301  out << v0 << ' ' << cell->level() << ' '
4302  << cell->material_id() << '\n'
4303  << v1 << ' ' << cell->level() << ' '
4304  << cell->material_id() << '\n';
4305 
4306  out << '\n' << '\n';
4307  }
4308  }
4309  }
4310  }
4311  }
4312 
4313  // make sure everything now gets to disk
4314  out.flush();
4315 
4316  AssertThrow(out, ExcIO());
4317  }
4318  } // namespace
4319 } // namespace internal
4320 
4321 
4322 
4323 template <int dim, int spacedim>
4324 void
4326  std::ostream & out,
4327  const Mapping<dim, spacedim> * mapping) const
4328 {
4329  internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4330 }
4331 
4332 
4333 
4334 namespace internal
4335 {
4336  namespace
4337  {
4338  struct LineEntry
4339  {
4342  bool colorize;
4343  unsigned int level;
4344  LineEntry(const Point<2> & f,
4345  const Point<2> & s,
4346  const bool c,
4347  const unsigned int l)
4348  : first(f)
4349  , second(s)
4350  , colorize(c)
4351  , level(l)
4352  {}
4353  };
4354 
4355 
4356  void
4357  write_eps(const ::Triangulation<1> &,
4358  std::ostream &,
4359  const Mapping<1> *,
4360  const GridOutFlags::Eps<2> &,
4361  const GridOutFlags::Eps<3> &)
4362  {
4363  Assert(false, ExcNotImplemented());
4364  }
4365 
4366  void
4367  write_eps(const ::Triangulation<1, 2> &,
4368  std::ostream &,
4369  const Mapping<1, 2> *,
4370  const GridOutFlags::Eps<2> &,
4371  const GridOutFlags::Eps<3> &)
4372  {
4373  Assert(false, ExcNotImplemented());
4374  }
4375 
4376  void
4377  write_eps(const ::Triangulation<1, 3> &,
4378  std::ostream &,
4379  const Mapping<1, 3> *,
4380  const GridOutFlags::Eps<2> &,
4381  const GridOutFlags::Eps<3> &)
4382  {
4383  Assert(false, ExcNotImplemented());
4384  }
4385 
4386  void
4387  write_eps(const ::Triangulation<2, 3> &,
4388  std::ostream &,
4389  const Mapping<2, 3> *,
4390  const GridOutFlags::Eps<2> &,
4391  const GridOutFlags::Eps<3> &)
4392  {
4393  Assert(false, ExcNotImplemented());
4394  }
4395 
4396 
4397 
4398  template <int dim, int spacedim>
4399  void
4400  write_eps(const ::Triangulation<dim, spacedim> &tria,
4401  std::ostream & out,
4402  const Mapping<dim, spacedim> * mapping,
4405  {
4406  using LineList = std::list<LineEntry>;
4407 
4408  // We should never get here in 1D since this function is overloaded for
4409  // all dim == 1 cases.
4410  Assert(dim == 2 || dim == 3, ExcInternalError());
4411 
4412  // Copy, with an object slice, something containing the flags common to
4413  // all dimensions in order to avoid the recurring distinctions between
4414  // the different eps_flags present.
4415  const GridOutFlags::EpsFlagsBase eps_flags_base =
4416  dim == 2 ?
4417  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4418  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4419 
4420  AssertThrow(out, ExcIO());
4421  const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4422 
4423  // make up a list of lines by which
4424  // we will construct the triangulation
4425  //
4426  // this part unfortunately is a bit
4427  // dimension dependent, so we have to
4428  // treat every dimension different.
4429  // however, by directly producing
4430  // the lines to be printed, i.e. their
4431  // 2d images, we can later do the
4432  // actual output dimension independent
4433  // again
4434  LineList line_list;
4435 
4436  switch (dim)
4437  {
4438  case 1:
4439  {
4440  Assert(false, ExcInternalError());
4441  break;
4442  }
4443 
4444  case 2:
4445  {
4446  for (const auto &cell : tria.active_cell_iterators())
4447  for (const unsigned int line_no : cell->line_indices())
4448  {
4449  typename ::Triangulation<dim, spacedim>::line_iterator
4450  line = cell->line(line_no);
4451 
4452  // first treat all
4453  // interior lines and
4454  // make up a list of
4455  // them. if curved
4456  // lines shall not be
4457  // supported (i.e. no
4458  // mapping is
4459  // provided), then also
4460  // treat all other
4461  // lines
4462  if (!line->has_children() &&
4463  (mapping == nullptr || !line->at_boundary()))
4464  // one would expect
4465  // make_pair(line->vertex(0),
4466  // line->vertex(1))
4467  // here, but that is
4468  // not dimension
4469  // independent, since
4470  // vertex(i) is
4471  // Point<dim>, but we
4472  // want a Point<2>.
4473  // in fact, whenever
4474  // we're here, the
4475  // vertex is a
4476  // Point<dim>, but
4477  // the compiler does
4478  // not know
4479  // this. hopefully,
4480  // the compiler will
4481  // optimize away this
4482  // little kludge
4483  line_list.emplace_back(
4484  Point<2>(line->vertex(0)(0), line->vertex(0)(1)),
4485  Point<2>(line->vertex(1)(0), line->vertex(1)(1)),
4486  line->user_flag_set(),
4487  cell->level());
4488  }
4489 
4490  // next if we are to treat
4491  // curved boundaries
4492  // specially, then add lines
4493  // to the list consisting of
4494  // pieces of the boundary
4495  // lines
4496  if (mapping != nullptr)
4497  {
4498  // to do so, first
4499  // generate a sequence of
4500  // points on a face and
4501  // project them onto the
4502  // faces of a unit cell
4503  std::vector<Point<dim - 1>> boundary_points(n_points);
4504 
4505  for (unsigned int i = 0; i < n_points; ++i)
4506  boundary_points[i](0) = 1. * (i + 1) / (n_points + 1);
4507 
4508  Quadrature<dim - 1> quadrature(boundary_points);
4509  Quadrature<dim> q_projector(
4511 
4512  // next loop over all
4513  // boundary faces and
4514  // generate the info from
4515  // them
4516  for (const auto &cell : tria.active_cell_iterators())
4517  for (const unsigned int face_no :
4519  {
4520  const typename ::Triangulation<dim, spacedim>::
4521  face_iterator face = cell->face(face_no);
4522 
4523  if (face->at_boundary())
4524  {
4525  Point<dim> p0_dim(face->vertex(0));
4526  Point<2> p0(p0_dim(0), p0_dim(1));
4527 
4528  // loop over
4529  // all pieces
4530  // of the line
4531  // and generate
4532  // line-lets
4533  const unsigned int offset = face_no * n_points;
4534  for (unsigned int i = 0; i < n_points; ++i)
4535  {
4536  const Point<dim> p1_dim(
4537  mapping->transform_unit_to_real_cell(
4538  cell, q_projector.point(offset + i)));
4539  const Point<2> p1(p1_dim(0), p1_dim(1));
4540 
4541  line_list.emplace_back(p0,
4542  p1,
4543  face->user_flag_set(),
4544  cell->level());
4545  p0 = p1;
4546  }
4547 
4548  // generate last piece
4549  const Point<dim> p1_dim(face->vertex(1));
4550  const Point<2> p1(p1_dim(0), p1_dim(1));
4551  line_list.emplace_back(p0,
4552  p1,
4553  face->user_flag_set(),
4554  cell->level());
4555  }
4556  }
4557  }
4558 
4559  break;
4560  }
4561 
4562  case 3:
4563  {
4564  // curved boundary output
4565  // presently not supported
4566  Assert(mapping == nullptr, ExcNotImplemented());
4567 
4568  // loop over all lines and compute their
4569  // projection on the plane perpendicular
4570  // to the direction of sight
4571 
4572  // direction of view equals the unit
4573  // vector of the position of the
4574  // spectator to the origin.
4575  //
4576  // we chose here the viewpoint as in
4577  // gnuplot as default.
4578  //
4579  // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4580  // GridOut
4581  // note: the following might be wrong
4582  // if one of the base vectors below
4583  // is in direction of the viewer, but
4584  // I am too tired at present to fix
4585  // this
4586  const double pi = numbers::PI;
4587  const double z_angle = eps_flags_3.azimut_angle;
4588  const double turn_angle = eps_flags_3.turn_angle;
4589  const Point<dim> view_direction(
4590  -std::sin(z_angle * 2. * pi / 360.) *
4591  std::sin(turn_angle * 2. * pi / 360.),
4592  +std::sin(z_angle * 2. * pi / 360.) *
4593  std::cos(turn_angle * 2. * pi / 360.),
4594  -std::cos(z_angle * 2. * pi / 360.));
4595 
4596  // decide about the two unit vectors
4597  // in this plane. we chose the first one
4598  // to be the projection of the z-axis
4599  // to this plane
4600  const Tensor<1, dim> vector1 =
4601  Point<dim>(0, 0, 1) -
4602  ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4603  const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4604 
4605  // now the third vector is fixed. we
4606  // chose the projection of a more or
4607  // less arbitrary vector to the plane
4608  // perpendicular to the first one
4609  const Tensor<1, dim> vector2 =
4610  (Point<dim>(1, 0, 0) -
4611  ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4612  ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4613  const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4614 
4615 
4616  for (const auto &cell : tria.active_cell_iterators())
4617  for (const unsigned int line_no : cell->line_indices())
4618  {
4619  typename ::Triangulation<dim, spacedim>::line_iterator
4620  line = cell->line(line_no);
4621  line_list.emplace_back(
4622  Point<2>(line->vertex(0) * unit_vector2,
4623  line->vertex(0) * unit_vector1),
4624  Point<2>(line->vertex(1) * unit_vector2,
4625  line->vertex(1) * unit_vector1),
4626  line->user_flag_set(),
4627  cell->level());
4628  }
4629 
4630  break;
4631  }
4632 
4633  default:
4634  Assert(false, ExcNotImplemented());
4635  }
4636 
4637 
4638 
4639  // find out minimum and maximum x and
4640  // y coordinates to compute offsets
4641  // and scaling factors
4642  double x_min = tria.begin_active()->vertex(0)(0);
4643  double x_max = x_min;
4644  double y_min = tria.begin_active()->vertex(0)(1);
4645  double y_max = y_min;
4646  unsigned int max_level = line_list.begin()->level;
4647 
4648  for (LineList::const_iterator line = line_list.begin();
4649  line != line_list.end();
4650  ++line)
4651  {
4652  x_min = std::min(x_min, line->first(0));
4653  x_min = std::min(x_min, line->second(0));
4654 
4655  x_max = std::max(x_max, line->first(0));
4656  x_max = std::max(x_max, line->second(0));
4657 
4658  y_min = std::min(y_min, line->first(1));
4659  y_min = std::min(y_min, line->second(1));
4660 
4661  y_max = std::max(y_max, line->first(1));
4662  y_max = std::max(y_max, line->second(1));
4663 
4664  max_level = std::max(max_level, line->level);
4665  }
4666 
4667  // scale in x-direction such that
4668  // in the output 0 <= x <= 300.
4669  // don't scale in y-direction to
4670  // preserve the shape of the
4671  // triangulation
4672  const double scale =
4673  (eps_flags_base.size /
4674  (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4675  x_max - x_min :
4676  y_min - y_max));
4677 
4678 
4679  // now write preamble
4680  {
4681  // block this to have local
4682  // variables destroyed after
4683  // use
4684  std::time_t time1 = std::time(nullptr);
4685  std::tm * time = std::localtime(&time1);
4686  out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4687  << "%%Title: deal.II Output" << '\n'
4688  << "%%Creator: the deal.II library" << '\n'
4689  << "%%Creation Date: " << time->tm_year + 1900 << "/"
4690  << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4691  << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4692  << std::setw(2) << time->tm_sec << '\n'
4693  << "%%BoundingBox: "
4694  // lower left corner
4695  << "0 0 "
4696  // upper right corner
4697  << static_cast<unsigned int>(
4698  std::floor(((x_max - x_min) * scale) + 1))
4699  << ' '
4700  << static_cast<unsigned int>(
4701  std::floor(((y_max - y_min) * scale) + 1))
4702  << '\n';
4703 
4704  // define some abbreviations to keep
4705  // the output small:
4706  // m=move turtle to
4707  // x=execute line stroke
4708  // b=black pen
4709  // r=red pen
4710  out << "/m {moveto} bind def" << '\n'
4711  << "/x {lineto stroke} bind def" << '\n'
4712  << "/b {0 0 0 setrgbcolor} def" << '\n'
4713  << "/r {1 0 0 setrgbcolor} def" << '\n';
4714 
4715  // calculate colors for level
4716  // coloring; level 0 is black,
4717  // other levels are blue
4718  // ... red
4719  if (eps_flags_base.color_lines_level)
4720  out << "/l { neg " << (max_level) << " add "
4721  << (0.66666 / std::max(1U, (max_level - 1)))
4722  << " mul 1 0.8 sethsbcolor} def" << '\n';
4723 
4724  // in 2d, we can also plot cell
4725  // and vertex numbers, but this
4726  // requires a somewhat more
4727  // lengthy preamble. please
4728  // don't ask me what most of
4729  // this means, it is reverse
4730  // engineered from what GNUPLOT
4731  // uses in its output
4732  if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
4733  eps_flags_2.write_vertex_numbers))
4734  {
4735  out
4736  << ("/R {rmoveto} bind def\n"
4737  "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
4738  "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
4739  "currentdict end definefont\n"
4740  "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
4741  "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
4742  "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
4743  "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
4744  "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
4745  "5 get stringwidth pop add}\n"
4746  "{pop} ifelse} forall} bind def\n"
4747  "/MCshow { currentpoint stroke m\n"
4748  "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
4749  << '\n';
4750  }
4751 
4752  out << "%%EndProlog" << '\n' << '\n';
4753 
4754  // set fine lines
4755  out << eps_flags_base.line_width << " setlinewidth" << '\n';
4756  }
4757 
4758  // now write the lines
4759  const Point<2> offset(x_min, y_min);
4760 
4761  for (LineList::const_iterator line = line_list.begin();
4762  line != line_list.end();
4763  ++line)
4764  if (eps_flags_base.color_lines_level && (line->level > 0))
4765  out << line->level << " l " << (line->first - offset) * scale << " m "
4766  << (line->second - offset) * scale << " x" << '\n';
4767  else
4768  out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
4769  "r " :
4770  "b ")
4771  << (line->first - offset) * scale << " m "
4772  << (line->second - offset) * scale << " x" << '\n';
4773 
4774  // finally write the cell numbers
4775  // in 2d, if that is desired
4776  if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
4777  {
4778  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
4779 
4780  for (const auto &cell : tria.active_cell_iterators())
4781  {
4782  out << (cell->center()(0) - offset(0)) * scale << ' '
4783  << (cell->center()(1) - offset(1)) * scale << " m" << '\n'
4784  << "[ [(Helvetica) 12.0 0.0 true true (";
4785  if (eps_flags_2.write_cell_number_level)
4786  out << cell;
4787  else
4788  out << cell->index();
4789 
4790  out << ")] "
4791  << "] -6 MCshow" << '\n';
4792  }
4793  }
4794 
4795  // and the vertex numbers
4796  if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
4797  {
4798  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
4799 
4800  // have a list of those
4801  // vertices which we have
4802  // already tracked, to avoid
4803  // doing this multiply
4804  std::set<unsigned int> treated_vertices;
4805  for (const auto &cell : tria.active_cell_iterators())
4806  for (const unsigned int vertex_no : cell->vertex_indices())
4807  if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
4808  treated_vertices.end())
4809  {
4810  treated_vertices.insert(cell->vertex_index(vertex_no));
4811 
4812  out << (cell->vertex(vertex_no)(0) - offset(0)) * scale << ' '
4813  << (cell->vertex(vertex_no)(1) - offset(1)) * scale
4814  << " m" << '\n'
4815  << "[ [(Helvetica) 10.0 0.0 true true ("
4816  << cell->vertex_index(vertex_no) << ")] "
4817  << "] -6 MCshow" << '\n';
4818  }
4819  }
4820 
4821  out << "showpage" << '\n';
4822 
4823  // make sure everything now gets to
4824  // disk
4825  out.flush();
4826 
4827  AssertThrow(out, ExcIO());
4828  }
4829  } // namespace
4830 } // namespace internal
4831 
4832 
4833 template <int dim, int spacedim>
4834 void
4836  std::ostream & out,
4837  const Mapping<dim, spacedim> * mapping) const
4838 {
4839  internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
4840 }
4841 
4842 
4843 template <int dim, int spacedim>
4844 void
4846  std::ostream & out,
4847  const OutputFormat output_format,
4848  const Mapping<dim, spacedim> * mapping) const
4849 {
4850  switch (output_format)
4851  {
4852  case none:
4853  return;
4854 
4855  case dx:
4856  write_dx(tria, out);
4857  return;
4858 
4859  case ucd:
4860  write_ucd(tria, out);
4861  return;
4862 
4863  case gnuplot:
4864  write_gnuplot(tria, out, mapping);
4865  return;
4866 
4867  case eps:
4868  write_eps(tria, out, mapping);
4869  return;
4870 
4871  case xfig:
4872  write_xfig(tria, out, mapping);
4873  return;
4874 
4875  case msh:
4876  write_msh(tria, out);
4877  return;
4878 
4879  case svg:
4880  write_svg(tria, out);
4881  return;
4882 
4883  case mathgl:
4884  write_mathgl(tria, out);
4885  return;
4886 
4887  case vtk:
4888  write_vtk(tria, out);
4889  return;
4890 
4891  case vtu:
4892  write_vtu(tria, out);
4893  return;
4894  }
4895 
4896  Assert(false, ExcInternalError());
4897 }
4898 
4899 
4900 template <int dim, int spacedim>
4901 void
4903  std::ostream & out,
4904  const Mapping<dim, spacedim> * mapping) const
4905 {
4906  write(tria, out, default_format, mapping);
4907 }
4908 
4909 
4910 // explicit instantiations
4911 #include "grid_out.inst"
4912 
4913 
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:3578
unsigned int n_active_cells() const
Definition: tria.cc:12036
long int get_integer(const std::string &entry_string) const
unsigned int n_used_vertices() const
Definition: tria.cc:12571
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:951
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:11546
OutputFormat
Definition: grid_out.h:1000
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3350
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:3397
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:11374
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:1576
bool write_diameter
Definition: grid_out.h:71
Point< 2 > second
Definition: grid_out.cc:4341
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:2801
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:11354
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:2441
unsigned int write_msh_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3698
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:4325
cell_iterator end() const
Definition: tria.cc:11440
unsigned int n_boundary_lines(const Triangulation< dim, spacedim > &tria) const
Definition: grid_out.cc:3594
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:11557
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:10665
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:1466
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:372
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:647
unsigned int level
Definition: grid_out.cc:4343
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:3131
IteratorRange< active_face_iterator > active_face_iterators() const
Definition: tria.cc:11645
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:339
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:11979
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4340
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:3859
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:12580
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
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:3901
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:3740
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:4835
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:4342
IteratorRange< cell_iterator > cell_iterators() const
Definition: tria.cc:11537
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:4845
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:12647
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