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