Reference documentation for deal.II version Git d3aed38b93 2021-10-28 13:33:27 +0200
\(\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 - 2021 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  else
969  {}
970 
971  // The wrapper
972  out << "object \"deal data\" class field" << '\n'
973  << "component \"positions\" value \"vertices\"" << '\n'
974  << "component \"connections\" value \"cells\"" << '\n';
975 
976  if (write_cells)
977  {
978  out << "object \"cell data\" class field" << '\n'
979  << "component \"positions\" value \"vertices\"" << '\n'
980  << "component \"connections\" value \"cells\"" << '\n';
981  out << "component \"material\" value \"material\"" << '\n';
982  out << "component \"level\" value \"level\"" << '\n';
984  out << "component \"measure\" value \"measure\"" << '\n';
986  out << "component \"diameter\" value \"diameter\"" << '\n';
987  }
988 
989  if (write_faces)
990  {
991  out << "object \"face data\" class field" << '\n'
992  << "component \"positions\" value \"vertices\"" << '\n'
993  << "component \"connections\" value \"faces\"" << '\n';
994  out << "component \"boundary\" value \"boundary\"" << '\n';
996  out << "component \"measure\" value \"face measure\"" << '\n';
998  out << "component \"diameter\" value \"face diameter\"" << '\n';
999  }
1000 
1001  out << '\n' << "object \"grid data\" class group" << '\n';
1002  if (write_cells)
1003  out << "member \"cells\" value \"cell data\"" << '\n';
1004  if (write_faces)
1005  out << "member \"faces\" value \"face data\"" << '\n';
1006  out << "end" << '\n';
1007 
1008  // make sure everything now gets to
1009  // disk
1010  out.flush();
1011 
1012  AssertThrow(out, ExcIO());
1013 }
1014 
1015 
1016 
1017 template <int dim, int spacedim>
1018 void
1020  std::ostream & out) const
1021 {
1022  AssertThrow(out, ExcIO());
1023 
1024  // get the positions of the
1025  // vertices and whether they are
1026  // used.
1027  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1028  const std::vector<bool> & vertex_used = tria.get_used_vertices();
1029 
1030  const unsigned int n_vertices = tria.n_used_vertices();
1031 
1032  // Write Header
1033  // The file format is:
1034  /*
1035 
1036 
1037  @f$NOD
1038  number-of-nodes
1039  node-number x-coord y-coord z-coord
1040  ...
1041  @f$ENDNOD
1042  @f$ELM
1043  number-of-elements
1044  elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
1045  ...
1046  @f$ENDELM
1047  */
1048  out << "@f$NOD" << '\n' << n_vertices << '\n';
1049 
1050  // actually write the vertices.
1051  // note that we shall number them
1052  // with first index 1 instead of 0
1053  for (unsigned int i = 0; i < vertices.size(); ++i)
1054  if (vertex_used[i])
1055  {
1056  out << i + 1 // vertex index
1057  << " " << vertices[i];
1058  for (unsigned int d = spacedim + 1; d <= 3; ++d)
1059  out << " 0"; // fill with zeroes
1060  out << '\n';
1061  }
1062 
1063  // Write cells preamble
1064  out << "@f$ENDNOD" << '\n'
1065  << "@f$ELM" << '\n'
1066  << tria.n_active_cells() +
1067  ((msh_flags.write_faces ? n_boundary_faces(tria) : 0) +
1068  (msh_flags.write_lines ? n_boundary_lines(tria) : 0))
1069  << '\n';
1070 
1071  /*
1072  elm-type
1073  defines the geometrical type of the n-th element:
1074  1
1075  Line (2 nodes).
1076  2
1077  Triangle (3 nodes).
1078  3
1079  Quadrangle (4 nodes).
1080  4
1081  Tetrahedron (4 nodes).
1082  5
1083  Hexahedron (8 nodes).
1084  6
1085  Prism (6 nodes).
1086  7
1087  Pyramid (5 nodes).
1088  8
1089  Second order line (3 nodes: 2 associated with the vertices and 1 with the
1090  edge).
1091  9
1092  Second order triangle (6 nodes: 3 associated with the vertices and 3 with
1093  the edges). 10 Second order quadrangle (9 nodes: 4 associated with the
1094  vertices, 4 with the edges and 1 with the face). 11 Second order tetrahedron
1095  (10 nodes: 4 associated with the vertices and 6 with the edges). 12 Second
1096  order hexahedron (27 nodes: 8 associated with the vertices, 12 with the
1097  edges, 6 with the faces and 1 with the volume). 13 Second order prism (18
1098  nodes: 6 associated with the vertices, 9 with the edges and 3 with the
1099  quadrangular faces). 14 Second order pyramid (14 nodes: 5 associated with
1100  the vertices, 8 with the edges and 1 with the quadrangular face). 15 Point
1101  (1 node).
1102  */
1103  unsigned int elm_type;
1104  switch (dim)
1105  {
1106  case 1:
1107  elm_type = 1;
1108  break;
1109  case 2:
1110  elm_type = 3;
1111  break;
1112  case 3:
1113  elm_type = 5;
1114  break;
1115  default:
1116  Assert(false, ExcNotImplemented());
1117  }
1118 
1119  // write cells. Enumerate cells
1120  // consecutively, starting with 1
1121  for (const auto &cell : tria.active_cell_iterators())
1122  {
1123  out << cell->active_cell_index() + 1 << ' ' << elm_type << ' '
1124  << cell->material_id() << ' ' << cell->subdomain_id() << ' '
1125  << cell->n_vertices() << ' ';
1126 
1127  // Vertex numbering follows UCD conventions.
1128 
1129  for (const unsigned int vertex : GeometryInfo<dim>::vertex_indices())
1130  out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) + 1
1131  << ' ';
1132  out << '\n';
1133  }
1134 
1135  // write faces and lines with non-zero boundary indicator
1136  unsigned int next_element_index = tria.n_active_cells() + 1;
1137  if (msh_flags.write_faces)
1138  {
1139  next_element_index = write_msh_faces(tria, next_element_index, out);
1140  }
1141  if (msh_flags.write_lines)
1142  {
1143  next_element_index = write_msh_lines(tria, next_element_index, out);
1144  }
1145 
1146  out << "@f$ENDELM\n";
1147 
1148  // make sure everything now gets to
1149  // disk
1150  out.flush();
1151 
1152  AssertThrow(out, ExcIO());
1153 }
1154 
1155 
1156 template <int dim, int spacedim>
1157 void
1159  std::ostream & out) const
1160 {
1161  AssertThrow(out, ExcIO());
1162 
1163  // get the positions of the
1164  // vertices and whether they are
1165  // used.
1166  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
1167  const std::vector<bool> & vertex_used = tria.get_used_vertices();
1168 
1169  const unsigned int n_vertices = tria.n_used_vertices();
1170 
1171  // write preamble
1173  {
1174  // block this to have local
1175  // variables destroyed after
1176  // use
1177  std::time_t time1 = std::time(nullptr);
1178  std::tm * time = std::localtime(&time1);
1179  out
1180  << "# This file was generated by the deal.II library." << '\n'
1181  << "# Date = " << time->tm_year + 1900 << "/" << time->tm_mon + 1
1182  << "/" << time->tm_mday << '\n'
1183  << "# Time = " << time->tm_hour << ":" << std::setw(2) << time->tm_min
1184  << ":" << std::setw(2) << time->tm_sec << '\n'
1185  << "#" << '\n'
1186  << "# For a description of the UCD format see the AVS Developer's guide."
1187  << '\n'
1188  << "#" << '\n';
1189  }
1190 
1191  // start with ucd data
1192  out << n_vertices << ' '
1193  << tria.n_active_cells() +
1194  ((ucd_flags.write_faces ? n_boundary_faces(tria) : 0) +
1195  (ucd_flags.write_lines ? n_boundary_lines(tria) : 0))
1196  << " 0 0 0" // no data
1197  << '\n';
1198 
1199  // actually write the vertices.
1200  // note that we shall number them
1201  // with first index 1 instead of 0
1202  for (unsigned int i = 0; i < vertices.size(); ++i)
1203  if (vertex_used[i])
1204  {
1205  out << i + 1 // vertex index
1206  << " " << vertices[i];
1207  for (unsigned int d = spacedim + 1; d <= 3; ++d)
1208  out << " 0"; // fill with zeroes
1209  out << '\n';
1210  }
1211 
1212  // write cells. Enumerate cells
1213  // consecutively, starting with 1
1214  for (const auto &cell : tria.active_cell_iterators())
1215  {
1216  out << cell->active_cell_index() + 1 << ' ' << cell->material_id() << ' ';
1217  switch (dim)
1218  {
1219  case 1:
1220  out << "line ";
1221  break;
1222  case 2:
1223  out << "quad ";
1224  break;
1225  case 3:
1226  out << "hex ";
1227  break;
1228  default:
1229  Assert(false, ExcNotImplemented());
1230  }
1231 
1232  // it follows a list of the
1233  // vertices of each cell. in 1d
1234  // this is simply a list of the
1235  // two vertices, in 2d its counter
1236  // clockwise, as usual in this
1237  // library. in 3d, the same applies
1238  // (special thanks to AVS for
1239  // numbering their vertices in a
1240  // way compatible to deal.II!)
1241  //
1242  // technical reference:
1243  // AVS Developer's Guide, Release 4,
1244  // May, 1992, p. E6
1245  //
1246  // note: vertex numbers are 1-base
1247  for (const unsigned int vertex : GeometryInfo<dim>::vertex_indices())
1248  out << cell->vertex_index(GeometryInfo<dim>::ucd_to_deal[vertex]) + 1
1249  << ' ';
1250  out << '\n';
1251  }
1252 
1253  // write faces and lines with non-zero boundary indicator
1254  unsigned int next_element_index = tria.n_active_cells() + 1;
1255  if (ucd_flags.write_faces)
1256  {
1257  next_element_index = write_ucd_faces(tria, next_element_index, out);
1258  }
1259  if (ucd_flags.write_lines)
1260  {
1261  next_element_index = write_ucd_lines(tria, next_element_index, out);
1262  }
1263 
1264  // make sure everything now gets to
1265  // disk
1266  out.flush();
1267 
1268  AssertThrow(out, ExcIO());
1269 }
1270 
1271 
1272 
1273 template <int dim, int spacedim>
1274 void
1276  std::ostream &,
1277  const Mapping<dim, spacedim> *) const
1278 {
1279  Assert(false, ExcNotImplemented());
1280 }
1281 
1282 
1283 // TODO:[GK] Obey parameters
1284 template <>
1285 void
1287  std::ostream & out,
1288  const Mapping<2> * /*mapping*/) const
1289 {
1290  const int dim = 2;
1291  const int spacedim = 2;
1292 
1293  const unsigned int nv = GeometryInfo<dim>::vertices_per_cell;
1294 
1295  // The following text was copied
1296  // from an existing XFig file.
1297  out << "#FIG 3.2\nLandscape\nCenter\nInches" << std::endl
1298  << "A4\n100.00\nSingle"
1299  << std::endl
1300  // Background is transparent
1301  << "-3" << std::endl
1302  << "# generated by deal.II GridOut class" << std::endl
1303  << "# reduce first number to scale up image" << std::endl
1304  << "1200 2" << std::endl;
1305  // Write custom palette
1306  // grey
1307  unsigned int colno = 32;
1308  out << "0 " << colno++ << " #ff0000" << std::endl;
1309  out << "0 " << colno++ << " #ff8000" << std::endl;
1310  out << "0 " << colno++ << " #ffd000" << std::endl;
1311  out << "0 " << colno++ << " #ffff00" << std::endl;
1312  out << "0 " << colno++ << " #c0ff00" << std::endl;
1313  out << "0 " << colno++ << " #80ff00" << std::endl;
1314  out << "0 " << colno++ << " #00f000" << std::endl;
1315  out << "0 " << colno++ << " #00f0c0" << std::endl;
1316  out << "0 " << colno++ << " #00f0ff" << std::endl;
1317  out << "0 " << colno++ << " #00c0ff" << std::endl;
1318  out << "0 " << colno++ << " #0080ff" << std::endl;
1319  out << "0 " << colno++ << " #0040ff" << std::endl;
1320  out << "0 " << colno++ << " #0000c0" << std::endl;
1321  out << "0 " << colno++ << " #5000ff" << std::endl;
1322  out << "0 " << colno++ << " #8000ff" << std::endl;
1323  out << "0 " << colno++ << " #b000ff" << std::endl;
1324  out << "0 " << colno++ << " #ff00ff" << std::endl;
1325  out << "0 " << colno++ << " #ff80ff" << std::endl;
1326  // grey
1327  for (unsigned int i = 0; i < 8; ++i)
1328  out << "0 " << colno++ << " #" << std::hex << 32 * i + 31 << 32 * i + 31
1329  << 32 * i + 31 << std::dec << std::endl;
1330  // green
1331  for (unsigned int i = 1; i < 16; ++i)
1332  out << "0 " << colno++ << " #00" << std::hex << 16 * i + 15 << std::dec
1333  << "00" << std::endl;
1334  // yellow
1335  for (unsigned int i = 1; i < 16; ++i)
1336  out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << 16 * i + 15
1337  << std::dec << "00" << std::endl;
1338  // red
1339  for (unsigned int i = 1; i < 16; ++i)
1340  out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << std::dec
1341  << "0000" << std::endl;
1342  // purple
1343  for (unsigned int i = 1; i < 16; ++i)
1344  out << "0 " << colno++ << " #" << std::hex << 16 * i + 15 << "00"
1345  << 16 * i + 15 << std::dec << std::endl;
1346  // blue
1347  for (unsigned int i = 1; i < 16; ++i)
1348  out << "0 " << colno++ << " #0000" << std::hex << 16 * i + 15 << std::dec
1349  << std::endl;
1350  // cyan
1351  for (unsigned int i = 1; i < 16; ++i)
1352  out << "0 " << colno++ << " #00" << std::hex << 16 * i + 15 << 16 * i + 15
1353  << std::dec << std::endl;
1354 
1355  // We write all cells and cells on
1356  // coarser levels are behind cells
1357  // on finer levels. Level 0
1358  // corresponds to a depth of 900,
1359  // each level subtracting 1
1360  for (const auto &cell : tria.cell_iterators())
1361  {
1362  // If depth is not encoded, write finest level only
1363  if (!xfig_flags.level_depth && !cell->is_active())
1364  continue;
1365  // Code for polygon
1366  out << "2 3 " << xfig_flags.line_style << ' '
1368  // with black line
1369  << " 0 ";
1370  // Fill color
1371  switch (xfig_flags.color_by)
1372  {
1373  // TODO[GK]: Simplify after deprecation period is over
1375  out << cell->material_id() + 32;
1376  break;
1378  out << cell->level() + 8;
1379  break;
1381  out << cell->subdomain_id() + 32;
1382  break;
1384  out << cell->level_subdomain_id() + 32;
1385  break;
1386  default:
1387  Assert(false, ExcInternalError());
1388  }
1389 
1390  // Depth, unused, fill
1391  out << ' '
1392  << (xfig_flags.level_depth ? (900 - cell->level()) :
1393  (900 + cell->material_id()))
1394  << " 0 " << xfig_flags.fill_style
1395  << " 0.0 "
1396  // some style parameters
1397  << " 0 0 -1 0 0 "
1398  // number of points
1399  << nv + 1 << std::endl;
1400 
1401  // For each point, write scaled
1402  // and shifted coordinates
1403  // multiplied by 1200
1404  // (dots/inch)
1405  for (unsigned int k = 0; k <= nv; ++k)
1406  {
1407  const Point<dim> &p =
1408  cell->vertex(GeometryInfo<dim>::ucd_to_deal[k % nv]);
1409  for (unsigned int d = 0; d < static_cast<unsigned int>(dim); ++d)
1410  {
1411  int val = static_cast<int>(1200 * xfig_flags.scaling(d) *
1412  (p(d) - xfig_flags.offset(d)));
1413  out << '\t' << ((d == 0) ? val : -val);
1414  }
1415  out << std::endl;
1416  }
1417  // Now write boundary edges
1418  static const unsigned int face_reorder[4] = {2, 1, 3, 0};
1420  for (const unsigned int f : face_reorder)
1421  {
1422  Triangulation<dim, spacedim>::face_iterator face = cell->face(f);
1423  const types::boundary_id bi = face->boundary_id();
1425  {
1426  // Code for polyline
1427  out << "2 1 "
1428  // with line style and thickness
1429  << xfig_flags.boundary_style << ' '
1430  << xfig_flags.boundary_thickness << ' ' << 1 + bi;
1431  // Fill color
1432  out << " -1 ";
1433  // Depth 100 less than cells
1434  out << (xfig_flags.level_depth ? (800 - cell->level()) :
1435  800 + bi)
1436  // unused, no fill
1437  << " 0 -1 0.0 "
1438  // some style parameters
1439  << " 0 0 -1 0 0 "
1440  // number of points
1441  << GeometryInfo<dim>::vertices_per_face << std::endl;
1442 
1443  // For each point, write scaled
1444  // and shifted coordinates
1445  // multiplied by 1200
1446  // (dots/inch)
1447 
1448  for (unsigned int k = 0;
1449  k < GeometryInfo<dim>::vertices_per_face;
1450  ++k)
1451  {
1452  const Point<dim> &p = face->vertex(k % nv);
1453  for (unsigned int d = 0; d < static_cast<unsigned int>(dim);
1454  ++d)
1455  {
1456  int val =
1457  static_cast<int>(1200 * xfig_flags.scaling(d) *
1458  (p(d) - xfig_flags.offset(d)));
1459  out << '\t' << ((d == 0) ? val : -val);
1460  }
1461  out << std::endl;
1462  }
1463  }
1464  }
1465  }
1466 
1467  // make sure everything now gets to
1468  // disk
1469  out.flush();
1470 
1471  AssertThrow(out, ExcIO());
1472 }
1473 
1474 
1475 
1476 #ifdef DEAL_II_GMSH_WITH_API
1477 template <int dim, int spacedim>
1478 void
1480  const std::string & filename) const
1481 {
1482  // mesh Type renumbering
1483  const std::array<int, 8> dealii_to_gmsh_type = {{15, 1, 2, 3, 4, 7, 6, 5}};
1484 
1485  // Vertex renumbering, by dealii type
1486  const std::array<std::vector<unsigned int>, 8> dealii_to_gmsh = {
1487  {{0},
1488  {{0, 1}},
1489  {{0, 1, 2}},
1490  {{0, 1, 3, 2}},
1491  {{0, 1, 2, 3}},
1492  {{0, 1, 3, 2, 4}},
1493  {{0, 1, 2, 3, 4, 5}},
1494  {{0, 1, 3, 2, 4, 5, 7, 6}}}};
1495 
1496  // Extract all vertices (nodes in gmsh terminology), and store their three
1497  // dimensional coordinates (regardless of dim).
1498  const auto & vertices = tria.get_vertices();
1499  std::vector<double> coords(3 * vertices.size());
1500  std::vector<std::size_t> nodes(vertices.size());
1501 
1502  // Each node has a strictly positive tag. We assign simply its index+1.
1503  std::size_t i = 0;
1504  for (const auto &p : vertices)
1505  {
1506  for (unsigned int d = 0; d < spacedim; ++d)
1507  coords[i * 3 + d] = p[d];
1508  nodes[i] = i + 1;
1509  ++i;
1510  }
1511 
1512  // Construct one entity tag per boundary and manifold id pair.
1513  // We need to be smart here, in order to save some disk space. All cells need
1514  // to be written, but only faces and lines that have non default boundary ids
1515  // and/or manifold ids. We collect them into pairs, and for each unique pair,
1516  // we create a gmsh entity where we store the elements. Pre-count all the
1517  // entities, and make sure we know which pair refers to what entity and
1518  // vice-versa.
1519  using IdPair = std::pair<types::material_id, types::manifold_id>;
1520  std::map<IdPair, int> id_pair_to_entity_tag;
1521  std::vector<IdPair> all_pairs;
1522  {
1523  std::set<IdPair> set_of_pairs;
1524  for (const auto &cell : tria.active_cell_iterators())
1525  {
1526  set_of_pairs.insert({cell->material_id(), cell->manifold_id()});
1527  for (const auto &f : cell->face_iterators())
1528  if (f->manifold_id() != numbers::flat_manifold_id ||
1529  (f->boundary_id() != 0 &&
1530  f->boundary_id() != numbers::internal_face_boundary_id))
1531  set_of_pairs.insert({f->boundary_id(), f->manifold_id()});
1532  if (dim > 2)
1533  for (const auto l : cell->line_indices())
1534  {
1535  const auto &f = cell->line(l);
1536  if (f->manifold_id() != numbers::flat_manifold_id ||
1537  (f->boundary_id() != 0 &&
1538  f->boundary_id() != numbers::internal_face_boundary_id))
1539  set_of_pairs.insert({f->boundary_id(), f->manifold_id()});
1540  }
1541  }
1542  all_pairs = {set_of_pairs.begin(), set_of_pairs.end()};
1543 
1544  int entity = 1;
1545  for (const auto &p : set_of_pairs)
1546  id_pair_to_entity_tag[p] = entity++;
1547  }
1548 
1549  const auto n_entity_tags = id_pair_to_entity_tag.size();
1550 
1551  // All elements in the mesh, by entity tag, and by dealii type.
1552  std::vector<std::vector<std::vector<std::size_t>>> element_ids(
1553  n_entity_tags, std::vector<std::vector<std::size_t>>(8));
1554  std::vector<std::vector<std::vector<std::size_t>>> element_nodes(
1555  n_entity_tags, std::vector<std::vector<std::size_t>>(8));
1556 
1557  // One element id counter for all dimensions.
1558  std::size_t element_id = 1;
1559 
1560  const auto add_element = [&](const auto &element, const int &entity_tag) {
1561  const auto type = element->reference_cell();
1562 
1563  Assert(entity_tag > 0, ExcInternalError());
1564  // Add all vertex ids. Make sure we renumber to gmsh, and we add 1 to the
1565  // global index.
1566  for (const auto v : element->vertex_indices())
1567  element_nodes[entity_tag - 1][type].emplace_back(
1568  element->vertex_index(dealii_to_gmsh[type][v]) + 1);
1569 
1570  // Save the element id.
1571  element_ids[entity_tag - 1][type].emplace_back(element_id);
1572  ++element_id;
1573  };
1574 
1575  // Will create a separate gmsh entity, only if it's a cell, or if the
1576  // boundary and/or the manifold ids are not the default ones.
1577  // In the meanwhile, also store each pair of dimension and entity tag that was
1578  // requested.
1579  std::set<std::pair<int, int>> dim_entity_tag;
1580 
1581  auto maybe_add_element =
1582  [&](const auto & element,
1583  const types::boundary_id &boundary_or_material_id) {
1584  const auto struct_dim = element->structure_dimension;
1585  const auto manifold_id = element->manifold_id();
1586 
1587  // Exclude default boundary/manifold id or invalid/flag
1588  const bool non_default_boundary_or_material_id =
1589  (boundary_or_material_id != 0 &&
1590  boundary_or_material_id != numbers::internal_face_boundary_id);
1591  const bool non_default_manifold =
1593  if (struct_dim == dim || non_default_boundary_or_material_id ||
1594  non_default_manifold)
1595  {
1596  const auto entity_tag =
1597  id_pair_to_entity_tag[{boundary_or_material_id, manifold_id}];
1598  add_element(element, entity_tag);
1599  dim_entity_tag.insert({struct_dim, entity_tag});
1600  }
1601  };
1602 
1603  // Loop recursively over all cells, faces, and possibly lines.
1604  for (const auto &cell : tria.active_cell_iterators())
1605  {
1606  maybe_add_element(cell, cell->material_id());
1607  for (const auto &face : cell->face_iterators())
1608  maybe_add_element(face, face->boundary_id());
1609  if (dim > 2)
1610  for (const auto l : cell->line_indices())
1611  maybe_add_element(cell->line(l), cell->line(l)->boundary_id());
1612  }
1613 
1614  // Now that we collected everything, plug them into gmsh
1615  gmsh::initialize();
1616  gmsh::option::setNumber("General.Verbosity", 0);
1617  gmsh::model::add("Grid generated in deal.II");
1618  for (const auto &p : dim_entity_tag)
1619  {
1620  gmsh::model::addDiscreteEntity(p.first, p.second);
1621  gmsh::model::mesh::addNodes(p.first, p.second, nodes, coords);
1622  }
1623 
1624  for (unsigned int entity_tag = 0; entity_tag < n_entity_tags; ++entity_tag)
1625  for (unsigned int t = 1; t < 8; ++t)
1626  {
1627  const auto all_element_ids = element_ids[entity_tag][t];
1628  const auto all_element_nodes = element_nodes[entity_tag][t];
1629  const auto gmsh_t = dealii_to_gmsh_type[t];
1630  if (all_element_ids.size() > 0)
1631  gmsh::model::mesh::addElementsByType(entity_tag + 1,
1632  gmsh_t,
1633  all_element_ids,
1634  all_element_nodes);
1635  }
1636 
1637  // Make sure nodes belong to the right entities.
1638  gmsh::model::mesh::reclassifyNodes();
1639  gmsh::model::mesh::removeDuplicateNodes();
1640 
1641  // Now for each individual pair of dim and entry, add a physical group, if
1642  // necessary
1643  for (const auto &it : dim_entity_tag)
1644  {
1645  const auto &d = it.first;
1646  const auto &entity_tag = it.second;
1647  const auto &boundary_id = all_pairs[entity_tag - 1].first;
1648  const auto &manifold_id = all_pairs[entity_tag - 1].second;
1649 
1650  std::string physical_name;
1651  if (d == dim && boundary_id != 0)
1652  physical_name += "MaterialID:" + Utilities::int_to_string(
1653  static_cast<int>(boundary_id));
1654  else if (d < dim && boundary_id != 0)
1655  physical_name +=
1656  "BoundaryID:" +
1658  "-1" :
1659  Utilities::int_to_string(static_cast<int>(boundary_id)));
1660 
1661  std::string sep = physical_name != "" ? ", " : "";
1663  physical_name +=
1664  sep + "ManifoldID:" +
1665  Utilities::int_to_string(static_cast<int>(manifold_id));
1666  const auto physical_tag =
1667  gmsh::model::addPhysicalGroup(d, {entity_tag}, -1);
1668  if (physical_name != "")
1669  gmsh::model::setPhysicalName(d, physical_tag, physical_name);
1670  }
1671 
1672 
1673  gmsh::write(filename);
1674  gmsh::clear();
1675  gmsh::finalize();
1676 }
1677 #endif
1678 
1679 
1680 
1681 namespace
1682 {
1693  Point<2>
1694  svg_project_point(const Point<3> & point,
1695  const Point<3> & camera_position,
1696  const Tensor<1, 3> &camera_direction,
1697  const Tensor<1, 3> &camera_horizontal,
1698  const float camera_focus)
1699  {
1700  const Tensor<1, 3> camera_vertical =
1701  cross_product_3d(camera_horizontal, camera_direction);
1702 
1703  const float phi =
1704  camera_focus / ((point - camera_position) * camera_direction);
1705 
1706  const Point<3> projection =
1707  camera_position + phi * (point - camera_position);
1708 
1709  return {(projection - camera_position - camera_focus * camera_direction) *
1710  camera_horizontal,
1711  (projection - camera_position - camera_focus * camera_direction) *
1712  camera_vertical};
1713  }
1714 } // namespace
1715 
1716 
1717 
1718 template <int dim, int spacedim>
1719 void
1721  std::ostream & /*out*/) const
1722 {
1723  Assert(false,
1724  ExcMessage("Mesh output in SVG format is not implemented for anything "
1725  "other than two-dimensional meshes in two-dimensional "
1726  "space. That's because three-dimensional meshes are best "
1727  "viewed in programs that allow changing the viewpoint, "
1728  "but SVG format does not allow this: It is an inherently "
1729  "2d format, and for three-dimensional meshes would "
1730  "require choosing one, fixed viewpoint."
1731  "\n\n"
1732  "You probably want to output your mesh in a format such "
1733  "as VTK, VTU, or gnuplot."));
1734 }
1735 
1736 
1737 void
1738 GridOut::write_svg(const Triangulation<2, 2> &tria, std::ostream &out) const
1739 {
1740  unsigned int n = 0;
1741 
1742  unsigned int min_level, max_level;
1743 
1744  // Svg files require an underlying drawing grid. The size of this
1745  // grid is provided in the parameters height and width. Each of them
1746  // may be zero, such that it is computed from the other. Obviously,
1747  // both of them zero does not produce reasonable output.
1748  unsigned int height = svg_flags.height;
1749  unsigned int width = svg_flags.width;
1750  Assert(height != 0 || width != 0,
1751  ExcMessage("You have to set at least one of width and height"));
1752 
1753  unsigned int margin_in_percent = 0;
1755  margin_in_percent = 8;
1756 
1757  // initial font size for cell labels
1758  unsigned int cell_label_font_size;
1759 
1760  // get date and time
1761  // time_t time_stamp;
1762  // tm *now;
1763  // time_stamp = time(0);
1764  // now = localtime(&time_stamp);
1765 
1766  float camera_focus;
1767 
1768  Point<3> point;
1769  Point<2> projection_decomposition;
1770 
1771  float x_max_perspective, x_min_perspective;
1772  float y_max_perspective, y_min_perspective;
1773 
1774  float x_dimension_perspective, y_dimension_perspective;
1775 
1776 
1777  // auxiliary variables for the bounding box and the range of cell levels
1778  double x_min = tria.begin()->vertex(0)[0];
1779  double x_max = x_min;
1780  double y_min = tria.begin()->vertex(0)[1];
1781  double y_max = y_min;
1782 
1783  double x_dimension, y_dimension;
1784 
1785  min_level = max_level = tria.begin()->level();
1786 
1787  // auxiliary set for the materials being used
1788  std::set<unsigned int> materials;
1789 
1790  // auxiliary set for the levels being used
1791  std::set<unsigned int> levels;
1792 
1793  // auxiliary set for the subdomains being used
1794  std::set<unsigned int> subdomains;
1795 
1796  // auxiliary set for the level subdomains being used
1797  std::set<int> level_subdomains;
1798 
1799  // We use an active cell iterator to determine the
1800  // bounding box of the given triangulation and check
1801  // the cells for material id, level number, subdomain id
1802  // (, and level subdomain id).
1803  for (const auto &cell : tria.cell_iterators())
1804  {
1805  for (unsigned int vertex_index = 0; vertex_index < cell->n_vertices();
1806  ++vertex_index)
1807  {
1808  if (cell->vertex(vertex_index)[0] < x_min)
1809  x_min = cell->vertex(vertex_index)[0];
1810  if (cell->vertex(vertex_index)[0] > x_max)
1811  x_max = cell->vertex(vertex_index)[0];
1812 
1813  if (cell->vertex(vertex_index)[1] < y_min)
1814  y_min = cell->vertex(vertex_index)[1];
1815  if (cell->vertex(vertex_index)[1] > y_max)
1816  y_max = cell->vertex(vertex_index)[1];
1817  }
1818 
1819  if (static_cast<unsigned int>(cell->level()) < min_level)
1820  min_level = cell->level();
1821  if (static_cast<unsigned int>(cell->level()) > max_level)
1822  max_level = cell->level();
1823 
1824  materials.insert(cell->material_id());
1825  levels.insert(cell->level());
1826  if (cell->is_active())
1827  subdomains.insert(cell->subdomain_id() + 2);
1828  level_subdomains.insert(cell->level_subdomain_id() + 2);
1829  }
1830 
1831  x_dimension = x_max - x_min;
1832  y_dimension = y_max - y_min;
1833 
1834  // count the materials being used
1835  const unsigned int n_materials = materials.size();
1836 
1837  // count the levels being used
1838  const unsigned int n_levels = levels.size();
1839 
1840  // count the subdomains being used
1841  const unsigned int n_subdomains = subdomains.size();
1842 
1843  // count the level subdomains being used
1844  const unsigned int n_level_subdomains = level_subdomains.size();
1845 
1846  switch (svg_flags.coloring)
1847  {
1849  n = n_materials;
1850  break;
1852  n = n_levels;
1853  break;
1855  n = n_subdomains;
1856  break;
1858  n = n_level_subdomains;
1859  break;
1860  default:
1861  break;
1862  }
1863 
1864  // set the camera position to top view, targeting at the origin
1865  // vectors and variables for the perspective view
1866  Point<3> camera_position;
1867  camera_position[0] = 0;
1868  camera_position[1] = 0;
1869  camera_position[2] = 2. * std::max(x_dimension, y_dimension);
1870 
1871  Tensor<1, 3> camera_direction;
1872  camera_direction[0] = 0;
1873  camera_direction[1] = 0;
1874  camera_direction[2] = -1;
1875 
1876  Tensor<1, 3> camera_horizontal;
1877  camera_horizontal[0] = 1;
1878  camera_horizontal[1] = 0;
1879  camera_horizontal[2] = 0;
1880 
1881  camera_focus = .5 * std::max(x_dimension, y_dimension);
1882 
1883  Point<3> camera_position_temp;
1884  Point<3> camera_direction_temp;
1885  Point<3> camera_horizontal_temp;
1886 
1887  const double angle_factor = 3.14159265 / 180.;
1888 
1889  // (I) rotate the camera to the chosen polar angle
1890  camera_position_temp[1] =
1891  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[1] -
1892  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[2];
1893  camera_position_temp[2] =
1894  std::sin(angle_factor * svg_flags.polar_angle) * camera_position[1] +
1895  std::cos(angle_factor * svg_flags.polar_angle) * camera_position[2];
1896 
1897  camera_direction_temp[1] =
1898  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[1] -
1899  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1900  camera_direction_temp[2] =
1901  std::sin(angle_factor * svg_flags.polar_angle) * camera_direction[1] +
1902  std::cos(angle_factor * svg_flags.polar_angle) * camera_direction[2];
1903 
1904  camera_horizontal_temp[1] =
1905  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] -
1906  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1907  camera_horizontal_temp[2] =
1908  std::sin(angle_factor * svg_flags.polar_angle) * camera_horizontal[1] +
1909  std::cos(angle_factor * svg_flags.polar_angle) * camera_horizontal[2];
1910 
1911  camera_position[1] = camera_position_temp[1];
1912  camera_position[2] = camera_position_temp[2];
1913 
1914  camera_direction[1] = camera_direction_temp[1];
1915  camera_direction[2] = camera_direction_temp[2];
1916 
1917  camera_horizontal[1] = camera_horizontal_temp[1];
1918  camera_horizontal[2] = camera_horizontal_temp[2];
1919 
1920  // (II) rotate the camera to the chosen azimuth angle
1921  camera_position_temp[0] =
1922  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[0] -
1923  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1924  camera_position_temp[1] =
1925  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_position[0] +
1926  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_position[1];
1927 
1928  camera_direction_temp[0] =
1929  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] -
1930  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1931  camera_direction_temp[1] =
1932  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_direction[0] +
1933  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_direction[1];
1934 
1935  camera_horizontal_temp[0] =
1936  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] -
1937  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1938  camera_horizontal_temp[1] =
1939  std::sin(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[0] +
1940  std::cos(angle_factor * svg_flags.azimuth_angle) * camera_horizontal[1];
1941 
1942  camera_position[0] = camera_position_temp[0];
1943  camera_position[1] = camera_position_temp[1];
1944 
1945  camera_direction[0] = camera_direction_temp[0];
1946  camera_direction[1] = camera_direction_temp[1];
1947 
1948  camera_horizontal[0] = camera_horizontal_temp[0];
1949  camera_horizontal[1] = camera_horizontal_temp[1];
1950 
1951  // translate the camera to the given triangulation
1952  camera_position[0] = x_min + .5 * x_dimension;
1953  camera_position[1] = y_min + .5 * y_dimension;
1954 
1955  camera_position[0] += 2. * std::max(x_dimension, y_dimension) *
1956  std::sin(angle_factor * svg_flags.polar_angle) *
1957  std::sin(angle_factor * svg_flags.azimuth_angle);
1958  camera_position[1] -= 2. * std::max(x_dimension, y_dimension) *
1959  std::sin(angle_factor * svg_flags.polar_angle) *
1960  std::cos(angle_factor * svg_flags.azimuth_angle);
1961 
1962 
1963  // determine the bounding box of the given triangulation on the projection
1964  // plane of the camera viewing system
1965  point[0] = tria.begin()->vertex(0)[0];
1966  point[1] = tria.begin()->vertex(0)[1];
1967  point[2] = 0;
1968 
1969  float min_level_min_vertex_distance = 0;
1970 
1972  {
1973  point[2] = svg_flags.level_height_factor *
1974  (static_cast<float>(tria.begin()->level()) /
1975  static_cast<float>(n_levels)) *
1976  std::max(x_dimension, y_dimension);
1977  }
1978 
1979  projection_decomposition = svg_project_point(
1980  point, camera_position, camera_direction, camera_horizontal, camera_focus);
1981 
1982  x_max_perspective = projection_decomposition[0];
1983  x_min_perspective = projection_decomposition[0];
1984 
1985  y_max_perspective = projection_decomposition[1];
1986  y_min_perspective = projection_decomposition[1];
1987 
1988  for (const auto &cell : tria.cell_iterators())
1989  {
1990  point[0] = cell->vertex(0)[0];
1991  point[1] = cell->vertex(0)[1];
1992  point[2] = 0;
1993 
1995  {
1996  point[2] =
1998  (static_cast<float>(cell->level()) / static_cast<float>(n_levels)) *
1999  std::max(x_dimension, y_dimension);
2000  }
2001 
2002  projection_decomposition = svg_project_point(point,
2003  camera_position,
2004  camera_direction,
2005  camera_horizontal,
2006  camera_focus);
2007 
2008  if (x_max_perspective < projection_decomposition[0])
2009  x_max_perspective = projection_decomposition[0];
2010  if (x_min_perspective > projection_decomposition[0])
2011  x_min_perspective = projection_decomposition[0];
2012 
2013  if (y_max_perspective < projection_decomposition[1])
2014  y_max_perspective = projection_decomposition[1];
2015  if (y_min_perspective > projection_decomposition[1])
2016  y_min_perspective = projection_decomposition[1];
2017 
2018  point[0] = cell->vertex(1)[0];
2019  point[1] = cell->vertex(1)[1];
2020 
2021  projection_decomposition = svg_project_point(point,
2022  camera_position,
2023  camera_direction,
2024  camera_horizontal,
2025  camera_focus);
2026 
2027  if (x_max_perspective < projection_decomposition[0])
2028  x_max_perspective = projection_decomposition[0];
2029  if (x_min_perspective > projection_decomposition[0])
2030  x_min_perspective = projection_decomposition[0];
2031 
2032  if (y_max_perspective < projection_decomposition[1])
2033  y_max_perspective = projection_decomposition[1];
2034  if (y_min_perspective > projection_decomposition[1])
2035  y_min_perspective = projection_decomposition[1];
2036 
2037  point[0] = cell->vertex(2)[0];
2038  point[1] = cell->vertex(2)[1];
2039 
2040  projection_decomposition = svg_project_point(point,
2041  camera_position,
2042  camera_direction,
2043  camera_horizontal,
2044  camera_focus);
2045 
2046  if (x_max_perspective < projection_decomposition[0])
2047  x_max_perspective = projection_decomposition[0];
2048  if (x_min_perspective > projection_decomposition[0])
2049  x_min_perspective = projection_decomposition[0];
2050 
2051  if (y_max_perspective < projection_decomposition[1])
2052  y_max_perspective = projection_decomposition[1];
2053  if (y_min_perspective > projection_decomposition[1])
2054  y_min_perspective = projection_decomposition[1];
2055 
2056  if (cell->n_vertices() == 4) // in case of quadrilateral
2057  {
2058  point[0] = cell->vertex(3)[0];
2059  point[1] = cell->vertex(3)[1];
2060 
2061  projection_decomposition = svg_project_point(point,
2062  camera_position,
2063  camera_direction,
2064  camera_horizontal,
2065  camera_focus);
2066 
2067  if (x_max_perspective < projection_decomposition[0])
2068  x_max_perspective = projection_decomposition[0];
2069  if (x_min_perspective > projection_decomposition[0])
2070  x_min_perspective = projection_decomposition[0];
2071 
2072  if (y_max_perspective < projection_decomposition[1])
2073  y_max_perspective = projection_decomposition[1];
2074  if (y_min_perspective > projection_decomposition[1])
2075  y_min_perspective = projection_decomposition[1];
2076  }
2077 
2078  if (static_cast<unsigned int>(cell->level()) == min_level)
2079  min_level_min_vertex_distance = cell->minimum_vertex_distance();
2080  }
2081 
2082  x_dimension_perspective = x_max_perspective - x_min_perspective;
2083  y_dimension_perspective = y_max_perspective - y_min_perspective;
2084 
2085  // create the svg file with an internal style sheet
2086  if (width == 0)
2087  width = static_cast<unsigned int>(
2088  .5 + height * (x_dimension_perspective / y_dimension_perspective));
2089  else if (height == 0)
2090  height = static_cast<unsigned int>(
2091  .5 + width * (y_dimension_perspective / x_dimension_perspective));
2092  unsigned int additional_width = 0;
2093  // font size for date, time, legend, and colorbar
2094  unsigned int font_size =
2095  static_cast<unsigned int>(.5 + (height / 100.) * 1.75);
2096  cell_label_font_size = static_cast<unsigned int>(
2097  .5 + (height * .15 * svg_flags.cell_font_scaling *
2098  min_level_min_vertex_distance / std::min(x_dimension, y_dimension)));
2099 
2100  if (svg_flags.draw_legend &&
2104  {
2105  additional_width = static_cast<unsigned int>(
2106  .5 + height * .4); // additional width for legend
2107  }
2109  {
2110  additional_width = static_cast<unsigned int>(
2111  .5 + height * .175); // additional width for colorbar
2112  }
2113 
2114  // out << "<!-- deal.ii GridOut " << now->tm_mday << '/' << now->tm_mon + 1 <<
2115  // '/' << now->tm_year + 1900
2116  // << ' ' << now->tm_hour << ':';
2117  //
2118  // if (now->tm_min < 10) out << '0';
2119  //
2120  // out << now->tm_min << " -->" << '\n';
2121 
2122  // basic svg header
2123  out << "<svg width=\"" << width + additional_width << "\" height=\"" << height
2124  << "\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">" << '\n'
2125  << '\n';
2126 
2127 
2129  {
2130  out
2131  << " <linearGradient id=\"background_gradient\" gradientUnits=\"userSpaceOnUse\" x1=\"0\" y1=\"0\" x2=\"0\" y2=\""
2132  << height << "\">" << '\n'
2133  << " <stop offset=\"0\" style=\"stop-color:white\"/>" << '\n'
2134  << " <stop offset=\"1\" style=\"stop-color:lightsteelblue\"/>" << '\n'
2135  << " </linearGradient>" << '\n';
2136  }
2137 
2138  out << '\n';
2139 
2140  // header for the internal style sheet
2141  out << "<!-- internal style sheet -->" << '\n'
2142  << "<style type=\"text/css\"><![CDATA[" << '\n';
2143 
2144  // set the background of the output graphic
2146  out << " rect.background{fill:url(#background_gradient)}" << '\n';
2148  out << " rect.background{fill:white}" << '\n';
2149  else
2150  out << " rect.background{fill:none}" << '\n';
2151 
2152  // basic svg graphic element styles
2153  out << " rect{fill:none; stroke:rgb(25,25,25); stroke-width:"
2154  << svg_flags.line_thickness << '}' << '\n'
2155  << " text{font-family:Helvetica; text-anchor:middle; fill:rgb(25,25,25)}"
2156  << '\n'
2157  << " line{stroke:rgb(25,25,25); stroke-width:"
2158  << svg_flags.boundary_line_thickness << '}' << '\n'
2159  << " path{fill:none; stroke:rgb(25,25,25); stroke-width:"
2160  << svg_flags.line_thickness << '}' << '\n'
2161  << " circle{fill:white; stroke:black; stroke-width:2}" << '\n'
2162  << '\n';
2163 
2164  // polygon styles with respect to the chosen cell coloring
2165  if (svg_flags.coloring)
2166  {
2167  unsigned int labeling_index = 0;
2168  auto materials_it = materials.begin();
2169  auto levels_it = levels.begin();
2170  auto subdomains_it = subdomains.begin();
2171  auto level_subdomains_it = level_subdomains.begin();
2172 
2173  for (unsigned int index = 0; index < n; ++index)
2174  {
2175  double h;
2176 
2177  if (n != 1)
2178  h = .6 - (index / (n - 1.)) * .6;
2179  else
2180  h = .6;
2181 
2182  unsigned int r = 0;
2183  unsigned int g = 0;
2184  unsigned int b = 0;
2185 
2186  unsigned int i = static_cast<unsigned int>(h * 6);
2187 
2188  double f = h * 6 - i;
2189  double q = 1 - f;
2190  double t = f;
2191 
2192  switch (i % 6)
2193  {
2194  case 0:
2195  r = 255, g = static_cast<unsigned int>(.5 + 255 * t);
2196  break;
2197  case 1:
2198  r = static_cast<unsigned int>(.5 + 255 * q), g = 255;
2199  break;
2200  case 2:
2201  g = 255, b = static_cast<unsigned int>(.5 + 255 * t);
2202  break;
2203  case 3:
2204  g = static_cast<unsigned int>(.5 + 255 * q), b = 255;
2205  break;
2206  case 4:
2207  r = static_cast<unsigned int>(.5 + 255 * t), b = 255;
2208  break;
2209  case 5:
2210  r = 255, b = static_cast<unsigned int>(.5 + 255 * q);
2211  break;
2212  default:
2213  break;
2214  }
2215 
2216  switch (svg_flags.coloring)
2217  {
2219  labeling_index = *materials_it++;
2220  break;
2222  labeling_index = *levels_it++;
2223  break;
2225  labeling_index = *subdomains_it++;
2226  break;
2228  labeling_index = *level_subdomains_it++;
2229  break;
2230  default:
2231  break;
2232  }
2233 
2234  out << " path.p" << labeling_index << "{fill:rgb(" << r << ',' << g
2235  << ',' << b << "); "
2236  << "stroke:rgb(25,25,25); stroke-width:"
2237  << svg_flags.line_thickness << '}' << '\n';
2238 
2239  out << " path.ps" << labeling_index << "{fill:rgb("
2240  << static_cast<unsigned int>(.5 + .75 * r) << ','
2241  << static_cast<unsigned int>(.5 + .75 * g) << ','
2242  << static_cast<unsigned int>(.5 + .75 * b) << "); "
2243  << "stroke:rgb(20,20,20); stroke-width:"
2244  << svg_flags.line_thickness << '}' << '\n';
2245 
2246  out << " rect.r" << labeling_index << "{fill:rgb(" << r << ',' << g
2247  << ',' << b << "); "
2248  << "stroke:rgb(25,25,25); stroke-width:"
2249  << svg_flags.line_thickness << '}' << '\n';
2250 
2251  labeling_index++;
2252  }
2253  }
2254 
2255  out << "]]></style>" << '\n' << '\n';
2256 
2257  // background rectangle
2258  out << " <rect class=\"background\" width=\"" << width << "\" height=\""
2259  << height << "\"/>" << '\n';
2260 
2262  {
2263  unsigned int x_offset = 0;
2264 
2265  if (svg_flags.margin)
2266  x_offset = static_cast<unsigned int>(.5 + (height / 100.) *
2267  (margin_in_percent / 2.));
2268  else
2269  x_offset = static_cast<unsigned int>(.5 + height * .025);
2270 
2271  out
2272  << " <text x=\"" << x_offset << "\" y=\""
2273  << static_cast<unsigned int>(.5 + height * .0525) << '\"'
2274  << " style=\"font-weight:100; fill:lightsteelblue; text-anchor:start; font-family:Courier; font-size:"
2275  << static_cast<unsigned int>(.5 + height * .045) << "px\">"
2276  << "deal.II"
2277  << "</text>" << '\n';
2278 
2279  // out << " <text x=\"" << x_offset + static_cast<unsigned int>(.5 +
2280  // height * .045 * 4.75) << "\" y=\"" << static_cast<unsigned int>(.5 +
2281  // height * .0525) << '\"'
2282  // << " style=\"fill:lightsteelblue; text-anchor:start; font-size:" <<
2283  // font_size << "\">"
2284  // << now->tm_mday << '/' << now->tm_mon + 1 << '/' << now->tm_year +
2285  // 1900
2286  // << " - " << now->tm_hour << ':';
2287  //
2288  // if(now->tm_min < 10) out << '0';
2289  //
2290  // out << now->tm_min
2291  // << "</text>"<< '\n' << '\n';
2292  }
2293 
2294  // draw the cells, starting out from the minimal level (in order to guaranty a
2295  // correct perspective view)
2296  out << " <!-- cells -->" << '\n';
2297 
2298  for (unsigned int level_index = min_level; level_index <= max_level;
2299  level_index++)
2300  {
2301  for (const auto &cell : tria.cell_iterators_on_level(level_index))
2302  {
2303  if (!svg_flags.convert_level_number_to_height && !cell->is_active())
2304  continue;
2305 
2306  // draw the current cell
2307  out << " <path";
2308 
2309  if (svg_flags.coloring)
2310  {
2311  out << " class=\"p";
2312 
2313  if (!cell->is_active() &&
2315  out << 's';
2316 
2317  switch (svg_flags.coloring)
2318  {
2320  out << cell->material_id();
2321  break;
2323  out << static_cast<unsigned int>(cell->level());
2324  break;
2326  if (cell->is_active())
2327  out << cell->subdomain_id() + 2;
2328  else
2329  out << 'X';
2330  break;
2332  out << cell->level_subdomain_id() + 2;
2333  break;
2334  default:
2335  break;
2336  }
2337 
2338  out << '\"';
2339  }
2340 
2341  out << " d=\"M ";
2342 
2343  point[0] = cell->vertex(0)[0];
2344  point[1] = cell->vertex(0)[1];
2345  point[2] = 0;
2346 
2348  {
2349  point[2] = svg_flags.level_height_factor *
2350  (static_cast<float>(cell->level()) /
2351  static_cast<float>(n_levels)) *
2352  std::max(x_dimension, y_dimension);
2353  }
2354 
2355  projection_decomposition = svg_project_point(point,
2356  camera_position,
2357  camera_direction,
2358  camera_horizontal,
2359  camera_focus);
2360 
2361  out << static_cast<unsigned int>(
2362  .5 +
2363  ((projection_decomposition[0] - x_min_perspective) /
2364  x_dimension_perspective) *
2365  (width - (width / 100.) * 2. * margin_in_percent) +
2366  ((width / 100.) * margin_in_percent))
2367  << ' '
2368  << static_cast<unsigned int>(
2369  .5 + height - (height / 100.) * margin_in_percent -
2370  ((projection_decomposition[1] - y_min_perspective) /
2371  y_dimension_perspective) *
2372  (height - (height / 100.) * 2. * margin_in_percent));
2373 
2374  out << " L ";
2375 
2376  point[0] = cell->vertex(1)[0];
2377  point[1] = cell->vertex(1)[1];
2378 
2379  projection_decomposition = svg_project_point(point,
2380  camera_position,
2381  camera_direction,
2382  camera_horizontal,
2383  camera_focus);
2384 
2385  out << static_cast<unsigned int>(
2386  .5 +
2387  ((projection_decomposition[0] - x_min_perspective) /
2388  x_dimension_perspective) *
2389  (width - (width / 100.) * 2. * margin_in_percent) +
2390  ((width / 100.) * margin_in_percent))
2391  << ' '
2392  << static_cast<unsigned int>(
2393  .5 + height - (height / 100.) * margin_in_percent -
2394  ((projection_decomposition[1] - y_min_perspective) /
2395  y_dimension_perspective) *
2396  (height - (height / 100.) * 2. * margin_in_percent));
2397 
2398  out << " L ";
2399 
2400  if (cell->n_vertices() == 4) // in case of quadrilateral
2401  {
2402  point[0] = cell->vertex(3)[0];
2403  point[1] = cell->vertex(3)[1];
2404 
2405  projection_decomposition = svg_project_point(point,
2406  camera_position,
2407  camera_direction,
2408  camera_horizontal,
2409  camera_focus);
2410 
2411  out << static_cast<unsigned int>(
2412  .5 +
2413  ((projection_decomposition[0] - x_min_perspective) /
2414  x_dimension_perspective) *
2415  (width - (width / 100.) * 2. * margin_in_percent) +
2416  ((width / 100.) * margin_in_percent))
2417  << ' '
2418  << static_cast<unsigned int>(
2419  .5 + height - (height / 100.) * margin_in_percent -
2420  ((projection_decomposition[1] - y_min_perspective) /
2421  y_dimension_perspective) *
2422  (height - (height / 100.) * 2. * margin_in_percent));
2423 
2424  out << " L ";
2425  }
2426 
2427  point[0] = cell->vertex(2)[0];
2428  point[1] = cell->vertex(2)[1];
2429 
2430  projection_decomposition = svg_project_point(point,
2431  camera_position,
2432  camera_direction,
2433  camera_horizontal,
2434  camera_focus);
2435 
2436  out << static_cast<unsigned int>(
2437  .5 +
2438  ((projection_decomposition[0] - x_min_perspective) /
2439  x_dimension_perspective) *
2440  (width - (width / 100.) * 2. * margin_in_percent) +
2441  ((width / 100.) * margin_in_percent))
2442  << ' '
2443  << static_cast<unsigned int>(
2444  .5 + height - (height / 100.) * margin_in_percent -
2445  ((projection_decomposition[1] - y_min_perspective) /
2446  y_dimension_perspective) *
2447  (height - (height / 100.) * 2. * margin_in_percent));
2448 
2449  out << " L ";
2450 
2451  point[0] = cell->vertex(0)[0];
2452  point[1] = cell->vertex(0)[1];
2453 
2454  projection_decomposition = svg_project_point(point,
2455  camera_position,
2456  camera_direction,
2457  camera_horizontal,
2458  camera_focus);
2459 
2460  out << static_cast<unsigned int>(
2461  .5 +
2462  ((projection_decomposition[0] - x_min_perspective) /
2463  x_dimension_perspective) *
2464  (width - (width / 100.) * 2. * margin_in_percent) +
2465  ((width / 100.) * margin_in_percent))
2466  << ' '
2467  << static_cast<unsigned int>(
2468  .5 + height - (height / 100.) * margin_in_percent -
2469  ((projection_decomposition[1] - y_min_perspective) /
2470  y_dimension_perspective) *
2471  (height - (height / 100.) * 2. * margin_in_percent));
2472 
2473  out << "\"/>" << '\n';
2474 
2475  // label the current cell
2479  {
2480  point[0] = cell->center()[0];
2481  point[1] = cell->center()[1];
2482  point[2] = 0;
2483 
2485  {
2486  point[2] = svg_flags.level_height_factor *
2487  (static_cast<float>(cell->level()) /
2488  static_cast<float>(n_levels)) *
2489  std::max(x_dimension, y_dimension);
2490  }
2491 
2492  const double distance_to_camera =
2493  std::sqrt(std::pow(point[0] - camera_position[0], 2.) +
2494  std::pow(point[1] - camera_position[1], 2.) +
2495  std::pow(point[2] - camera_position[2], 2.));
2496  const double distance_factor =
2497  distance_to_camera / (2. * std::max(x_dimension, y_dimension));
2498 
2499  projection_decomposition = svg_project_point(point,
2500  camera_position,
2501  camera_direction,
2502  camera_horizontal,
2503  camera_focus);
2504 
2505  const unsigned int font_size_this_cell =
2506  static_cast<unsigned int>(
2507  .5 +
2508  cell_label_font_size *
2509  std::pow(.5, cell->level() - 4. + 3.5 * distance_factor));
2510 
2511  out << " <text"
2512  << " x=\""
2513  << static_cast<unsigned int>(
2514  .5 +
2515  ((projection_decomposition[0] - x_min_perspective) /
2516  x_dimension_perspective) *
2517  (width - (width / 100.) * 2. * margin_in_percent) +
2518  ((width / 100.) * margin_in_percent))
2519  << "\" y=\""
2520  << static_cast<unsigned int>(
2521  .5 + height - (height / 100.) * margin_in_percent -
2522  ((projection_decomposition[1] - y_min_perspective) /
2523  y_dimension_perspective) *
2524  (height - (height / 100.) * 2. * margin_in_percent) +
2525  0.5 * font_size_this_cell)
2526  << "\" style=\"font-size:" << font_size_this_cell << "px\">";
2527 
2529  {
2530  out << cell->level();
2531  }
2532 
2534  {
2536  out << '.';
2537  out << cell->index();
2538  }
2539 
2541  {
2544  out << ',';
2545  out
2546  << static_cast<std::make_signed<types::material_id>::type>(
2547  cell->material_id());
2548  }
2549 
2551  {
2554  out << ',';
2555  if (cell->is_active())
2556  out << static_cast<
2557  std::make_signed<types::subdomain_id>::type>(
2558  cell->subdomain_id());
2559  else
2560  out << 'X';
2561  }
2562 
2564  {
2569  out << ',';
2570  out
2571  << static_cast<std::make_signed<types::subdomain_id>::type>(
2572  cell->level_subdomain_id());
2573  }
2574 
2575  out << "</text>" << '\n';
2576  }
2577 
2578  // if the current cell lies at the boundary of the triangulation, draw
2579  // the additional boundary line
2581  {
2582  for (auto faceIndex : cell->face_indices())
2583  {
2584  if (cell->at_boundary(faceIndex))
2585  {
2586  point[0] = cell->face(faceIndex)->vertex(0)[0];
2587  point[1] = cell->face(faceIndex)->vertex(0)[1];
2588  point[2] = 0;
2589 
2591  {
2592  point[2] = svg_flags.level_height_factor *
2593  (static_cast<float>(cell->level()) /
2594  static_cast<float>(n_levels)) *
2595  std::max(x_dimension, y_dimension);
2596  }
2597 
2598  projection_decomposition =
2599  svg_project_point(point,
2600  camera_position,
2601  camera_direction,
2602  camera_horizontal,
2603  camera_focus);
2604 
2605  out << " <line x1=\""
2606  << static_cast<unsigned int>(
2607  .5 +
2608  ((projection_decomposition[0] -
2609  x_min_perspective) /
2610  x_dimension_perspective) *
2611  (width -
2612  (width / 100.) * 2. * margin_in_percent) +
2613  ((width / 100.) * margin_in_percent))
2614  << "\" y1=\""
2615  << static_cast<unsigned int>(
2616  .5 + height -
2617  (height / 100.) * margin_in_percent -
2618  ((projection_decomposition[1] -
2619  y_min_perspective) /
2620  y_dimension_perspective) *
2621  (height -
2622  (height / 100.) * 2. * margin_in_percent));
2623 
2624  point[0] = cell->face(faceIndex)->vertex(1)[0];
2625  point[1] = cell->face(faceIndex)->vertex(1)[1];
2626  point[2] = 0;
2627 
2629  {
2630  point[2] = svg_flags.level_height_factor *
2631  (static_cast<float>(cell->level()) /
2632  static_cast<float>(n_levels)) *
2633  std::max(x_dimension, y_dimension);
2634  }
2635 
2636  projection_decomposition =
2637  svg_project_point(point,
2638  camera_position,
2639  camera_direction,
2640  camera_horizontal,
2641  camera_focus);
2642 
2643  out << "\" x2=\""
2644  << static_cast<unsigned int>(
2645  .5 +
2646  ((projection_decomposition[0] -
2647  x_min_perspective) /
2648  x_dimension_perspective) *
2649  (width -
2650  (width / 100.) * 2. * margin_in_percent) +
2651  ((width / 100.) * margin_in_percent))
2652  << "\" y2=\""
2653  << static_cast<unsigned int>(
2654  .5 + height -
2655  (height / 100.) * margin_in_percent -
2656  ((projection_decomposition[1] -
2657  y_min_perspective) /
2658  y_dimension_perspective) *
2659  (height -
2660  (height / 100.) * 2. * margin_in_percent))
2661  << "\"/>" << '\n';
2662 
2663 
2665  {
2666  const double distance_to_camera = std::sqrt(
2667  std::pow(point[0] - camera_position[0], 2.) +
2668  std::pow(point[1] - camera_position[1], 2.) +
2669  std::pow(point[2] - camera_position[2], 2.));
2670  const double distance_factor =
2671  distance_to_camera /
2672  (2. * std::max(x_dimension, y_dimension));
2673 
2674  const unsigned int font_size_this_edge =
2675  static_cast<unsigned int>(
2676  .5 + .5 * cell_label_font_size *
2677  std::pow(.5,
2678  cell->level() - 4. +
2679  3.5 * distance_factor));
2680 
2681  point[0] = cell->face(faceIndex)->center()[0];
2682  point[1] = cell->face(faceIndex)->center()[1];
2683  point[2] = 0;
2684 
2686  {
2687  point[2] = svg_flags.level_height_factor *
2688  (static_cast<float>(cell->level()) /
2689  static_cast<float>(n_levels)) *
2690  std::max(x_dimension, y_dimension);
2691  }
2692 
2693  projection_decomposition =
2694  svg_project_point(point,
2695  camera_position,
2696  camera_direction,
2697  camera_horizontal,
2698  camera_focus);
2699 
2700  const unsigned int xc = static_cast<unsigned int>(
2701  .5 +
2702  ((projection_decomposition[0] - x_min_perspective) /
2703  x_dimension_perspective) *
2704  (width -
2705  (width / 100.) * 2. * margin_in_percent) +
2706  ((width / 100.) * margin_in_percent));
2707  const unsigned int yc = static_cast<unsigned int>(
2708  .5 + height - (height / 100.) * margin_in_percent -
2709  ((projection_decomposition[1] - y_min_perspective) /
2710  y_dimension_perspective) *
2711  (height -
2712  (height / 100.) * 2. * margin_in_percent));
2713 
2714  out << " <circle cx=\"" << xc << "\" cy=\"" << yc
2715  << "\" r=\"" << font_size_this_edge << "\" />"
2716  << '\n';
2717 
2718  out << " <text x=\"" << xc << "\" y=\"" << yc
2719  << "\" style=\"font-size:" << font_size_this_edge
2720  << "px\" dominant-baseline=\"middle\">"
2721  << static_cast<int>(
2722  cell->face(faceIndex)->boundary_id())
2723  << "</text>" << '\n';
2724  }
2725  }
2726  }
2727  }
2728  }
2729  }
2730 
2731 
2732 
2733  // draw the legend
2734  if (svg_flags.draw_legend)
2735  out << '\n' << " <!-- legend -->" << '\n';
2736 
2737  additional_width = 0;
2738  if (!svg_flags.margin)
2739  additional_width = static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2740 
2741  // explanation of the cell labeling
2742  if (svg_flags.draw_legend &&
2746  {
2747  unsigned int line_offset = 0;
2748  out << " <rect x=\"" << width + additional_width << "\" y=\""
2749  << static_cast<unsigned int>(.5 + (height / 100.) * margin_in_percent)
2750  << "\" width=\""
2751  << static_cast<unsigned int>(.5 + (height / 100.) *
2752  (40. - margin_in_percent))
2753  << "\" height=\"" << static_cast<unsigned int>(.5 + height * .215)
2754  << "\"/>" << '\n';
2755 
2756  out << " <text x=\""
2757  << width + additional_width +
2758  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2759  << "\" y=\""
2760  << static_cast<unsigned int>(.5 +
2761  (height / 100.) * margin_in_percent +
2762  (++line_offset) * 1.5 * font_size)
2763  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2764  << font_size << "px\">"
2765  << "cell label"
2766  << "</text>" << '\n';
2767 
2769  {
2770  out << " <text x=\""
2771  << width + additional_width +
2772  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2773  << "\" y=\""
2774  << static_cast<unsigned int>(.5 +
2775  (height / 100.) * margin_in_percent +
2776  (++line_offset) * 1.5 * font_size)
2777  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2778  << font_size << "px\">"
2779  << "cell_level";
2780 
2784  out << '.';
2785 
2786  out << "</text>" << '\n';
2787  }
2788 
2790  {
2791  out << " <text x=\""
2792  << width + additional_width +
2793  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2794  << "\" y=\""
2795  << static_cast<unsigned int>(.5 +
2796  (height / 100.) * margin_in_percent +
2797  (++line_offset) * 1.5 * font_size)
2798  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2799  << font_size << "px\">"
2800  << "cell_index";
2801 
2804  out << ',';
2805 
2806  out << "</text>" << '\n';
2807  }
2808 
2810  {
2811  out << " <text x=\""
2812  << width + additional_width +
2813  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2814  << "\" y=\""
2815  << static_cast<unsigned int>(.5 +
2816  (height / 100.) * margin_in_percent +
2817  (++line_offset) * 1.5 * font_size)
2818  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2819  << font_size << "px\">"
2820  << "material_id";
2821 
2824  out << ',';
2825 
2826  out << "</text>" << '\n';
2827  }
2828 
2830  {
2831  out << " <text x= \""
2832  << width + additional_width +
2833  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2834  << "\" y=\""
2835  << static_cast<unsigned int>(.5 +
2836  (height / 100.) * margin_in_percent +
2837  (++line_offset) * 1.5 * font_size)
2838  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2839  << font_size << "px\">"
2840  << "subdomain_id";
2841 
2843  out << ',';
2844 
2845  out << "</text>" << '\n';
2846  }
2847 
2849  {
2850  out << " <text x= \""
2851  << width + additional_width +
2852  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2853  << "\" y=\""
2854  << static_cast<unsigned int>(.5 +
2855  (height / 100.) * margin_in_percent +
2856  (++line_offset) * 1.5 * font_size)
2857  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2858  << font_size << "px\">"
2859  << "level_subdomain_id"
2860  << "</text>" << '\n';
2861  }
2862 
2864  {
2865  out << " <text x=\""
2866  << width + additional_width +
2867  static_cast<unsigned int>(.5 + (height / 100.) * 1.25)
2868  << "\" y=\""
2869  << static_cast<unsigned int>(.5 +
2870  (height / 100.) * margin_in_percent +
2871  (++line_offset) * 1.5 * font_size)
2872  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2873  << font_size << "px\">"
2874  << "edge label"
2875  << "</text>" << '\n';
2876 
2877  out << " <text x= \""
2878  << width + additional_width +
2879  static_cast<unsigned int>(.5 + (height / 100.) * 2.)
2880  << "\" y=\""
2881  << static_cast<unsigned int>(.5 +
2882  (height / 100.) * margin_in_percent +
2883  (++line_offset) * 1.5 * font_size)
2884  << "\" style=\"text-anchor:start; font-style:oblique; font-size:"
2885  << font_size << "px\">"
2886  << "boundary_id"
2887  << "</text>" << '\n';
2888  }
2889  }
2890 
2891  // show azimuth angle and polar angle as text below the explanation of the
2892  // cell labeling
2893  if (svg_flags.draw_legend)
2894  {
2895  out << " <text x=\"" << width + additional_width << "\" y=\""
2896  << static_cast<unsigned int>(
2897  .5 + (height / 100.) * margin_in_percent + 13.75 * font_size)
2898  << "\" style=\"text-anchor:start; font-size:" << font_size << "px\">"
2899  << "azimuth: " << svg_flags.azimuth_angle
2900  << "°, polar: " << svg_flags.polar_angle << "°</text>" << '\n';
2901  }
2902 
2903 
2904  // draw the colorbar
2906  {
2907  out << '\n' << " <!-- colorbar -->" << '\n';
2908 
2909  out << " <text x=\"" << width + additional_width << "\" y=\""
2910  << static_cast<unsigned int>(
2911  .5 + (height / 100.) * (margin_in_percent + 29.) -
2912  (font_size / 1.25))
2913  << "\" style=\"text-anchor:start; font-weight:bold; font-size:"
2914  << font_size << "px\">";
2915 
2916  switch (svg_flags.coloring)
2917  {
2918  case 1:
2919  out << "material_id";
2920  break;
2921  case 2:
2922  out << "level_number";
2923  break;
2924  case 3:
2925  out << "subdomain_id";
2926  break;
2927  case 4:
2928  out << "level_subdomain_id";
2929  break;
2930  default:
2931  break;
2932  }
2933 
2934  out << "</text>" << '\n';
2935 
2936  unsigned int element_height = static_cast<unsigned int>(
2937  ((height / 100.) * (71. - 2. * margin_in_percent)) / n);
2938  unsigned int element_width =
2939  static_cast<unsigned int>(.5 + (height / 100.) * 2.5);
2940 
2941  int labeling_index = 0;
2942  auto materials_it = materials.begin();
2943  auto levels_it = levels.begin();
2944  auto subdomains_it = subdomains.begin();
2945  auto level_subdomains_it = level_subdomains.begin();
2946 
2947  for (unsigned int index = 0; index < n; ++index)
2948  {
2949  switch (svg_flags.coloring)
2950  {
2952  labeling_index = *materials_it++;
2953  break;
2955  labeling_index = *levels_it++;
2956  break;
2958  labeling_index = *subdomains_it++;
2959  break;
2961  labeling_index = *level_subdomains_it++;
2962  break;
2963  default:
2964  break;
2965  }
2966 
2967  out << " <rect class=\"r" << labeling_index << "\" x=\""
2968  << width + additional_width << "\" y=\""
2969  << static_cast<unsigned int>(.5 + (height / 100.) *
2970  (margin_in_percent + 29)) +
2971  (n - index - 1) * element_height
2972  << "\" width=\"" << element_width << "\" height=\""
2973  << element_height << "\"/>" << '\n';
2974 
2975  out << " <text x=\""
2976  << width + additional_width + 1.5 * element_width << "\" y=\""
2977  << static_cast<unsigned int>(.5 + (height / 100.) *
2978  (margin_in_percent + 29)) +
2979  (n - index - 1 + .5) * element_height +
2980  static_cast<unsigned int>(.5 + font_size * .35)
2981  << "\""
2982  << " style=\"text-anchor:start; font-size:"
2983  << static_cast<unsigned int>(.5 + font_size) << "px";
2984 
2985  if (index == 0 || index == n - 1)
2986  out << "; font-weight:bold";
2987 
2988  out << "\">" << labeling_index;
2989 
2990  if (index == n - 1)
2991  out << " max";
2992  if (index == 0)
2993  out << " min";
2994 
2995  out << "</text>" << '\n';
2996 
2997  labeling_index++;
2998  }
2999  }
3000 
3001 
3002  // finalize the svg file
3003  out << '\n' << "</svg>";
3004  out.flush();
3005 }
3006 
3007 
3008 
3009 template <>
3010 void
3011 GridOut::write_mathgl(const Triangulation<1> &, std::ostream &) const
3012 {
3013  // 1d specialization not done yet
3014  Assert(false, ExcNotImplemented());
3015 }
3016 
3017 
3018 
3019 template <int dim, int spacedim>
3020 void
3022  std::ostream & out) const
3023 {
3024  AssertThrow(out, ExcIO());
3025 
3026  // (i) write header
3027  {
3028  // block this to have local variables destroyed after use
3029  const std::time_t time1 = std::time(nullptr);
3030  const std::tm * time = std::localtime(&time1);
3031 
3032  out
3033  << "\n#"
3034  << "\n# This file was generated by the deal.II library."
3035  << "\n# Date = " << time->tm_year + 1900 << "/" << std::setfill('0')
3036  << std::setw(2) << time->tm_mon + 1 << "/" << std::setfill('0')
3037  << std::setw(2) << time->tm_mday << "\n# Time = " << std::setfill('0')
3038  << std::setw(2) << time->tm_hour << ":" << std::setfill('0')
3039  << std::setw(2) << time->tm_min << ":" << std::setfill('0')
3040  << std::setw(2) << time->tm_sec << "\n#"
3041  << "\n# For a description of the MathGL script format see the MathGL manual. "
3042  << "\n#"
3043  << "\n# Note: This file is understood by MathGL v2.1 and higher only, and can "
3044  << "\n# be quickly viewed in a graphical environment using \'mglview\'. "
3045  << "\n#"
3046  << "\n";
3047  }
3048 
3049  // define a helper to keep loops approximately dim-independent
3050  // since MathGL labels axes as x, y, z
3051  const std::string axes = "xyz";
3052 
3053  // (ii) write preamble and graphing tweaks
3054  out << "\n#"
3055  << "\n# Preamble."
3056  << "\n#"
3057  << "\n";
3058 
3060  out << "\nbox";
3061 
3062  // deal with dimension dependent preamble; eg. default sizes and
3063  // views for MathGL (cf. gnuplot).
3064  switch (dim)
3065  {
3066  case 2:
3067  out << "\nsetsize 800 800";
3068  out << "\nrotate 0 0";
3069  break;
3070  case 3:
3071  out << "\nsetsize 800 800";
3072  out << "\nrotate 60 40";
3073  break;
3074  default:
3075  Assert(false, ExcNotImplemented());
3076  }
3077  out << "\n";
3078 
3079  // (iii) write vertex ordering
3080  out << "\n#"
3081  << "\n# Vertex ordering."
3082  << "\n# list <vertex order> <vertex indices>"
3083  << "\n#"
3084  << "\n";
3085 
3086  // todo: This denotes the natural ordering of vertices, but it needs
3087  // to check this is really always true for a given grid (it's not
3088  // true in @ref step_1 "step-1" grid-2 for instance).
3089  switch (dim)
3090  {
3091  case 2:
3092  out << "\nlist f 0 1 2 3"
3093  << "\n";
3094  break;
3095  case 3:
3096  out
3097  << "\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"
3098  << "\n";
3099  break;
3100  default:
3101  Assert(false, ExcNotImplemented());
3102  }
3103 
3104  // (iv) write a list of vertices of cells
3105  out << "\n#"
3106  << "\n# List of vertices."
3107  << "\n# list <id> <vertices>"
3108  << "\n#"
3109  << "\n";
3110 
3111  // run over all active cells and write out a list of
3112  // xyz-coordinates that correspond to vertices
3113  // No global indices in deal.II, so we make one up here.
3114  for (const auto &cell : tria.active_cell_iterators())
3115  {
3116  for (unsigned int i = 0; i < dim; ++i)
3117  {
3118  // if (cell->direction_flag ()==true)
3119  // out << "\ntrue";
3120  // else
3121  // out << "\nfalse";
3122 
3123  out << "\nlist " << axes[i] << cell->active_cell_index() << " ";
3124  for (const unsigned int j : GeometryInfo<dim>::vertex_indices())
3125  out << cell->vertex(j)[i] << " ";
3126  }
3127  out << '\n';
3128  }
3129 
3130  // (v) write out cells to plot as quadplot objects
3131  out << "\n#"
3132  << "\n# List of cells to quadplot."
3133  << "\n# quadplot <vertex order> <id> <style>"
3134  << "\n#"
3135  << "\n";
3136  for (unsigned int i = 0; i < tria.n_active_cells(); ++i)
3137  {
3138  out << "\nquadplot f ";
3139  for (unsigned int j = 0; j < dim; ++j)
3140  out << axes[j] << i << " ";
3141  out << "\'k#\'";
3142  }
3143  out << "\n";
3144 
3145  // (vi) write footer
3146  out << "\n#"
3147  << "\n#"
3148  << "\n#"
3149  << "\n";
3150 
3151  // make sure everything now gets to the output stream
3152  out.flush();
3153  AssertThrow(out, ExcIO());
3154 }
3155 
3156 
3157 
3158 namespace
3159 {
3166  template <int dim, int spacedim, typename ITERATOR, typename END>
3167  void
3168  generate_triangulation_patches(
3169  std::vector<DataOutBase::Patch<dim, spacedim>> &patches,
3170  ITERATOR cell,
3171  END end)
3172  {
3173  // convert each of the active cells into a patch
3174  for (; cell != end; ++cell)
3175  {
3177  patch.reference_cell = cell->reference_cell();
3178  patch.n_subdivisions = 1;
3179  patch.data.reinit(5, cell->n_vertices());
3180 
3181  for (const unsigned int v : cell->vertex_indices())
3182  {
3183  patch.vertices[v] = cell->vertex(v);
3184  patch.data(0, v) = cell->level();
3185  patch.data(1, v) =
3186  static_cast<std::make_signed<types::manifold_id>::type>(
3187  cell->manifold_id());
3188  patch.data(2, v) =
3189  static_cast<std::make_signed<types::material_id>::type>(
3190  cell->material_id());
3191  if (cell->is_active())
3192  patch.data(3, v) =
3193  static_cast<std::make_signed<types::subdomain_id>::type>(
3194  cell->subdomain_id());
3195  else
3196  patch.data(3, v) = -1;
3197  patch.data(4, v) =
3198  static_cast<std::make_signed<types::subdomain_id>::type>(
3199  cell->level_subdomain_id());
3200  }
3201  patches.push_back(patch);
3202  }
3203  }
3204 
3205 
3206 
3207  std::vector<std::string>
3208  triangulation_patch_data_names()
3209  {
3210  std::vector<std::string> v(5);
3211  v[0] = "level";
3212  v[1] = "manifold";
3213  v[2] = "material";
3214  v[3] = "subdomain";
3215  v[4] = "level_subdomain";
3216  return v;
3217  }
3218 
3222  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3223  get_boundary_edge_iterators(const Triangulation<3, 3> &tria)
3224  {
3225  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3226 
3227  std::vector<bool> flags;
3228  tria.save_user_flags_line(flags);
3229  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3230 
3231  for (auto face : tria.active_face_iterators())
3232  for (const auto l : face->line_indices())
3233  {
3234  const auto line = face->line(l);
3235  if (line->user_flag_set() || line->has_children())
3236  continue;
3237  else
3238  line->set_user_flag();
3239  if (line->at_boundary())
3240  res.emplace_back(line);
3241  }
3242  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3243  return res;
3244  }
3245 
3246 
3247 
3251  template <int dim, int spacedim>
3252  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3253  get_boundary_edge_iterators(const Triangulation<dim, spacedim> &)
3254  {
3255  return {};
3256  }
3257 
3258 
3259 
3264  std::vector<typename Triangulation<3, 3>::active_line_iterator>
3265  get_relevant_edge_iterators(const Triangulation<3, 3> &tria)
3266  {
3267  std::vector<typename Triangulation<3, 3>::active_line_iterator> res;
3268 
3269  std::vector<bool> flags;
3270  tria.save_user_flags_line(flags);
3271  const_cast<Triangulation<3, 3> &>(tria).clear_user_flags_line();
3272 
3273  for (auto face : tria.active_face_iterators())
3274  for (const auto l : face->line_indices())
3275  {
3276  const auto line = face->line(l);
3277  if (line->user_flag_set() || line->has_children())
3278  continue;
3279  else
3280  line->set_user_flag();
3281  if (line->manifold_id() != numbers::flat_manifold_id ||
3282  (line->boundary_id() != 0 &&
3283  line->boundary_id() != numbers::invalid_boundary_id))
3284  res.emplace_back(line);
3285  }
3286  const_cast<Triangulation<3, 3> &>(tria).load_user_flags_line(flags);
3287  return res;
3288  }
3289 
3290 
3294  template <int dim, int spacedim>
3295  std::vector<typename Triangulation<dim, spacedim>::active_line_iterator>
3296  get_relevant_edge_iterators(const Triangulation<dim, spacedim> &)
3297  {
3298  return {};
3299  }
3300 
3301 
3302 
3306  template <int dim, int spacedim>
3307  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3308  get_boundary_face_iterators(const Triangulation<dim, spacedim> &tria)
3309  {
3310  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3311  res;
3312  if (dim == 1)
3313  return res;
3314  for (auto face : tria.active_face_iterators())
3315  {
3316  if (face->boundary_id() != numbers::invalid_boundary_id)
3317  res.push_back(face);
3318  }
3319  return res;
3320  }
3321 
3322 
3323 
3328  template <int dim, int spacedim>
3329  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3330  get_relevant_face_iterators(const Triangulation<dim, spacedim> &tria)
3331  {
3332  std::vector<typename Triangulation<dim, spacedim>::active_face_iterator>
3333  res;
3334  if (dim == 1)
3335  return res;
3336  for (auto face : tria.active_face_iterators())
3337  {
3338  if (face->manifold_id() != numbers::flat_manifold_id ||
3339  (face->boundary_id() != 0 &&
3340  face->boundary_id() != numbers::invalid_boundary_id))
3341  res.push_back(face);
3342  }
3343  return res;
3344  }
3345 } // namespace
3346 
3347 
3348 
3349 template <int dim, int spacedim>
3350 void
3352  std::ostream & out) const
3353 {
3354  AssertThrow(out, ExcIO());
3355 
3356  // get the positions of the vertices
3357  const std::vector<Point<spacedim>> &vertices = tria.get_vertices();
3358 
3359  const auto n_vertices = vertices.size();
3360 
3361  out << "# vtk DataFile Version 3.0\n"
3362  << "Triangulation generated with deal.II\n"
3363  << "ASCII\n"
3364  << "DATASET UNSTRUCTURED_GRID\n"
3365  << "POINTS " << n_vertices << " double\n";
3366 
3367  // actually write the vertices.
3368  for (const auto &v : vertices)
3369  {
3370  out << v;
3371  for (unsigned int d = spacedim + 1; d <= 3; ++d)
3372  out << " 0"; // fill with zeroes
3373  out << '\n';
3374  }
3375 
3376  const auto faces = vtk_flags.output_only_relevant ?
3377  get_relevant_face_iterators(tria) :
3378  get_boundary_face_iterators(tria);
3379  const auto edges = vtk_flags.output_only_relevant ?
3380  get_relevant_edge_iterators(tria) :
3381  get_boundary_edge_iterators(tria);
3382 
3383  AssertThrow(
3384  vtk_flags.output_cells || (dim >= 2 && vtk_flags.output_faces) ||
3385  (dim >= 3 && vtk_flags.output_edges),
3386  ExcMessage(
3387  "At least one of the flags (output_cells, output_faces, output_edges) has to be enabled!"));
3388 
3389  // Write cells preamble
3390  const int n_cells = (vtk_flags.output_cells ? tria.n_active_cells() : 0) +
3391  (vtk_flags.output_faces ? faces.size() : 0) +
3392  (vtk_flags.output_edges ? edges.size() : 0);
3393 
3394  // VTK now expects a number telling the total storage requirement to read all
3395  // cell connectivity information. The connectivity information is read cell by
3396  // cell, first specifying how many vertices are required to describe the cell,
3397  // and then specifying the index of every vertex. This means that for every
3398  // deal.II object type, we always need n_vertices + 1 integer per cell.
3399  // Compute the total number here.
3400  int cells_size = 0;
3401 
3402  if (vtk_flags.output_cells)
3403  for (const auto &cell : tria.active_cell_iterators())
3404  cells_size += cell->n_vertices() + 1;
3405 
3406  if (vtk_flags.output_faces)
3407  for (const auto &face : faces)
3408  cells_size += face->n_vertices() + 1;
3409 
3410  if (vtk_flags.output_edges)
3411  for (const auto &edge : edges)
3412  cells_size += edge->n_vertices() + 1;
3413 
3414  AssertThrow(cells_size > 0, ExcMessage("No cells given to be output!"));
3415 
3416  out << "\nCELLS " << n_cells << ' ' << cells_size << '\n';
3417  /*
3418  * VTK cells:
3419  *
3420  * 1 VTK_VERTEX
3421  * 3 VTK_LINE
3422  * 5 VTK_TRIANGLE
3423  * 9 VTK_QUAD
3424  * 10 VTK_TETRA
3425  * 14 VTK_PYRAMID
3426  * 13 VTK_WEDGE
3427  * 12 VTK_HEXAHEDRON
3428  *
3429  * see also: https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
3430  */
3431  static const std::array<int, 8> deal_to_vtk_cell_type = {
3432  {1, 3, 5, 9, 10, 14, 13, 12}};
3433  static const std::array<unsigned int, 8> vtk_to_deal_hypercube = {
3434  {0, 1, 3, 2, 4, 5, 7, 6}};
3435 
3436  // write cells.
3437  if (vtk_flags.output_cells)
3438  for (const auto &cell : tria.active_cell_iterators())
3439  {
3440  out << cell->n_vertices();
3441  for (const unsigned int i : cell->vertex_indices())
3442  {
3443  out << ' ';
3444  const auto reference_cell = cell->reference_cell();
3445 
3450  out << cell->vertex_index(vtk_to_deal_hypercube[i]);
3451  else if ((reference_cell == ReferenceCells::Triangle) ||
3454  out << cell->vertex_index(i);
3456  {
3457  static const std::array<unsigned int, 5> permutation_table{
3458  {0, 1, 3, 2, 4}};
3459  out << cell->vertex_index(permutation_table[i]);
3460  }
3461  else
3462  Assert(false, ExcNotImplemented());
3463  }
3464  out << '\n';
3465  }
3466  if (vtk_flags.output_faces)
3467  for (const auto &face : faces)
3468  {
3469  out << face->n_vertices();
3470  for (const unsigned int i : face->vertex_indices())
3471  {
3472  out << ' '
3473  << face->vertex_index(GeometryInfo<dim>::vertices_per_face ==
3474  face->n_vertices() ?
3475  vtk_to_deal_hypercube[i] :
3476  i);
3477  }
3478  out << '\n';
3479  }
3480  if (vtk_flags.output_edges)
3481  for (const auto &edge : edges)
3482  {
3483  out << 2;
3484  for (const unsigned int i : edge->vertex_indices())
3485  out << ' ' << edge->vertex_index(i);
3486  out << '\n';
3487  }
3488 
3489  // write cell types
3490  out << "\nCELL_TYPES " << n_cells << '\n';
3491  if (vtk_flags.output_cells)
3492  {
3493  for (const auto &cell : tria.active_cell_iterators())
3494  out << deal_to_vtk_cell_type[static_cast<int>(cell->reference_cell())]
3495  << ' ';
3496  out << '\n';
3497  }
3498  if (vtk_flags.output_faces)
3499  {
3500  for (const auto &face : faces)
3501  out << deal_to_vtk_cell_type[static_cast<int>(face->reference_cell())]
3502  << ' ';
3503  out << '\n';
3504  }
3505  if (vtk_flags.output_edges)
3506  {
3507  for (const auto &edge : edges)
3508  out << deal_to_vtk_cell_type[static_cast<int>(edge->reference_cell())]
3509  << ' ';
3510  }
3511  out << "\n\nCELL_DATA " << n_cells << '\n'
3512  << "SCALARS MaterialID int 1\n"
3513  << "LOOKUP_TABLE default\n";
3514 
3515  // Now material id and boundary id
3516  if (vtk_flags.output_cells)
3517  {
3518  for (const auto &cell : tria.active_cell_iterators())
3519  {
3520  out << static_cast<std::make_signed<types::material_id>::type>(
3521  cell->material_id())
3522  << ' ';
3523  }
3524  out << '\n';
3525  }
3526  if (vtk_flags.output_faces)
3527  {
3528  for (const auto &face : faces)
3529  {
3530  out << static_cast<std::make_signed<types::boundary_id>::type>(
3531  face->boundary_id())
3532  << ' ';
3533  }
3534  out << '\n';
3535  }
3536  if (vtk_flags.output_edges)
3537  {
3538  for (const auto &edge : edges)
3539  {
3540  out << static_cast<std::make_signed<types::boundary_id>::type>(
3541  edge->boundary_id())
3542  << ' ';
3543  }
3544  }
3545 
3546  out << "\n\nSCALARS ManifoldID int 1\n"
3547  << "LOOKUP_TABLE default\n";
3548 
3549  // Now manifold id
3550  if (vtk_flags.output_cells)
3551  {
3552  for (const auto &cell : tria.active_cell_iterators())
3553  {
3554  out << static_cast<std::make_signed<types::manifold_id>::type>(
3555  cell->manifold_id())
3556  << ' ';
3557  }
3558  out << '\n';
3559  }
3560  if (vtk_flags.output_faces)
3561  {
3562  for (const auto &face : faces)
3563  {
3564  out << static_cast<std::make_signed<types::manifold_id>::type>(
3565  face->manifold_id())
3566  << ' ';
3567  }
3568  out << '\n';
3569  }
3570  if (vtk_flags.output_edges)
3571  {
3572  for (const auto &edge : edges)
3573  {
3574  out << static_cast<std::make_signed<types::manifold_id>::type>(
3575  edge->manifold_id())
3576  << ' ';
3577  }
3578  out << '\n';
3579  }
3580 
3581  out.flush();
3582 
3583  AssertThrow(out, ExcIO());
3584 }
3585 
3586 
3587 
3588 template <int dim, int spacedim>
3589 void
3591  std::ostream & out) const
3592 {
3593  AssertThrow(out, ExcIO());
3594 
3595  // convert the cells of the triangulation into a set of patches
3596  // and then have them output. since there is no data attached to
3597  // the geometry, we also do not have to provide any names, identifying
3598  // information, etc.
3599  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3600  patches.reserve(tria.n_active_cells());
3601  generate_triangulation_patches(patches, tria.begin_active(), tria.end());
3602 
3605  patches,
3606  triangulation_patch_data_names(),
3607  std::vector<
3608  std::tuple<unsigned int,
3609  unsigned int,
3610  std::string,
3612  vtu_flags,
3613  out);
3615  {
3616  out << " </UnstructuredGrid>\n";
3617  out << "<dealiiData encoding=\"base64\">";
3618  std::stringstream outstring;
3619  boost::archive::binary_oarchive ia(outstring);
3620  tria.save(ia, 0);
3621  const auto compressed = Utilities::compress(outstring.str());
3622  out << Utilities::encode_base64({compressed.begin(), compressed.end()});
3623  out << "\n</dealiiData>\n";
3624  out << "</VTKFile>\n";
3625  }
3626  else
3628 
3629  out << std::flush;
3630  AssertThrow(out, ExcIO());
3631 }
3632 
3633 
3634 
3635 template <int dim, int spacedim>
3636 void
3639  const std::string & filename_without_extension,
3640  const bool view_levels,
3641  const bool include_artificial) const
3642 {
3643  std::vector<DataOutBase::Patch<dim, spacedim>> patches;
3644  const unsigned int n_datasets = 4;
3645  std::vector<std::string> data_names;
3646  data_names.emplace_back("level");
3647  data_names.emplace_back("subdomain");
3648  data_names.emplace_back("level_subdomain");
3649  data_names.emplace_back("proc_writing");
3650 
3651  const auto reference_cells = tria.get_reference_cells();
3652 
3653  AssertDimension(reference_cells.size(), 1);
3654 
3655  const auto &reference_cell = reference_cells[0];
3656 
3657  const unsigned int n_q_points = reference_cell.n_vertices();
3658 
3659  for (const auto &cell : tria.cell_iterators())
3660  {
3661  if (!view_levels)
3662  {
3663  if (cell->has_children())
3664  continue;
3665  if (!include_artificial &&
3666  cell->subdomain_id() == numbers::artificial_subdomain_id)
3667  continue;
3668  }
3669  else if (!include_artificial)
3670  {
3671  if (cell->has_children() &&
3672  cell->level_subdomain_id() == numbers::artificial_subdomain_id)
3673  continue;
3674  else if (cell->is_active() &&
3675  cell->level_subdomain_id() ==
3677  cell->subdomain_id() == numbers::artificial_subdomain_id)
3678  continue;
3679  }
3680 
3682  patch.data.reinit(n_datasets, n_q_points);
3683  patch.points_are_available = false;
3685 
3686  for (unsigned int vertex = 0; vertex < n_q_points; ++vertex)
3687  {
3688  patch.vertices[vertex] = cell->vertex(vertex);
3689  patch.data(0, vertex) = cell->level();
3690  if (cell->is_active())
3691  patch.data(1, vertex) = static_cast<double>(
3692  static_cast<std::make_signed<types::subdomain_id>::type>(
3693  cell->subdomain_id()));
3694  else
3695  patch.data(1, vertex) = -1.0;
3696  patch.data(2, vertex) = static_cast<double>(
3697  static_cast<std::make_signed<types::subdomain_id>::type>(
3698  cell->level_subdomain_id()));
3699  patch.data(3, vertex) = tria.locally_owned_subdomain();
3700  }
3701 
3702  for (auto f : reference_cell.face_indices())
3704  patches.push_back(patch);
3705  }
3706 
3707  // only create .pvtu file if running in parallel
3708  // if not, just create a .vtu file with no reference
3709  // to the processor number
3710  std::string new_file = filename_without_extension + ".vtu";
3712  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
3713  {
3714  new_file = filename_without_extension + ".proc" +
3715  Utilities::int_to_string(tr->locally_owned_subdomain(), 4) +
3716  ".vtu";
3717 
3718  // create .pvtu record
3719  if (tr->locally_owned_subdomain() == 0)
3720  {
3721  std::vector<std::string> filenames;
3722 
3723  // .pvtu needs to reference the files without a relative path because
3724  // it will be written in the same directory. For this, remove any
3725  // paths from filename.
3726  std::size_t pos = filename_without_extension.find_last_of('/');
3727  if (pos == std::string::npos)
3728  pos = 0;
3729  else
3730  pos += 1;
3731  const unsigned int n_procs =
3732  Utilities::MPI::n_mpi_processes(tr->get_communicator());
3733  for (unsigned int i = 0; i < n_procs; ++i)
3734  filenames.push_back(filename_without_extension.substr(pos) +
3735  ".proc" + Utilities::int_to_string(i, 4) +
3736  ".vtu");
3737 
3738  const std::string pvtu_filename =
3739  (filename_without_extension + ".pvtu");
3740  std::ofstream pvtu_output(pvtu_filename.c_str());
3741 
3742  DataOut<dim, spacedim> data_out;
3743  data_out.attach_triangulation(*tr);
3744 
3745  // We need a dummy vector with the names of the data values in the
3746  // .vtu files in order that the .pvtu contains reference these values
3747  Vector<float> dummy_vector(tr->n_active_cells());
3748  data_out.add_data_vector(dummy_vector, "level");
3749  data_out.add_data_vector(dummy_vector, "subdomain");
3750  data_out.add_data_vector(dummy_vector, "level_subdomain");
3751  data_out.add_data_vector(dummy_vector, "proc_writing");
3752 
3753  data_out.build_patches();
3754 
3755  data_out.write_pvtu_record(pvtu_output, filenames);
3756  }
3757  }
3758 
3759  std::ofstream out(new_file.c_str());
3760  std::vector<
3761  std::tuple<unsigned int,
3762  unsigned int,
3763  std::string,
3765  vector_data_ranges;
3766  DataOutBase::VtkFlags flags;
3767  DataOutBase::write_vtu(patches, data_names, vector_data_ranges, flags, out);
3768 }
3769 
3770 
3771 
3772 unsigned int
3774 {
3775  return 0;
3776 }
3777 
3778 unsigned int
3780 {
3781  return 0;
3782 }
3783 
3784 
3785 unsigned int
3787 {
3788  return 0;
3789 }
3790 
3791 unsigned int
3793 {
3794  return 0;
3795 }
3796 
3797 unsigned int
3799 {
3800  return 0;
3801 }
3802 
3803 unsigned int
3805 {
3806  return 0;
3807 }
3808 
3809 unsigned int
3811 {
3812  return 0;
3813 }
3814 
3815 unsigned int
3817 {
3818  return 0;
3819 }
3820 
3821 
3822 
3823 template <int dim, int spacedim>
3824 unsigned int
3826 {
3828  unsigned int n_faces = 0;
3829 
3830  for (const auto &face : tria.active_face_iterators())
3831  if ((face->at_boundary()) && (face->boundary_id() != 0))
3832  n_faces++;
3833 
3834  return n_faces;
3835 }
3836 
3837 
3838 
3839 template <int dim, int spacedim>
3840 unsigned int
3842 {
3843  // save the user flags for lines so
3844  // we can use these flags to track
3845  // which ones we've already counted
3846  std::vector<bool> line_flags;
3847  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3848  line_flags);
3849  const_cast<::Triangulation<dim, spacedim> &>(tria)
3850  .clear_user_flags_line();
3851 
3852  unsigned int n_lines = 0;
3853 
3854  for (const auto &cell : tria.active_cell_iterators())
3855  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
3856  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
3857  (cell->line(l)->user_flag_set() == false))
3858  {
3859  ++n_lines;
3860  cell->line(l)->set_user_flag();
3861  }
3862 
3863  // at the end, restore the user
3864  // flags for the lines
3865  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
3866  line_flags);
3867 
3868  return n_lines;
3869 }
3870 
3871 
3872 
3873 unsigned int
3875  const unsigned int next_element_index,
3876  std::ostream &) const
3877 {
3878  return next_element_index;
3879 }
3880 
3881 
3882 unsigned int
3884  const unsigned int next_element_index,
3885  std::ostream &) const
3886 {
3887  return next_element_index;
3888 }
3889 
3890 unsigned int
3892  const unsigned int next_element_index,
3893  std::ostream &) const
3894 {
3895  return next_element_index;
3896 }
3897 
3898 
3899 unsigned int
3901  const unsigned int next_element_index,
3902  std::ostream &) const
3903 {
3904  return next_element_index;
3905 }
3906 
3907 unsigned int
3909  const unsigned int next_element_index,
3910  std::ostream &) const
3911 {
3912  return next_element_index;
3913 }
3914 
3915 
3916 unsigned int
3918  const unsigned int next_element_index,
3919  std::ostream &) const
3920 {
3921  return next_element_index;
3922 }
3923 
3924 
3925 unsigned int
3927  const unsigned int next_element_index,
3928  std::ostream &) const
3929 {
3930  return next_element_index;
3931 }
3932 
3933 unsigned int
3935  const unsigned int next_element_index,
3936  std::ostream &) const
3937 {
3938  return next_element_index;
3939 }
3940 
3941 
3942 
3943 template <int dim, int spacedim>
3944 unsigned int
3946  const unsigned int next_element_index,
3947  std::ostream & out) const
3948 {
3949  unsigned int current_element_index = next_element_index;
3950 
3951  for (const auto &face : tria.active_face_iterators())
3952  if (face->at_boundary() && (face->boundary_id() != 0))
3953  {
3954  out << current_element_index << ' ';
3955  switch (dim)
3956  {
3957  case 2:
3958  out << 1 << ' ';
3959  break;
3960  case 3:
3961  out << 3 << ' ';
3962  break;
3963  default:
3964  Assert(false, ExcNotImplemented());
3965  }
3966  out << static_cast<unsigned int>(face->boundary_id()) << ' '
3967  << static_cast<unsigned int>(face->boundary_id()) << ' '
3969  // note: vertex numbers are 1-base
3970  for (unsigned int vertex = 0;
3971  vertex < GeometryInfo<dim>::vertices_per_face;
3972  ++vertex)
3973  out << ' '
3974  << face->vertex_index(
3976  1;
3977  out << '\n';
3978 
3979  ++current_element_index;
3980  }
3981  return current_element_index;
3982 }
3983 
3984 
3985 template <int dim, int spacedim>
3986 unsigned int
3988  const unsigned int next_element_index,
3989  std::ostream & out) const
3990 {
3991  unsigned int current_element_index = next_element_index;
3992  // save the user flags for lines so
3993  // we can use these flags to track
3994  // which ones we've already taken
3995  // care of
3996  std::vector<bool> line_flags;
3997  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
3998  line_flags);
3999  const_cast<::Triangulation<dim, spacedim> &>(tria)
4000  .clear_user_flags_line();
4001 
4002  for (const auto &cell : tria.active_cell_iterators())
4003  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4004  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
4005  (cell->line(l)->user_flag_set() == false))
4006  {
4007  out << next_element_index << " 1 ";
4008  out << static_cast<unsigned int>(cell->line(l)->boundary_id()) << ' '
4009  << static_cast<unsigned int>(cell->line(l)->boundary_id())
4010  << " 2 ";
4011  // note: vertex numbers are 1-base
4012  for (unsigned int vertex = 0; vertex < 2; ++vertex)
4013  out << ' '
4014  << cell->line(l)->vertex_index(
4016  1;
4017  out << '\n';
4018 
4019  // move on to the next line
4020  // but mark the current one
4021  // as taken care of
4022  ++current_element_index;
4023  cell->line(l)->set_user_flag();
4024  }
4025 
4026  // at the end, restore the user
4027  // flags for the lines
4028  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
4029  line_flags);
4030 
4031  return current_element_index;
4032 }
4033 
4034 
4035 
4036 unsigned int
4038  const unsigned int next_element_index,
4039  std::ostream &) const
4040 {
4041  return next_element_index;
4042 }
4043 
4044 unsigned int
4046  const unsigned int next_element_index,
4047  std::ostream &) const
4048 {
4049  return next_element_index;
4050 }
4051 
4052 unsigned int
4054  const unsigned int next_element_index,
4055  std::ostream &) const
4056 {
4057  return next_element_index;
4058 }
4059 
4060 unsigned int
4062  const unsigned int next_element_index,
4063  std::ostream &) const
4064 {
4065  return next_element_index;
4066 }
4067 
4068 unsigned int
4070  const unsigned int next_element_index,
4071  std::ostream &) const
4072 {
4073  return next_element_index;
4074 }
4075 
4076 
4077 unsigned int
4079  const unsigned int next_element_index,
4080  std::ostream &) const
4081 {
4082  return next_element_index;
4083 }
4084 
4085 
4086 unsigned int
4088  const unsigned int next_element_index,
4089  std::ostream &) const
4090 {
4091  return next_element_index;
4092 }
4093 
4094 unsigned int
4096  const unsigned int next_element_index,
4097  std::ostream &) const
4098 {
4099  return next_element_index;
4100 }
4101 
4102 
4103 
4104 template <int dim, int spacedim>
4105 unsigned int
4107  const unsigned int next_element_index,
4108  std::ostream & out) const
4109 {
4110  unsigned int current_element_index = next_element_index;
4112 
4113  for (const auto &face : tria.active_face_iterators())
4114  if (face->at_boundary() && (face->boundary_id() != 0))
4115  {
4116  out << current_element_index << " "
4117  << static_cast<unsigned int>(face->boundary_id()) << " ";
4118  switch (dim)
4119  {
4120  case 2:
4121  out << "line ";
4122  break;
4123  case 3:
4124  out << "quad ";
4125  break;
4126  default:
4127  Assert(false, ExcNotImplemented());
4128  }
4129  // note: vertex numbers are 1-base
4130  for (unsigned int vertex = 0;
4131  vertex < GeometryInfo<dim>::vertices_per_face;
4132  ++vertex)
4133  out << face->vertex_index(
4135  1
4136  << ' ';
4137  out << '\n';
4138 
4139  ++current_element_index;
4140  }
4141  return current_element_index;
4142 }
4143 
4144 
4145 
4146 template <int dim, int spacedim>
4147 unsigned int
4149  const unsigned int next_element_index,
4150  std::ostream & out) const
4151 {
4152  unsigned int current_element_index = next_element_index;
4153  // save the user flags for lines so
4154  // we can use these flags to track
4155  // which ones we've already taken
4156  // care of
4157  std::vector<bool> line_flags;
4158  const_cast<::Triangulation<dim, spacedim> &>(tria).save_user_flags_line(
4159  line_flags);
4160  const_cast<::Triangulation<dim, spacedim> &>(tria)
4161  .clear_user_flags_line();
4162 
4163  for (const auto &cell : tria.active_cell_iterators())
4164  for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
4165  if (cell->line(l)->at_boundary() && (cell->line(l)->boundary_id() != 0) &&
4166  (cell->line(l)->user_flag_set() == false))
4167  {
4168  out << current_element_index << " "
4169  << static_cast<unsigned int>(cell->line(l)->boundary_id())
4170  << " line ";
4171  // note: vertex numbers in ucd format are 1-base
4172  for (unsigned int vertex = 0; vertex < 2; ++vertex)
4173  out << cell->line(l)->vertex_index(
4175  1
4176  << ' ';
4177  out << '\n';
4178 
4179  // move on to the next line
4180  // but mark the current one
4181  // as taken care of
4182  ++current_element_index;
4183  cell->line(l)->set_user_flag();
4184  }
4185 
4186  // at the end, restore the user
4187  // flags for the lines
4188  const_cast<::Triangulation<dim, spacedim> &>(tria).load_user_flags_line(
4189  line_flags);
4190  return current_element_index;
4191 }
4192 
4193 
4194 
4195 namespace internal
4196 {
4197  namespace
4198  {
4207  template <int spacedim>
4208  void
4209  remove_colinear_points(std::vector<Point<spacedim>> &points)
4210  {
4211  while (points.size() > 2)
4212  {
4213  Tensor<1, spacedim> first_difference = points[1] - points[0];
4214  first_difference /= first_difference.norm();
4215  Tensor<1, spacedim> second_difference = points[2] - points[1];
4216  second_difference /= second_difference.norm();
4217  // If the three points are colinear then remove the middle one.
4218  if ((first_difference - second_difference).norm() < 1e-10)
4219  points.erase(points.begin() + 1);
4220  else
4221  break;
4222  }
4223  }
4224 
4225 
4226 
4227  template <int spacedim>
4228  void
4229  write_gnuplot(const ::Triangulation<1, spacedim> &tria,
4230  std::ostream & out,
4231  const Mapping<1, spacedim> *,
4233  {
4234  AssertThrow(out, ExcIO());
4235 
4236  for (const auto &cell : tria.active_cell_iterators())
4237  {
4238  if (gnuplot_flags.write_cell_numbers)
4239  out << "# cell " << cell << '\n';
4240 
4241  out << cell->vertex(0) << ' ' << cell->level() << ' '
4242  << cell->material_id() << '\n'
4243  << cell->vertex(1) << ' ' << cell->level() << ' '
4244  << cell->material_id() << '\n'
4245  << "\n\n";
4246  }
4247 
4248  // make sure everything now gets to
4249  // disk
4250  out.flush();
4251 
4252  AssertThrow(out, ExcIO());
4253  }
4254 
4255 
4256 
4257  template <int spacedim>
4258  void
4259  write_gnuplot(const ::Triangulation<2, spacedim> &tria,
4260  std::ostream & out,
4261  const Mapping<2, spacedim> * mapping,
4262  const GridOutFlags::Gnuplot & gnuplot_flags)
4263  {
4264  AssertThrow(out, ExcIO());
4265 
4266  const int dim = 2;
4267 
4268  const unsigned int n_additional_points =
4269  gnuplot_flags.n_extra_curved_line_points;
4270  const unsigned int n_points = 2 + n_additional_points;
4271 
4272  // If we need to plot curved lines then generate a quadrature formula to
4273  // place points via the mapping
4274  Quadrature<dim> q_projector;
4275  std::vector<Point<dim - 1>> boundary_points;
4276  if (mapping != nullptr)
4277  {
4278  boundary_points.resize(n_points);
4279  boundary_points[0][0] = 0;
4280  boundary_points[n_points - 1][0] = 1;
4281  for (unsigned int i = 1; i < n_points - 1; ++i)
4282  boundary_points[i](0) = 1. * i / (n_points - 1);
4283 
4284  std::vector<double> dummy_weights(n_points, 1. / n_points);
4285  Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
4286 
4288  ::ReferenceCells::Quadrilateral, quadrature);
4289  }
4290 
4291  for (const auto &cell : tria.active_cell_iterators())
4292  {
4293  if (gnuplot_flags.write_cell_numbers)
4294  out << "# cell " << cell << '\n';
4295 
4296  if (mapping == nullptr ||
4297  (dim == spacedim ?
4298  (!cell->at_boundary() && !gnuplot_flags.curved_inner_cells) :
4299  // ignore checking for boundary or interior cells in the codim
4300  // 1 case: 'or false' is a no-op
4301  false))
4302  {
4303  // write out the four sides of this cell by putting the four
4304  // points (+ the initial point again) in a row and lifting the
4305  // drawing pencil at the end
4306  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
4307  out << cell->vertex(GeometryInfo<dim>::ucd_to_deal[i]) << ' '
4308  << cell->level() << ' ' << cell->material_id() << '\n';
4309  out << cell->vertex(0) << ' ' << cell->level() << ' '
4310  << cell->material_id() << '\n'
4311  << '\n' // double new line for gnuplot 3d plots
4312  << '\n';
4313  }
4314  else
4315  // cell is at boundary and we are to treat curved boundaries. so
4316  // loop over all faces and draw them as small pieces of lines
4317  {
4318  for (const unsigned int face_no :
4320  {
4321  const typename ::Triangulation<dim,
4322  spacedim>::face_iterator
4323  face = cell->face(face_no);
4324  if (dim != spacedim || face->at_boundary() ||
4325  gnuplot_flags.curved_inner_cells)
4326  {
4327  // Save the points on each face to a vector and then try
4328  // to remove colinear points that won't show up in the
4329  // generated plot.
4330  std::vector<Point<spacedim>> line_points;
4331  // compute offset of quadrature points within set of
4332  // projected points
4333  const unsigned int offset = face_no * n_points;
4334  for (unsigned int i = 0; i < n_points; ++i)
4335  line_points.push_back(
4336  mapping->transform_unit_to_real_cell(
4337  cell, q_projector.point(offset + i)));
4338  internal::remove_colinear_points(line_points);
4339 
4340  for (const Point<spacedim> &point : line_points)
4341  out << point << ' ' << cell->level() << ' '
4342  << cell->material_id() << '\n';
4343 
4344  out << '\n' << '\n';
4345  }
4346  else
4347  {
4348  // if, however, the face is not at the boundary and we
4349  // don't want to curve anything, then draw it as usual
4350  out << face->vertex(0) << ' ' << cell->level() << ' '
4351  << cell->material_id() << '\n'
4352  << face->vertex(1) << ' ' << cell->level() << ' '
4353  << cell->material_id() << '\n'
4354  << '\n'
4355  << '\n';
4356  }
4357  }
4358  }
4359  }
4360 
4361  // make sure everything now gets to disk
4362  out.flush();
4363 
4364  AssertThrow(out, ExcIO());
4365  }
4366 
4367 
4368 
4369  template <int spacedim>
4370  void
4371  write_gnuplot(const ::Triangulation<3, spacedim> &tria,
4372  std::ostream & out,
4373  const Mapping<3, spacedim> * mapping,
4374  const GridOutFlags::Gnuplot & gnuplot_flags)
4375  {
4376  AssertThrow(out, ExcIO());
4377 
4378  const int dim = 3;
4379 
4380  const unsigned int n_additional_points =
4381  gnuplot_flags.n_extra_curved_line_points;
4382  const unsigned int n_points = 2 + n_additional_points;
4383 
4384  // If we need to plot curved lines then generate a quadrature formula to
4385  // place points via the mapping
4386  std::unique_ptr<Quadrature<dim>> q_projector;
4387  std::vector<Point<1>> boundary_points;
4388  if (mapping != nullptr)
4389  {
4390  boundary_points.resize(n_points);
4391  boundary_points[0][0] = 0;
4392  boundary_points[n_points - 1][0] = 1;
4393  for (unsigned int i = 1; i < n_points - 1; ++i)
4394  boundary_points[i](0) = 1. * i / (n_points - 1);
4395 
4396  std::vector<double> dummy_weights(n_points, 1. / n_points);
4397  Quadrature<1> quadrature1d(boundary_points, dummy_weights);
4398 
4399  // tensor product of points, only one copy
4400  QIterated<dim - 1> quadrature(quadrature1d, 1);
4401  q_projector = std::make_unique<Quadrature<dim>>(
4403  }
4404 
4405  for (const auto &cell : tria.active_cell_iterators())
4406  {
4407  if (gnuplot_flags.write_cell_numbers)
4408  out << "# cell " << cell << '\n';
4409 
4410  if (mapping == nullptr || n_points == 2 ||
4411  (!cell->has_boundary_lines() &&
4412  !gnuplot_flags.curved_inner_cells))
4413  {
4414  // front face
4415  out << cell->vertex(0) << ' ' << cell->level() << ' '
4416  << cell->material_id() << '\n'
4417  << cell->vertex(1) << ' ' << cell->level() << ' '
4418  << cell->material_id() << '\n'
4419  << cell->vertex(5) << ' ' << cell->level() << ' '
4420  << cell->material_id() << '\n'
4421  << cell->vertex(4) << ' ' << cell->level() << ' '
4422  << cell->material_id() << '\n'
4423  << cell->vertex(0) << ' ' << cell->level() << ' '
4424  << cell->material_id() << '\n'
4425  << '\n';
4426  // back face
4427  out << cell->vertex(2) << ' ' << cell->level() << ' '
4428  << cell->material_id() << '\n'
4429  << cell->vertex(3) << ' ' << cell->level() << ' '
4430  << cell->material_id() << '\n'
4431  << cell->vertex(7) << ' ' << cell->level() << ' '
4432  << cell->material_id() << '\n'
4433  << cell->vertex(6) << ' ' << cell->level() << ' '
4434  << cell->material_id() << '\n'
4435  << cell->vertex(2) << ' ' << cell->level() << ' '
4436  << cell->material_id() << '\n'
4437  << '\n';
4438 
4439  // now for the four connecting lines
4440  out << cell->vertex(0) << ' ' << cell->level() << ' '
4441  << cell->material_id() << '\n'
4442  << cell->vertex(2) << ' ' << cell->level() << ' '
4443  << cell->material_id() << '\n'
4444  << '\n';
4445  out << cell->vertex(1) << ' ' << cell->level() << ' '
4446  << cell->material_id() << '\n'
4447  << cell->vertex(3) << ' ' << cell->level() << ' '
4448  << cell->material_id() << '\n'
4449  << '\n';
4450  out << cell->vertex(5) << ' ' << cell->level() << ' '
4451  << cell->material_id() << '\n'
4452  << cell->vertex(7) << ' ' << cell->level() << ' '
4453  << cell->material_id() << '\n'
4454  << '\n';
4455  out << cell->vertex(4) << ' ' << cell->level() << ' '
4456  << cell->material_id() << '\n'
4457  << cell->vertex(6) << ' ' << cell->level() << ' '
4458  << cell->material_id() << '\n'
4459  << '\n';
4460  }
4461  else
4462  {
4463  for (const unsigned int face_no :
4465  {
4466  const typename ::Triangulation<dim,
4467  spacedim>::face_iterator
4468  face = cell->face(face_no);
4469 
4470  if (face->at_boundary() &&
4471  gnuplot_flags.write_additional_boundary_lines)
4472  {
4473  const unsigned int offset = face_no * n_points * n_points;
4474  for (unsigned int i = 0; i < n_points - 1; ++i)
4475  for (unsigned int j = 0; j < n_points - 1; ++j)
4476  {
4477  const Point<spacedim> p0 =
4478  mapping->transform_unit_to_real_cell(
4479  cell,
4480  q_projector->point(offset + i * n_points + j));
4481  out << p0 << ' ' << cell->level() << ' '
4482  << cell->material_id() << '\n';
4483  out << (mapping->transform_unit_to_real_cell(
4484  cell,
4485  q_projector->point(
4486  offset + (i + 1) * n_points + j)))
4487  << ' ' << cell->level() << ' '
4488  << cell->material_id() << '\n';
4489  out << (mapping->transform_unit_to_real_cell(
4490  cell,
4491  q_projector->point(
4492  offset + (i + 1) * n_points + j + 1)))
4493  << ' ' << cell->level() << ' '
4494  << cell->material_id() << '\n';
4495  out << (mapping->transform_unit_to_real_cell(
4496  cell,
4497  q_projector->point(offset + i * n_points +
4498  j + 1)))
4499  << ' ' << cell->level() << ' '
4500  << cell->material_id() << '\n';
4501  // and the first point again
4502  out << p0 << ' ' << cell->level() << ' '
4503  << cell->material_id() << '\n';
4504  out << '\n' << '\n';
4505  }
4506  }
4507  else
4508  {
4509  for (unsigned int l = 0;
4510  l < GeometryInfo<dim>::lines_per_face;
4511  ++l)
4512  {
4513  const typename ::Triangulation<dim, spacedim>::
4514  line_iterator line = face->line(l);
4515 
4516  const Point<spacedim> &v0 = line->vertex(0),
4517  &v1 = line->vertex(1);
4518  if (line->at_boundary() ||
4519  gnuplot_flags.curved_inner_cells)
4520  {
4521  // Save the points on each face to a vector and
4522  // then try to remove colinear points that won't
4523  // show up in the generated plot.
4524  std::vector<Point<spacedim>> line_points;
4525  // transform_real_to_unit_cell could be replaced
4526  // by using QProjector<dim>::project_to_line
4527  // which is not yet implemented
4528  const Point<spacedim>
4529  u0 = mapping->transform_real_to_unit_cell(cell,
4530  v0),
4531  u1 = mapping->transform_real_to_unit_cell(cell,
4532  v1);
4533  for (unsigned int i = 0; i < n_points; ++i)
4534  line_points.push_back(
4535  mapping->transform_unit_to_real_cell(
4536  cell,
4537  (1 - boundary_points[i][0]) * u0 +
4538  boundary_points[i][0] * u1));
4539  internal::remove_colinear_points(line_points);
4540  for (const Point<spacedim> &point : line_points)
4541  out << point << ' ' << cell->level() << ' '
4542  << static_cast<unsigned int>(
4543  cell->material_id())
4544  << '\n';
4545  }
4546  else
4547  out << v0 << ' ' << cell->level() << ' '
4548  << cell->material_id() << '\n'
4549  << v1 << ' ' << cell->level() << ' '
4550  << cell->material_id() << '\n';
4551 
4552  out << '\n' << '\n';
4553  }
4554  }
4555  }
4556  }
4557  }
4558 
4559  // make sure everything now gets to disk
4560  out.flush();
4561 
4562  AssertThrow(out, ExcIO());
4563  }
4564  } // namespace
4565 } // namespace internal
4566 
4567 
4568 
4569 template <int dim, int spacedim>
4570 void
4572  std::ostream & out,
4573  const Mapping<dim, spacedim> * mapping) const
4574 {
4575  internal::write_gnuplot(tria, out, mapping, gnuplot_flags);
4576 }
4577 
4578 
4579 
4580 namespace internal
4581 {
4582  namespace
4583  {
4584  struct LineEntry
4585  {
4588  bool colorize;
4589  unsigned int level;
4590  LineEntry(const Point<2> & f,
4591  const Point<2> & s,
4592  const bool c,
4593  const unsigned int l)
4594  : first(f)
4595  , second(s)
4596  , colorize(c)
4597  , level(l)
4598  {}
4599  };
4600 
4601 
4602  void
4603  write_eps(const ::Triangulation<1> &,
4604  std::ostream &,
4605  const Mapping<1> *,
4606  const GridOutFlags::Eps<2> &,
4607  const GridOutFlags::Eps<3> &)
4608  {
4609  Assert(false, ExcNotImplemented());
4610  }
4611 
4612  void
4613  write_eps(const ::Triangulation<1, 2> &,
4614  std::ostream &,
4615  const Mapping<1, 2> *,
4616  const GridOutFlags::Eps<2> &,
4617  const GridOutFlags::Eps<3> &)
4618  {
4619  Assert(false, ExcNotImplemented());
4620  }
4621 
4622  void
4623  write_eps(const ::Triangulation<1, 3> &,
4624  std::ostream &,
4625  const Mapping<1, 3> *,
4626  const GridOutFlags::Eps<2> &,
4627  const GridOutFlags::Eps<3> &)
4628  {
4629  Assert(false, ExcNotImplemented());
4630  }
4631 
4632  void
4633  write_eps(const ::Triangulation<2, 3> &,
4634  std::ostream &,
4635  const Mapping<2, 3> *,
4636  const GridOutFlags::Eps<2> &,
4637  const GridOutFlags::Eps<3> &)
4638  {
4639  Assert(false, ExcNotImplemented());
4640  }
4641 
4642 
4643 
4644  template <int dim, int spacedim>
4645  void
4646  write_eps(const ::Triangulation<dim, spacedim> &tria,
4647  std::ostream & out,
4648  const Mapping<dim, spacedim> * mapping,
4651  {
4652  using LineList = std::list<LineEntry>;
4653 
4654  // We should never get here in 1D since this function is overloaded for
4655  // all dim == 1 cases.
4656  Assert(dim == 2 || dim == 3, ExcInternalError());
4657 
4658  // Copy, with an object slice, something containing the flags common to
4659  // all dimensions in order to avoid the recurring distinctions between
4660  // the different eps_flags present.
4661  const GridOutFlags::EpsFlagsBase eps_flags_base =
4662  dim == 2 ?
4663  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_2) :
4664  static_cast<const GridOutFlags::EpsFlagsBase &>(eps_flags_3);
4665 
4666  AssertThrow(out, ExcIO());
4667  const unsigned int n_points = eps_flags_base.n_boundary_face_points;
4668 
4669  // make up a list of lines by which
4670  // we will construct the triangulation
4671  //
4672  // this part unfortunately is a bit
4673  // dimension dependent, so we have to
4674  // treat every dimension different.
4675  // however, by directly producing
4676  // the lines to be printed, i.e. their
4677  // 2d images, we can later do the
4678  // actual output dimension independent
4679  // again
4680  LineList line_list;
4681 
4682  switch (dim)
4683  {
4684  case 1:
4685  {
4686  Assert(false, ExcInternalError());
4687  break;
4688  }
4689 
4690  case 2:
4691  {
4692  for (const auto &cell : tria.active_cell_iterators())
4693  for (const unsigned int line_no : cell->line_indices())
4694  {
4695  typename ::Triangulation<dim, spacedim>::line_iterator
4696  line = cell->line(line_no);
4697 
4698  // first treat all
4699  // interior lines and
4700  // make up a list of
4701  // them. if curved
4702  // lines shall not be
4703  // supported (i.e. no
4704  // mapping is
4705  // provided), then also
4706  // treat all other
4707  // lines
4708  if (!line->has_children() &&
4709  (mapping == nullptr || !line->at_boundary()))
4710  // one would expect
4711  // make_pair(line->vertex(0),
4712  // line->vertex(1))
4713  // here, but that is
4714  // not dimension
4715  // independent, since
4716  // vertex(i) is
4717  // Point<dim>, but we
4718  // want a Point<2>.
4719  // in fact, whenever
4720  // we're here, the
4721  // vertex is a
4722  // Point<dim>, but
4723  // the compiler does
4724  // not know
4725  // this. hopefully,
4726  // the compiler will
4727  // optimize away this
4728  // little kludge
4729  line_list.emplace_back(
4730  Point<2>(line->vertex(0)(0), line->vertex(0)(1)),
4731  Point<2>(line->vertex(1)(0), line->vertex(1)(1)),
4732  line->user_flag_set(),
4733  cell->level());
4734  }
4735 
4736  // next if we are to treat
4737  // curved boundaries
4738  // specially, then add lines
4739  // to the list consisting of
4740  // pieces of the boundary
4741  // lines
4742  if (mapping != nullptr)
4743  {
4744  // to do so, first
4745  // generate a sequence of
4746  // points on a face and
4747  // project them onto the
4748  // faces of a unit cell
4749  std::vector<Point<dim - 1>> boundary_points(n_points);
4750 
4751  for (unsigned int i = 0; i < n_points; ++i)
4752  boundary_points[i](0) = 1. * (i + 1) / (n_points + 1);
4753 
4754  Quadrature<dim - 1> quadrature(boundary_points);
4755  Quadrature<dim> q_projector(
4757 
4758  // next loop over all
4759  // boundary faces and
4760  // generate the info from
4761  // them
4762  for (const auto &cell : tria.active_cell_iterators())
4763  for (const unsigned int face_no :
4765  {
4766  const typename ::Triangulation<dim, spacedim>::
4767  face_iterator face = cell->face(face_no);
4768 
4769  if (face->at_boundary())
4770  {
4771  Point<dim> p0_dim(face->vertex(0));
4772  Point<2> p0(p0_dim(0), p0_dim(1));
4773 
4774  // loop over
4775  // all pieces
4776  // of the line
4777  // and generate
4778  // line-lets
4779  const unsigned int offset = face_no * n_points;
4780  for (unsigned int i = 0; i < n_points; ++i)
4781  {
4782  const Point<dim> p1_dim(
4783  mapping->transform_unit_to_real_cell(
4784  cell, q_projector.point(offset + i)));
4785  const Point<2> p1(p1_dim(0), p1_dim(1));
4786 
4787  line_list.emplace_back(p0,
4788  p1,
4789  face->user_flag_set(),
4790  cell->level());
4791  p0 = p1;
4792  }
4793 
4794  // generate last piece
4795  const Point<dim> p1_dim(face->vertex(1));
4796  const Point<2> p1(p1_dim(0), p1_dim(1));
4797  line_list.emplace_back(p0,
4798  p1,
4799  face->user_flag_set(),
4800  cell->level());
4801  }
4802  }
4803  }
4804 
4805  break;
4806  }
4807 
4808  case 3:
4809  {
4810  // curved boundary output
4811  // presently not supported
4812  Assert(mapping == nullptr, ExcNotImplemented());
4813 
4814  // loop over all lines and compute their
4815  // projection on the plane perpendicular
4816  // to the direction of sight
4817 
4818  // direction of view equals the unit
4819  // vector of the position of the
4820  // spectator to the origin.
4821  //
4822  // we chose here the viewpoint as in
4823  // gnuplot as default.
4824  //
4825  // TODO:[WB] Fix a potential problem with viewing angles in 3d Eps
4826  // GridOut
4827  // note: the following might be wrong
4828  // if one of the base vectors below
4829  // is in direction of the viewer, but
4830  // I am too tired at present to fix
4831  // this
4832  const double pi = numbers::PI;
4833  const double z_angle = eps_flags_3.azimut_angle;
4834  const double turn_angle = eps_flags_3.turn_angle;
4835  const Point<dim> view_direction(
4836  -std::sin(z_angle * 2. * pi / 360.) *
4837  std::sin(turn_angle * 2. * pi / 360.),
4838  +std::sin(z_angle * 2. * pi / 360.) *
4839  std::cos(turn_angle * 2. * pi / 360.),
4840  -std::cos(z_angle * 2. * pi / 360.));
4841 
4842  // decide about the two unit vectors
4843  // in this plane. we chose the first one
4844  // to be the projection of the z-axis
4845  // to this plane
4846  const Tensor<1, dim> vector1 =
4847  Point<dim>(0, 0, 1) -
4848  ((Point<dim>(0, 0, 1) * view_direction) * view_direction);
4849  const Tensor<1, dim> unit_vector1 = vector1 / vector1.norm();
4850 
4851  // now the third vector is fixed. we
4852  // chose the projection of a more or
4853  // less arbitrary vector to the plane
4854  // perpendicular to the first one
4855  const Tensor<1, dim> vector2 =
4856  (Point<dim>(1, 0, 0) -
4857  ((Point<dim>(1, 0, 0) * view_direction) * view_direction) -
4858  ((Point<dim>(1, 0, 0) * unit_vector1) * unit_vector1));
4859  const Tensor<1, dim> unit_vector2 = vector2 / vector2.norm();
4860 
4861 
4862  for (const auto &cell : tria.active_cell_iterators())
4863  for (const unsigned int line_no : cell->line_indices())
4864  {
4865  typename ::Triangulation<dim, spacedim>::line_iterator
4866  line = cell->line(line_no);
4867  line_list.emplace_back(
4868  Point<2>(line->vertex(0) * unit_vector2,
4869  line->vertex(0) * unit_vector1),
4870  Point<2>(line->vertex(1) * unit_vector2,
4871  line->vertex(1) * unit_vector1),
4872  line->user_flag_set(),
4873  cell->level());
4874  }
4875 
4876  break;
4877  }
4878 
4879  default:
4880  Assert(false, ExcNotImplemented());
4881  }
4882 
4883 
4884 
4885  // find out minimum and maximum x and
4886  // y coordinates to compute offsets
4887  // and scaling factors
4888  double x_min = tria.begin_active()->vertex(0)(0);
4889  double x_max = x_min;
4890  double y_min = tria.begin_active()->vertex(0)(1);
4891  double y_max = y_min;
4892  unsigned int max_level = line_list.begin()->level;
4893 
4894  for (LineList::const_iterator line = line_list.begin();
4895  line != line_list.end();
4896  ++line)
4897  {
4898  x_min = std::min(x_min, line->first(0));
4899  x_min = std::min(x_min, line->second(0));
4900 
4901  x_max = std::max(x_max, line->first(0));
4902  x_max = std::max(x_max, line->second(0));
4903 
4904  y_min = std::min(y_min, line->first(1));
4905  y_min = std::min(y_min, line->second(1));
4906 
4907  y_max = std::max(y_max, line->first(1));
4908  y_max = std::max(y_max, line->second(1));
4909 
4910  max_level = std::max(max_level, line->level);
4911  }
4912 
4913  // scale in x-direction such that
4914  // in the output 0 <= x <= 300.
4915  // don't scale in y-direction to
4916  // preserve the shape of the
4917  // triangulation
4918  const double scale =
4919  (eps_flags_base.size /
4920  (eps_flags_base.size_type == GridOutFlags::EpsFlagsBase::width ?
4921  x_max - x_min :
4922  y_min - y_max));
4923 
4924 
4925  // now write preamble
4926  {
4927  // block this to have local
4928  // variables destroyed after
4929  // use
4930  std::time_t time1 = std::time(nullptr);
4931  std::tm * time = std::localtime(&time1);
4932  out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
4933  << "%%Title: deal.II Output" << '\n'
4934  << "%%Creator: the deal.II library" << '\n'
4935  << "%%Creation Date: " << time->tm_year + 1900 << "/"
4936  << time->tm_mon + 1 << "/" << time->tm_mday << " - "
4937  << time->tm_hour << ":" << std::setw(2) << time->tm_min << ":"
4938  << std::setw(2) << time->tm_sec << '\n'
4939  << "%%BoundingBox: "
4940  // lower left corner
4941  << "0 0 "
4942  // upper right corner
4943  << static_cast<unsigned int>(
4944  std::floor(((x_max - x_min) * scale) + 1))
4945  << ' '
4946  << static_cast<unsigned int>(
4947  std::floor(((y_max - y_min) * scale) + 1))
4948  << '\n';
4949 
4950  // define some abbreviations to keep
4951  // the output small:
4952  // m=move turtle to
4953  // x=execute line stroke
4954  // b=black pen
4955  // r=red pen
4956  out << "/m {moveto} bind def" << '\n'
4957  << "/x {lineto stroke} bind def" << '\n'
4958  << "/b {0 0 0 setrgbcolor} def" << '\n'
4959  << "/r {1 0 0 setrgbcolor} def" << '\n';
4960 
4961  // calculate colors for level
4962  // coloring; level 0 is black,
4963  // other levels are blue
4964  // ... red
4965  if (eps_flags_base.color_lines_level)
4966  out << "/l { neg " << (max_level) << " add "
4967  << (0.66666 / std::max(1U, (max_level - 1)))
4968  << " mul 1 0.8 sethsbcolor} def" << '\n';
4969 
4970  // in 2d, we can also plot cell
4971  // and vertex numbers, but this
4972  // requires a somewhat more
4973  // lengthy preamble. please
4974  // don't ask me what most of
4975  // this means, it is reverse
4976  // engineered from what GNUPLOT
4977  // uses in its output
4978  if ((dim == 2) && (eps_flags_2.write_cell_numbers ||
4979  eps_flags_2.write_vertex_numbers))
4980  {
4981  out
4982  << ("/R {rmoveto} bind def\n"
4983  "/Symbol-Oblique /Symbol findfont [1 0 .167 1 0 0] makefont\n"
4984  "dup length dict begin {1 index /FID eq {pop pop} {def} ifelse} forall\n"
4985  "currentdict end definefont\n"
4986  "/MFshow {{dup dup 0 get findfont exch 1 get scalefont setfont\n"
4987  "[ currentpoint ] exch dup 2 get 0 exch rmoveto dup dup 5 get exch 4 get\n"
4988  "{show} {stringwidth pop 0 rmoveto}ifelse dup 3 get\n"
4989  "{2 get neg 0 exch rmoveto pop} {pop aload pop moveto}ifelse} forall} bind def\n"
4990  "/MFwidth {0 exch {dup 3 get{dup dup 0 get findfont exch 1 get scalefont setfont\n"
4991  "5 get stringwidth pop add}\n"
4992  "{pop} ifelse} forall} bind def\n"
4993  "/MCshow { currentpoint stroke m\n"
4994  "exch dup MFwidth -2 div 3 -1 roll R MFshow } def\n")
4995  << '\n';
4996  }
4997 
4998  out << "%%EndProlog" << '\n' << '\n';
4999 
5000  // set fine lines
5001  out << eps_flags_base.line_width << " setlinewidth" << '\n';
5002  }
5003 
5004  // now write the lines
5005  const Point<2> offset(x_min, y_min);
5006 
5007  for (LineList::const_iterator line = line_list.begin();
5008  line != line_list.end();
5009  ++line)
5010  if (eps_flags_base.color_lines_level && (line->level > 0))
5011  out << line->level << " l " << (line->first - offset) * scale << " m "
5012  << (line->second - offset) * scale << " x" << '\n';
5013  else
5014  out << ((line->colorize && eps_flags_base.color_lines_on_user_flag) ?
5015  "r " :
5016  "b ")
5017  << (line->first - offset) * scale << " m "
5018  << (line->second - offset) * scale << " x" << '\n';
5019 
5020  // finally write the cell numbers
5021  // in 2d, if that is desired
5022  if ((dim == 2) && (eps_flags_2.write_cell_numbers == true))
5023  {
5024  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5025 
5026  for (const auto &cell : tria.active_cell_iterators())
5027  {
5028  out << (cell->center()(0) - offset(0)) * scale << ' '
5029  << (cell->center()(1) - offset(1)) * scale << " m" << '\n'
5030  << "[ [(Helvetica) 12.0 0.0 true true (";
5031  if (eps_flags_2.write_cell_number_level)
5032  out << cell;
5033  else
5034  out << cell->index();
5035 
5036  out << ")] "
5037  << "] -6 MCshow" << '\n';
5038  }
5039  }
5040 
5041  // and the vertex numbers
5042  if ((dim == 2) && (eps_flags_2.write_vertex_numbers == true))
5043  {
5044  out << "(Helvetica) findfont 140 scalefont setfont" << '\n';
5045 
5046  // have a list of those
5047  // vertices which we have
5048  // already tracked, to avoid
5049  // doing this multiply
5050  std::set<unsigned int> treated_vertices;
5051  for (const auto &cell : tria.active_cell_iterators())
5052  for (const unsigned int vertex_no : cell->vertex_indices())
5053  if (treated_vertices.find(cell->vertex_index(vertex_no)) ==
5054  treated_vertices.end())
5055  {
5056  treated_vertices.insert(cell->vertex_index(vertex_no));
5057 
5058  out << (cell->vertex(vertex_no)(0) - offset(0)) * scale << ' '
5059  << (cell->vertex(vertex_no)(1) - offset(1)) * scale
5060  << " m" << '\n'
5061  << "[ [(Helvetica) 10.0 0.0 true true ("
5062  << cell->vertex_index(vertex_no) << ")] "
5063  << "] -6 MCshow" << '\n';
5064  }
5065  }
5066 
5067  out << "showpage" << '\n';
5068 
5069  // make sure everything now gets to
5070  // disk
5071  out.flush();
5072 
5073  AssertThrow(out, ExcIO());
5074  }
5075  } // namespace
5076 } // namespace internal
5077 
5078 
5079 template <int dim, int spacedim>
5080 void
5082  std::ostream & out,
5083  const Mapping<dim, spacedim> * mapping) const
5084 {
5085  internal::write_eps(tria, out, mapping, eps_flags_2, eps_flags_3);
5086 }
5087 
5088 
5089 template <int dim, int spacedim>
5090 void
5092  std::ostream & out,
5093  const OutputFormat output_format,
5094  const Mapping<dim, spacedim> * mapping) const
5095 {
5096  switch (output_format)
5097  {
5098  case none:
5099  return;
5100 
5101  case dx:
5102  write_dx(tria, out);
5103  return;
5104 
5105  case ucd:
5106  write_ucd(tria, out);
5107  return;
5108 
5109  case gnuplot:
5110  write_gnuplot(tria, out, mapping);
5111  return;
5112 
5113  case eps:
5114  write_eps(tria, out, mapping);
5115  return;
5116 
5117  case xfig:
5118  write_xfig(tria, out, mapping);
5119  return;
5120 
5121  case msh:
5122  write_msh(tria, out);
5123  return;
5124 
5125  case svg:
5126  write_svg(tria, out);
5127  return;
5128 
5129  case mathgl:
5130  write_mathgl(tria, out);
5131  return;
5132 
5133  case vtk:
5134  write_vtk(tria, out);
5135  return;
5136 
5137  case vtu:
5138  write_vtu(tria, out);
5139  return;
5140  }
5141 
5142  Assert(false, ExcInternalError());
5143 }
5144 
5145 
5146 template <int dim, int spacedim>
5147 void
5149  std::ostream & out,
5150  const Mapping<dim, spacedim> * mapping) const
5151 {
5152  write(tria, out, default_format, mapping);
5153 }
5154 
5155 
5156 // explicit instantiations
5157 #include "grid_out.inst"
5158 
5159 
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:3825
unsigned int n_active_cells() const
Definition: tria.cc:13754
long int get_integer(const std::string &entry_string) const
const std::vector< ReferenceCell > & get_reference_cells() const
Definition: tria.cc:14618
unsigned int n_used_vertices() const
Definition: tria.cc:14294
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:1655
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
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:1275
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:2047
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:13264
OutputFormat
Definition: grid_out.h:1000
const ::Triangulation< dim, spacedim > & tria
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3590
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:3637
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:963
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:13067
static const char U
unsigned int line_thickness
Definition: grid_out.h:672
const unsigned int v1
Definition: grid_tools.cc:963
bool convert_level_number_to_height
Definition: grid_out.h:744
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
bool write_diameter
Definition: grid_out.h:71
Point< 2 > second
Definition: grid_out.cc:4587
void write_mathgl(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:3021
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:13047
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:2650
unsigned int write_msh_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:3945
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:4571
cell_iterator end() const
Definition: tria.cc:13158
unsigned int n_boundary_lines(const Triangulation< dim, spacedim > &tria) const
Definition: grid_out.cc:3841
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:1738
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:13275
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:12358
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:1461
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:401
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:651
unsigned int level
Definition: grid_out.cc:4589
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:3351
IteratorRange< active_face_iterator > active_face_iterators() const
Definition: tria.cc:13363
GridOutFlags::Ucd ucd_flags
Definition: grid_out.h:1566
static void declare_parameters(ParameterHandler &param)
Definition: grid_out.cc:136
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:185
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:13697
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Point< 2 > first
Definition: grid_out.cc:4586
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:506
unsigned int write_ucd_faces(const Triangulation< dim, spacedim > &tria, const unsigned int next_element_index, std::ostream &out) const
Definition: grid_out.cc:4106
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:14303
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
void write_ucd(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition: grid_out.cc:1158
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:4148
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:3987
GridOutFlags::Svg svg_flags
Definition: grid_out.h:1600
void attach_triangulation(const Triangulation< dim, spacedim > &)
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:5081
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:4588
IteratorRange< cell_iterator > cell_iterators() const
Definition: tria.cc:13255
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:5091
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:14370
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:1019
Expression floor(const Expression &x)
std::string default_suffix() const
Definition: grid_out.cc:593