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