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