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