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