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