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